issue 384 netcdf add scale and offset bands metadata

openeogeotrellis/geopysparkdatacube.py

@@ -2160,6 +2160,9 @@ def add_gdalinfo_objects(assets_original):
             nodata = max_level.layer_metadata.no_data_value
             global_metadata = format_options.get("file_metadata",{})
             zlevel = format_options.get("ZLEVEL", 6)
+            for band_name, band_metadata in bands_metadata.items():
+                for tag, value in band_metadata.items():
+                    bands_metadata[band_name][tag] = str(value)
 
             if batch_mode and sample_by_feature:
                 _log.info("Output one netCDF file per feature.")
@@ -2178,6 +2181,7 @@ def add_gdalinfo_objects(assets_original):
                         band_names,
                         dim_names,
                         global_metadata,
+                        bands_metadata,
                         filename_prefix,
                     )
                 else:
@@ -2189,6 +2193,7 @@ def add_gdalinfo_objects(assets_original):
                         band_names,
                         dim_names,
                         global_metadata,
+                        bands_metadata,
                         filename_prefix,
                     )
 
@@ -2203,6 +2208,7 @@ def add_gdalinfo_objects(assets_original):
                         options.setBandNames(band_names)
                         options.setDimensionNames(dim_names)
                         options.setAttributes(global_metadata)
+                        options.setBandsMetadata(bands_metadata)
                         options.setZLevel(zlevel)
                         options.setCropBounds(crop_extent)
                         asset_paths = get_jvm().org.openeo.geotrellis.netcdf.NetCDFRDDWriter.writeRasters(
@@ -2214,14 +2220,14 @@ def add_gdalinfo_objects(assets_original):
                             asset_paths = get_jvm().org.openeo.geotrellis.netcdf.NetCDFRDDWriter.saveSingleNetCDF(max_level.srdd.rdd(),
                                 filename,
                                 band_names,
-                                dim_names,global_metadata,zlevel
+                                dim_names,global_metadata, bands_metadata,zlevel
                             )
                         else:
                             asset_paths = get_jvm().org.openeo.geotrellis.netcdf.NetCDFRDDWriter.saveSingleNetCDFSpatial(
                                 max_level.srdd.rdd(),
                                 filename,
                                 band_names,
-                                dim_names, global_metadata, zlevel
+                                dim_names, global_metadata, bands_metadata, zlevel
                                 )
                     return return_netcdf_assets(asset_paths, bands, nodata)
 

tests/test_api_result.py

@@ -5000,6 +5000,61 @@ def test_custom_geotiff_tags(api110, tmp_path):
     assert second_band.GetScale() == 1.0
     assert second_band.GetOffset() == 4.56
 
+@pytest.mark.parametrize("strict_cropping", [True, False])
+def test_custom_netcdf_tags(api110, tmp_path, strict_cropping):
+    process_graph = {
+            "loadcollection1": {
+                "process_id": "load_collection",
+                "arguments": {
+                    "id": "TestCollection-LonLat16x16",
+                    "temporal_extent": ["2021-01-05", "2021-01-06"],
+                    "spatial_extent": {"west": 0.0, "south": 50.0, "east": 5.0, "north": 55.0},
+                    "bands": ["Flat:1", "Flat:2"],
+                },
+            },
+            "saveresult1": {
+                "process_id": "save_result",
+                "arguments": {
+                    "data": {"from_node": "loadcollection1"},
+                    "format": "NETCDF",
+                    "options": {
+                        "bands_metadata": {
+                            "Flat:1": {
+                                "SCALE": 1.23,
+                                "ARBITRARY": "value",
+                            },
+                            "Flat:2": {
+                                "OFFSET": 4.56,
+                            },
+                            "Flat:3": {
+                                "SCALE": 7.89,
+                                "OFFSET": 10.11
+                            }
+                        },
+                        "strict_cropping": strict_cropping,
+                        "geometries":{"type": "Polygon", "coordinates": [[[0.1, 0.1], [1.8, 0.1], [1.1, 1.8], [0.1, 0.1]]]},
+                    },
+                },
+                "result": True,
+            },
+        }
+    res = api110.result(process_graph).assert_status_code(200)
+    res_path = tmp_path / "res.nc"
+    res_path.write_bytes(res.data)
+    ds = xarray.load_dataset(res_path)
+    assert ds.dims == {"t": 1, "x": 80, "y": 80}
+    assert numpy.datetime_as_string(ds.coords["t"].values, unit="D").tolist() == ["2021-01-05"]
+    assert list(ds.data_vars.keys())[1:] == ["Flat:1", "Flat:2"]
+    assert "long_name" in ds["Flat:1"].attrs
+    assert "units" in ds["Flat:1"].attrs
+    # assert "scale_factor" in ds["Flat:1"].attrs
+    # assert ds["Flat:1"].attrs["scale_factor"] == 1.23
+    assert "long_name" in ds["Flat:2"].attrs
+    assert "units" in ds["Flat:2"].attrs
+    assert "grid_mapping" in ds["Flat:2"].attrs
+    # assert "add_offset" in ds["Flat:2"].attrs
+    # assert ds["Flat:2"].attrs["add_offset"] == 4.56
+
 
 def test_reduce_bands_to_geotiff(api110, tmp_path):
     response = api110.check_result({

tests/test_download.py

@@ -234,10 +234,11 @@ def test_write_assets_samples_tile_grid_batch(self, tmp_path):
     @pytest.mark.parametrize("catalog", [False, True])
     @pytest.mark.parametrize("sample_by_feature", [False, True])
     @pytest.mark.parametrize("format_arg", ["netCDF"])  # "GTIFF" behaves different from "netCDF", so not testing now
+    @pytest.mark.parametrize("bands_metadata", [{},"s","so"])
     def test_write_assets_parameterize_batch(self, tmp_path, imagecollection_with_two_bands_and_three_dates,
                                              imagecollection_with_two_bands_spatial_only,
                                              format_arg, sample_by_feature, catalog, stitch, space_type,
-                                             tile_grid, filename_prefix):
+                                             tile_grid, filename_prefix,bands_metadata):
         d = locals()
         d = {i: d[i] for i in d if i != 'self' and i != "tmp_path" and i != "d"}
         test_name = "-".join(map(str, list(d.values())))  # a bit like how pytest names it
@@ -247,6 +248,10 @@ def test_write_assets_parameterize_batch(self, tmp_path, imagecollection_with_tw
         else:
             imagecollection = imagecollection_with_two_bands_spatial_only
 
+        if bands_metadata == "s":
+            bands_metadata = {"red":{"SCALE":1.23}}
+        elif bands_metadata == "so":
+            bands_metadata = {"red": {"SCALE": 1.23,"OFFSET":4.56}}
         geometries = geojson_to_geometry(self.features)
         assets = imagecollection.write_assets(
             str(tmp_path / "ignored<\0>.extension"),  # null byte to cause error if filename would be written to fs
@@ -259,6 +264,7 @@ def test_write_assets_parameterize_batch(self, tmp_path, imagecollection_with_tw
                 "filename_prefix": filename_prefix,
                 "stitch": stitch,
                 "tile_grid": tile_grid,
+                "bands_metadata": bands_metadata
             }
         )
         with open(self.test_write_assets_parameterize_batch_path + test_name + ".json", 'w') as fp: