Moknapsack2 (#27)

* raster.write_raster fixes: nodata, geotransform
   - gdal Only north-up non rotated geotransform supported

* knapsack: plot fixes, config options, prints -> logs
    - knapsack aplot args show=False
    - unit testing init

---------

Co-authored-by: fdobad <[email protected]>

Author: fdobad

src/fire2a/knapsack.py

@@ -408,8 +408,8 @@ def write_raster(
     outfile="output.tif",
     driver_name="GTiff",
     authid="EPSG:3857",
-    geotransform=(0, 1, 0, 0, 0, 1),
-    nodata=None,
+    geotransform=(0, 1, 0, 0, 0, -1),
+    nodata: int | None = None,
     feedback=None,
     logger=None,  # logger default ?
 ):
@@ -450,13 +450,15 @@ def write_raster(
     ds.SetGeoTransform(geotransform)
     ds.SetProjection(authid)
     band = ds.GetRasterBand(1)
-    if nodata:
-        if 0 != band.SetNoDataValue(nodata):
-            fprint("Set NoData failed", level="warning", feedback=feedback, logger=logger)
-            return False
     if 0 != band.WriteArray(data):
-        fprint(f"WriteArray failed for Burn Probability {burn_prob}", level="warning", feedback=feedback, logger=logger)
+        fprint("WriteArray failed", level="warning", feedback=feedback, logger=logger)
         return False
+    if nodata and data[data == nodata].size > 0:
+        band.SetNoDataValue(nodata)
+        # TBD : always returns 1?
+        # if 0 != band.SetNoDataValue(nodata):
+        #     fprint("Set NoData failed", level="warning", feedback=feedback, logger=logger)
+        #     return False
     ds.FlushCache()
     ds = None
     return True

src/fire2a/raster.py

@@ -0,0 +1,19 @@
+["/home/fdo/source/fire/mok-tmp/v0.tif"]
+value_weight = 1.0
+
+["/home/fdo/source/fire/mok-tmp/v1.tif"]
+value_weight = 1.25
+
+["/home/fdo/source/fire/mok-tmp/v2.tif"]
+value_weight = 1.25
+
+["/home/fdo/source/fire/mok-tmp/v3.tif"]
+value_weight = 1.0
+
+["/home/fdo/source/fire/mok-tmp/w0.tif"]
+capacity_sense = "<="
+capacity_ratio = 0.25
+
+["/home/fdo/source/fire/mok-tmp/w1.tif"]
+capacity_sense = "<="
+capacity_ratio = 0.25

tests/mok/OUT_LAYER.gpkg

@@ -0,0 +1,77 @@
+#!python3
+"""
+test instance for multi objective knapsack problem
+"""
+import sys
+
+import numpy as np
+from fire2a.raster import write_raster
+from matplotlib import pyplot as plt
+
+# arguments
+show, save, write = False, False, False
+if "show" in sys.argv:
+    show = True
+if "save" in sys.argv:
+    save = True
+if "write" in sys.argv:
+    write = True
+
+width = 80
+height = 60
+
+x, y = np.meshgrid(np.arange(width), np.arange(height))
+
+v = np.zeros((height, width, 4))
+v[:, :, 0] = np.sin(x * 2 * np.pi / width)
+v[:, :, 1] = -np.sin(x * 2 * np.pi / width)
+v[:, :, 2] = np.sin(y * 2 * np.pi / height)
+v[:, :, 3] = -np.sin(y * 2 * np.pi / height)
+
+if show or save:
+    fig, ax = plt.subplots(2, 2)
+    fig.suptitle("values")
+    for i, (j, k) in enumerate(np.indices((2, 2)).reshape(2, -1).T):
+        ax[j, k].imshow(v[:, :, i])
+        ax[j, k].set_title("v" + str(i))
+    if show:
+        plt.show()
+    if save:
+        plt.savefig("v.png")
+    plt.close()
+
+# all layers nodata
+v[1, 1, :] = -9999
+# one layer nodata
+v[0, 0, 0] = -9999
+
+if write:
+    for i in range(4):
+        write_raster(v[:, :, i], "v" + str(i) + ".tif", nodata=-9999)
+
+w = np.ones((height, width, 2))
+w[:, :, 0] += 1
+w[:, :, 1] += np.triu(w[:, :, 1])
+w[:, :, 1] += np.triu(w[:, :, 1], width // 3)
+w[:, :, 1] += np.triu(w[:, :, 1], width * 2 // 3)
+
+if show or save:
+    fig, ax = plt.subplots(1, 2)
+    fig.suptitle("weights")
+    for i in range(2):
+        ax[i].imshow(w[:, :, i])
+        ax[i].set_title("w" + str(i))
+    if show:
+        plt.show()
+    if save:
+        plt.savefig("w.png")
+    plt.close()
+
+# all layers nodata
+w[1, 1, :] = -9999
+# one layer nodata
+w[2, 2, 0] = -9999
+
+if write:
+    for i in range(2):
+        write_raster(w[:, :, i], "w" + str(i) + ".tif", nodata=-9999)

tests/mok/config.toml

@@ -0,0 +1 @@
+<span style="color:#777">QGIS version: 3.40.2-Bratislava</span><br/><span style="color:#777">QGIS code revision: 14826ca1e4a</span><br/><span style="color:#777">Qt version: 5.15.8</span><br/><span style="color:#777">Python version: 3.11.2</span><br/><span style="color:#777">GDAL version: 3.6.2</span><br/><span style="color:#777">GEOS version: 3.11.1-CAPI-1.17.1</span><br/><span style="color:#777">PROJ version: Rel. 9.1.1, December 1st, 2022</span><br/>Algorithm started at: 2025-01-06T00:39:16<br/>Input parameters:<br/><code>{ 'CUSTOM_OPTIONS_STRING' : '', 'DISPLAY_MODEL' : False, 'EXECUTABLE' : '', 'InputRasters' : ['/home/fdo/source/fire/mok-tmp/v0.tif','/home/fdo/source/fire/mok-tmp/v1.tif','/home/fdo/source/fire/mok-tmp/v2.tif','/home/fdo/source/fire/mok-tmp/v3.tif','/home/fdo/source/fire/mok-tmp/w0.tif','/home/fdo/source/fire/mok-tmp/w1.tif'], 'Matrix' : ['','1','','','','1.25','','','','1.25','','','','1','','','','','&lt;=','0.25','','','&lt;=','0.25'], 'NEOS_CUSTOM_OPTIONS_STRING' : '', 'NEOS_EMAIL' : '', 'NEOS_SOLVER' : 'cplex', 'OUT_LAYER' : 'TEMPORARY_OUTPUT', 'PLOTS' : True, 'SOLVER' : 'cplex: mipgap=0.005 timelimit=300' }</code><br/><br/>Solver unavailability:<br><br><br/><span style="color:#777">Input rasters names: ['v0', 'v1', 'v2', 'v3', 'w0', 'w1']<br></span><br/><span style="color:#777">/home/fdo/source/fire/mok-tmp/v0.tif : {'value_weight': 1.0}<br></span><br/><span style="color:#777">/home/fdo/source/fire/mok-tmp/v1.tif : {'value_weight': 1.25}<br></span><br/><span style="color:#777">/home/fdo/source/fire/mok-tmp/v2.tif : {'value_weight': 1.25}<br></span><br/><span style="color:#777">/home/fdo/source/fire/mok-tmp/v3.tif : {'value_weight': 1.0}<br></span><br/><span style="color:#777">/home/fdo/source/fire/mok-tmp/w0.tif : {'capacity_sense': '&lt;=', 'capacity_ratio': 0.25}<br></span><br/><span style="color:#777">/home/fdo/source/fire/mok-tmp/w1.tif : {'capacity_sense': '&lt;=', 'capacity_ratio': 0.25}<br></span><br/><span style="color:#b85a20;">fire2a-lib cli alternative:<br>python -m fire2a.knapsack -vv --authid EPSG:3857 --geotransform (0.0, 1.0, 0.0, 0.0, 0.0, -1.0) --plots --output_raster /tmp/processing_MoaEFQ/0ef1ab26b5dc4c0c9c463ca45fe8752c/OUT_LAYER.gpkg /tmp/processing_MoaEFQ/0ef1ab26b5dc4c0c9c463ca45fe8752c/tmptlh3t589.toml<br>Depending on the terminal the geotransform might need quotes &quot;(0,0,1,0,1)&quot; around it to be read correctly</span><br/><br/><span style="color:#777">options_string: mipgap=0.005 timelimit=300<br></span><br/>Solver object created 33%<br/><code style="color:#777"></code><br/><code style="color:#777"><br>Welcome to IBM(R) ILOG(R) CPLEX(R) Interactive Optimizer 22.1.1.0<br>  with Simplex, Mixed Integer &amp; Barrier Optimizers<br>5725-A06 5725-A29 5724-Y48 5724-Y49 5724-Y54 5724-Y55 5655-Y21<br>Copyright IBM Corp. 1988, 2022.  All Rights Reserved.<br>Type 'help' for a list of available commands.<br>Type 'help' followed by a command name for more<br>information on commands.<br></code><br/><code style="color:#777"><br>CPLEX&gt; Logfile 'cplex.log' closed.<br>Logfile '/tmp/tmp2799_93f.cplex.log' open.<br></code><br/><code style="color:#777">CPLEX&gt; New value for mixed integer optimality gap tolerance: 0.005<br></code><br/><code style="color:#777">CPLEX&gt; New value for time limit in seconds: 300<br></code><br/><code style="color:#777">CPLEX&gt; Problem '/tmp/tmpng6qoa01.pyomo.lp' read.<br>Read time = 0.01 sec. (0.26 ticks)<br></code><br/><code style="color:#777">CPLEX&gt; Problem name         : /tmp/tmpng6qoa01.pyomo.lp<br>Objective sense      : Maximize<br>Variables            :    4799  [Binary: 4799]<br>Objective nonzeros   :    4799<br>Linear constraints   :       2  [Less: 2]<br>  Nonzeros           :    9597<br>  RHS nonzeros       :       2<br>Variables            : Min LB: 0.000000         Max UB: 1.000000       <br>Objective nonzeros   : Min   : 2.000000         Max   : 2.500000       <br>Linear constraints   :<br>  Nonzeros           : Min   : 1.000000         Max   : 8.000000       <br>  RHS nonzeros       : Min   : 2399.000         Max   : 3077.500       <br></code><br/><code style="color:#777">CPLEX&gt; Version identifier: 22.1.1.0 | 2022-11-28 | 9160aff4d<br>CPXPARAM_TimeLimit                               300<br>CPXPARAM_MIP_Tolerances_MIPGap                   0.0050000000000000001<br>Found incumbent of value 0.000000 after 0.00 sec. (0.12 ticks)<br>Tried aggregator 1 time.<br>MIP Presolve eliminated 0 rows and 2358 columns.<br>MIP Presolve modified 3 coefficients.<br>Reduced MIP has 2 rows, 2441 columns, and 4881 nonzeros.<br>Reduced MIP has 1093 binaries, 1348 generals, 0 SOSs, and 0 indicators.<br>Presolve time = 0.01 sec. (4.64 ticks)<br></code><br/><code style="color:#777">Tried aggregator 1 time.<br>Detecting symmetries...<br>Reduced MIP has 2 rows, 2441 columns, and 4881 nonzeros.<br>Reduced MIP has 1093 binaries, 1348 generals, 0 SOSs, and 0 indicators.<br>Presolve time = 0.01 sec. (2.45 ticks)<br></code><br/><code style="color:#777">Probing time = 0.00 sec. (0.19 ticks)<br></code><br/><code style="color:#777">MIP emphasis: balance optimality and feasibility.<br>MIP search method: dynamic search.<br>Parallel mode: deterministic, using up to 4 threads.<br></code><br/><code style="color:#777">Root relaxation solution time = 0.01 sec. (1.72 ticks)<br></code><br/><code style="color:#777"><br>        Nodes                                         Cuts/<br>   Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap<br>*     0+    0                            0.0000    10797.7467              --- <br></code><br/><code style="color:#777">*     0     0      integral     0     2848.8583     2848.8583        6    0.00%<br>Elapsed time = 0.04 sec. (11.93 ticks, tree = 0.00 MB, solutions = 2)<br></code><br/><code style="color:#777"><br>Root node processing (before b&amp;c):<br>  Real time             =    0.04 sec. (12.11 ticks)<br>Parallel b&amp;c, 4 threads:<br>  Real time             =    0.00 sec. (0.00 ticks)<br>  Sync time (average)   =    0.00 sec.<br>  Wait time (average)   =    0.00 sec.<br>                          ------------<br>Total (root+branch&amp;cut) =    0.04 sec. (12.11 ticks)<br>Solution pool: 2 solutions saved.<br>MIP - Integer optimal solution:  Objective =  2.8488582500e+03<br>Solution time =    0.04 sec.  Iterations = 6  Nodes = 0<br>Deterministic time = 12.12 ticks  (299.14 ticks/sec)<br></code><br/><code style="color:#777">CPLEX&gt; Incumbent solution written to file '/tmp/tmpm5v51crs.cplex.sol'.<br></code><br/><code style="color:#777">CPLEX&gt; </code><br/><code style="color:#777"></code><br/><code style="color:#777"></code><br/>Solver finished!<br/><code style="color:#777">Solver status=&lt;SolverStatus.ok: 'ok'&gt; and termination_condition=&lt;TerminationCondition.optimal: 'optimal'&gt;<br>good enough</code><br/>Results: {'OUTPUT': 'OUT_LAYER_cceff504_0761_479d_8c76_59259f7b19b8'}<br/><span style="color:#777">Showing layer /tmp/processing_MoaEFQ/0ef1ab26b5dc4c0c9c463ca45fe8752c/OUT_LAYER.gpkg</span><br/><span style="color:#777">Plots:</span><br/><span style="color:#777">file:///tmp/processing_MoaEFQ/0ef1ab26b5dc4c0c9c463ca45fe8752c/observations.png</span><br/><img src="/tmp/processing_MoaEFQ/0ef1ab26b5dc4c0c9c463ca45fe8752c/observations.png"><br/><span style="color:#777">file:///tmp/processing_MoaEFQ/0ef1ab26b5dc4c0c9c463ca45fe8752c/scaled.png</span><br/><img src="/tmp/processing_MoaEFQ/0ef1ab26b5dc4c0c9c463ca45fe8752c/scaled.png"><br/><span style="color:#777">file:///tmp/processing_MoaEFQ/0ef1ab26b5dc4c0c9c463ca45fe8752c/scaled_weighted.png</span><br/><img src="/tmp/processing_MoaEFQ/0ef1ab26b5dc4c0c9c463ca45fe8752c/scaled_weighted.png"><br/><span style="color:#777">file:///tmp/processing_MoaEFQ/0ef1ab26b5dc4c0c9c463ca45fe8752c/solution.png</span><br/><img src="/tmp/processing_MoaEFQ/0ef1ab26b5dc4c0c9c463ca45fe8752c/solution.png"><br/>Results:<br/><code>&nbsp;&nbsp;OUT_LAYER: /tmp/processing_MoaEFQ/0ef1ab26b5dc4c0c9c463ca45fe8752c/OUT_LAYER.gpkg</code><br/>this output is written to: /tmp/algorithm_multirasterknapsack_log_3tlmt28n.html<br/>

tests/mok/instance.py

@@ -0,0 +1,21 @@
+<PAMDataset>
+  <PAMRasterBand band="1">
+    <Histograms>
+      <HistItem>
+        <HistMin>0.9965000000000001</HistMin>
+        <HistMax>8.003500000000001</HistMax>
+        <BucketCount>1000</BucketCount>
+        <IncludeOutOfRange>0</IncludeOutOfRange>
+        <Approximate>0</Approximate>
+        <HistCounts>1770|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|1544|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|1107|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|378</HistCounts>
+      </HistItem>
+    </Histograms>
+    <Metadata>
+      <MDI key="STATISTICS_MAXIMUM">8</MDI>
+      <MDI key="STATISTICS_MEAN">2.565117732861</MDI>
+      <MDI key="STATISTICS_MINIMUM">1</MDI>
+      <MDI key="STATISTICS_STDDEV">1.9513461411422</MDI>
+      <MDI key="STATISTICS_VALID_PERCENT">99.98</MDI>
+    </Metadata>
+  </PAMRasterBand>
+</PAMDataset>