Additional export to GeoTIFF

documentation/conf.py.in

@@ -33,7 +33,8 @@ release = '${campo_VERSION}'
 extensions = [
     'sphinx.ext.autodoc',
     'sphinx.ext.ifconfig',
-    'sphinx.ext.autosummary'
+    'sphinx.ext.autosummary',
+    'sphinx_autodoc_typehints'
 ]
 
 autosummary_generate=True

documentation/operations/dataframe.rst

@@ -168,3 +168,4 @@ Reference
    to_csv
    to_gpkg
    to_tiff
+   to_geotiff

environment/configuration/conda_environment.yaml

@@ -3,10 +3,11 @@ channels:
   - conda-forge
   - nodefaults
 dependencies:
+  - cmake
   - python=3.11
   - numpy
   - sphinx=5.3
-  - cmake
+  - sphinx-autodoc-typehints
   - pandas
   - lue>0.3.6
   - hpx=1.9.0

source/campo/op_experimental/export.py

@@ -262,7 +262,7 @@ def to_gpkg(dataframe, filename, crs='', timestep=None):
 
 
 
-def to_tiff(dataframe, crs='', directory='', timestep=None):
+def to_tiff(dataframe, crs="", directory="", timestep=None):
   """ Exports field agent property to a set of GeoTIFF outputs
 
   :param dataframe: Input dataframe from LUE dataset
@@ -555,3 +555,61 @@ def create_field_pdf(frame, filename):
 
   subprocess.check_call(cmd, shell=True, stdout=subprocess.DEVNULL)
   shutil.rmtree(tmpdir)
+
+
+
+
+
+
+def _gdal_datatype(data_type):
+    """ return corresponding GDAL datatype """
+
+    if data_type == 'bool':
+        return gdal.GDT_Byte
+    elif data_type == 'int32':
+        return gdal.GDT_Int32
+    elif data_type == 'int64':
+        return gdal.GDT_Int64
+    elif data_type == 'float32':
+        return gdal.GDT_Float32
+    elif data_type == 'float64':
+        return gdal.GDT_Float64
+    else:
+        raise ValueError(f"Data type '{data_type}' non-suported")
+
+
+def to_geotiff(dataframe, path: str, crs: str) -> None:
+    """ Exports field agent property to a GeoTIFF
+
+    :param dataframe: Input dataframe from LUE dataset for particular timestep
+    :param path: Output path
+    :param crs: Coordinate Reference System, e.g. "EPSG:4326"
+    """
+
+    if crs != "":
+        aut, code = crs.split(":")
+        if aut != "EPSG":
+            msg = _color_message('Provide CRS like "EPSG:4326"')
+            raise TypeError(msg)
+
+    rows = dataframe.values.shape[0]
+    cols = dataframe.values.shape[1]
+    cellsize_x = math.fabs(dataframe.xcoord[1].values - dataframe.xcoord[0].values)
+    cellsize_y = math.fabs(dataframe.ycoord[1].values - dataframe.ycoord[0].values)
+
+    data = dataframe.data
+
+    xmin = dataframe.xcoord[0].values.item()
+    # last ycoordinate is the lower of the topmost row, we need the upper ycoordinate for the extent
+    ymax = dataframe.ycoord[-1].values.item() + cellsize_y
+    geotransform = (xmin, cellsize_x, 0, ymax, 0, -cellsize_y)
+
+    out_ds = gdal.GetDriverByName('GTiff').Create(str(path), cols, rows, 1, _gdal_datatype(data.dtype))
+    out_ds.SetGeoTransform(geotransform)
+
+    srs = osr.SpatialReference()
+    srs.ImportFromEPSG(int(code))
+    out_ds.SetProjection(srs.ExportToWkt())
+
+    out_ds.GetRasterBand(1).WriteArray(data)
+    out_ds = None

source/test/CMakeLists.txt

@@ -4,6 +4,7 @@ set(SOURCES
   extract_dyn_diff.py
   extract_dyn_same.py
   unit_tests.py
+  test_dynamic_model.py
   test_phenomenon.py
   test_propertyset.py
   test_property.py

source/test/test_dataframe.py

@@ -1,3 +1,4 @@
+import pathlib
 import unittest
 
 import numpy as np
@@ -46,3 +47,20 @@ def test_2(self):
             campo.mobile_points_to_gpkg(coords, tmp_df, f"tmp_{timestep}.gpkg", 'EPSG:28992')
 
 
+    def test_3(self):
+        """ Writing GeoTiff stationary dynamic field agent """
+
+        dataset = ldm.open_dataset("TestDynamicModel_test_1.lue", "r")
+        df = campo.dataframe.select(dataset.phen, property_names=['fdata'])
+
+        agent_id = 0
+        crs = "EPSG:4326"
+        directory = "TestDataframe"
+
+        if not pathlib.Path(directory).exists():
+            pathlib.Path(directory).mkdir()
+
+        for timestep in range(1, 6):
+            raster = df["phen"]["field"]["fdata"][agent_id][timestep - 1]
+            filename = pathlib.Path(directory, f"fdata_{agent_id}_{timestep}.tiff")
+            campo.to_geotiff(raster, filename, crs)

source/test/test_dynamic_model.py

@@ -0,0 +1,69 @@
+import datetime
+import unittest
+
+import numpy as np
+
+import campo
+
+import lue.data_model as ldm
+
+
+class TestDynamicModel(unittest.TestCase):
+
+    @classmethod
+    def tearDownClass(self):
+        pass
+
+    @classmethod
+    def setUpClass(self):
+
+        # Some dummy spatial domain
+        with open('locations.csv', 'w') as content:
+            content.write('1,2\n3,4\n5,6\n7,8')
+
+        with open('extent.csv', 'w') as content:
+            content.write("0,0,30,20,2,3\n")
+            content.write("0,0,30,20,2,3\n")
+            content.write("0,0,30,20,2,3\n")
+            content.write("0,0,30,20,2,3\n")
+
+        self.ds = campo.Campo(seed=13)
+
+        self.phen = self.ds.add_phenomenon('phen')
+        self.phen.add_property_set('point', 'locations.csv')
+        self.phen.add_property_set('field', 'extent.csv')
+
+
+
+    def test_1(self):
+        """ Dynamic model, stationary points and fields """
+
+        # initial section
+        self.phen.point.pdata = 500
+        self.phen.field.fdata = 300
+
+        # dummy time
+        date = datetime.date(2000, 1, 1)
+        time = datetime.time(0, 0)
+        start = datetime.datetime.combine(date, time)
+        unit = campo.TimeUnit.month
+        stepsize = 1
+        timesteps = 5
+
+        # create the output lue data set
+        self.ds.create_dataset("TestDynamicModel_test_1.lue")
+        self.ds.set_time(start, unit, stepsize, timesteps)
+
+
+        self.phen.point.pdata.is_dynamic = True
+        self.phen.field.fdata.is_dynamic = True
+
+        # write after modifying coordinates
+        self.ds.write()
+
+        # dynamic section
+        for timestep in range(1, timesteps + 1):
+            self.phen.point.pdata += 200
+            self.phen.field.fdata += 100
+
+            self.ds.write(timestep)

source/test/unit_tests.py

@@ -10,6 +10,7 @@
 import test_property
 import test_mobile_agents
 import test_dataframe
+import test_dynamic_model
 
 
 if __name__ == "__main__":
@@ -26,6 +27,7 @@
     # suite.addTest(unittest.TestLoader().loadTestsFromModule(extract_dyn_diff))
 
     suite.addTest(unittest.TestLoader().loadTestsFromModule(test_mobile_agents))
+    suite.addTest(unittest.TestLoader().loadTestsFromModule(test_dynamic_model))
     suite.addTest(unittest.TestLoader().loadTestsFromModule(test_dataframe))
 
     result = unittest.TextTestRunner(verbosity=3).run(suite)