nodata bug fix

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "maturin"
 
 [project]
 name = "NLR-reVRt"
-version = "0.7.0"
+version = "0.7.1"
 description = "National Laboratory of the Rockies' (NLR's) Transmission for reV tool"
 readme = "README.rst"
 authors = [

revrt/costs/cli.py

@@ -226,7 +226,7 @@ def build_routing_layer_file(  # noqa: PLR0917, PLR0913
             ),
             log_runtime("Building routing layers"),
         ):
-            _build_routing_layer_file(
+            out_layer_fp = _build_routing_layer_file(
                 config,
                 lock,
                 validate_masks=validate_masks,
@@ -237,6 +237,8 @@ def build_routing_layer_file(  # noqa: PLR0917, PLR0913
         if client is not None:
             close_dask_client(client)
 
+    return out_layer_fp
+
 
 def _build_routing_layer_file(
     config, lock, validate_masks=False, create_kwargs=None
@@ -262,6 +264,8 @@ def _build_routing_layer_file(
     if config.merge_friction_and_barriers is not None:
         _combine_friction_and_barriers(config, lf_handler, lock=lock)
 
+    return str(lf_handler.fp)
+
 
 def _validated_config(**config_dict):
     """Validate use config inputs"""

revrt/costs/layer_creator.py

@@ -302,8 +302,12 @@ def _process_raster_layer(self, fname, config, tiff_chunks="file"):
             tiff_chunks=tiff_chunks,
             layer_dirs=[self.input_layer_dir, self.output_tiff_dir],
             band_index=0,
+            fillna=None,  # filled below _after_ processing
         )
-        return self._process_raster_data(data, config)
+        data = self._process_raster_data(data, config)
+        if config.na_fill is not None:
+            data = da.where(da.isnan(data), config.na_fill, data)
+        return data
 
     def _process_raster_data(self, data, config):
         """Create the desired layer from the data array"""
@@ -389,6 +393,7 @@ def _process_vector_layer(self, fname, config, tiff_chunks="file"):
             tiff_chunks=tiff_chunks,
             layer_dirs=[self.output_tiff_dir],
             band_index=0,
+            fillna=config.na_fill,
         )
 
     def _rasterize_vector_layer(self, fname, config, out_tiff):
@@ -415,7 +420,7 @@ def _rasterize_vector_layer(self, fname, config, out_tiff):
             all_touched=config.rasterize.all_touched,
             tile_size=config.rasterize.tile_size,
             simply_before_rasterize=config.rasterize.simply_before_rasterize,
-            fill=config.rasterize.fill,
+            fill=config.na_fill,
             dtype=vector_dtype,
             **kwargs,
         )
@@ -427,9 +432,7 @@ def _rasterize_vector_layer(self, fname, config, out_tiff):
 
         logger.debug("Writing rasterized vector '%s' to '%s'", fname, out_tiff)
         save_data_using_layer_file_profile(
-            layer_fp=self._io_handler.fp,
-            data=temp,
-            geotiff=out_tiff,
+            layer_fp=self._io_handler.fp, data=temp, geotiff=out_tiff
         )
 
     def _apply_mask(self, config, data):

revrt/models/cost_layers.py

@@ -140,12 +140,6 @@ class Rasterize(BaseModel, extra="forbid"):
     rasterization. By default, ``False``.
     """
 
-    fill: float | int | None = 0
-    """Value used to fill raster cells not burned by vector
-
-    If None, uses np.nan. By default, ``0``.
-    """
-
 
 class LayerBuildConfig(BaseModel, extra="forbid"):
     """Friction and barrier layers config model
@@ -197,6 +191,9 @@ class LayerBuildConfig(BaseModel, extra="forbid"):
     inclusions are allowed.
     """
 
+    na_fill: float | int | None = 0
+    """Value to fill NA cells with after processing"""
+
 
 def parse_barrier_values(barrier_values):
     """Parse barrier comparison text into an operator and threshold"""

revrt/utilities/base.py

@@ -25,8 +25,7 @@
 
 logger = logging.getLogger(__name__)
 _NUM_GEOTIFF_DIMS = 3  # (band, y, x)
-TRANSFORM_ATOL = 0.01
-"""Tolerance in transform comparison when checking GeoTIFFs"""
+_TRANSFORM_ATOL = 0.01
 
 
 def buffer_routes(
@@ -247,7 +246,13 @@ def file_full_path(file_name, *layer_dirs):
 
 
 def load_data_using_layer_file_profile(
-    layer_fp, geotiff, tiff_chunks="file", layer_dirs=None, band_index=None
+    layer_fp,
+    geotiff,
+    tiff_chunks="file",
+    layer_dirs=None,
+    band_index=None,
+    check_tiff=True,
+    fillna=0,
 ):
     """Load GeoTIFF data, reprojecting to LayeredFile CRS if needed
 
@@ -265,10 +270,22 @@ def load_data_using_layer_file_profile(
         Directories to search for `geotiff` in, if not found in current
         directory. By default, ``None``, which means only the current
         directory is searched.
+    check_tiff : bool, default=True
+        Whether to check that the GeoTIFF profile matches the template
+        defined by the layered file. If ``True``, the GeoTIFF profile is
+        compared to the template, and if they don't match, the GeoTIFF
+        is re-projected to match the template. If ``False``, no checks
+        are performed and the GeoTIFF is loaded as-is.
+        By default, ``True``.
     band_index : int, optional
         Optional index of band to load from the GeoTIFF. If provided,
         only that band will be returned. By default, ``None``, which
         means all bands will be returned.
+    fillna : int or float, optional
+        Value to fill NA cells with after loading the GeoTIFF. If
+        ``None``, NA cells will not be filled (and will typically be
+        represented as ``np.nan`` in the output array).
+        By default, ``0``.
 
     Returns
     -------
@@ -295,20 +312,24 @@ def load_data_using_layer_file_profile(
         "Using the following chunks to open '%s': %r", geotiff, tiff_chunks
     )
 
-    tif = rioxarray.open_rasterio(geotiff, chunks=tiff_chunks)
-    try:
-        check_geotiff(layer_fp, geotiff, transform_atol=TRANSFORM_ATOL)
-    except revrtProfileCheckError:
-        logger.info(
-            "Profile of '%s' does not match template, reprojecting...",
-            geotiff,
-        )
-        geo_box = odc.geo.geobox.GeoBox(
-            shape=(height, width), affine=transform, crs=crs
-        )
-        tif = tif.odc.reproject(
-            how=geo_box, resampling=Resampling.nearest, INIT_DEST=0
-        )
+    tif = rioxarray.open_rasterio(geotiff, chunks=tiff_chunks, masked=True)
+    if check_tiff:
+        try:
+            check_geotiff(layer_fp, geotiff, transform_atol=_TRANSFORM_ATOL)
+        except revrtProfileCheckError:
+            logger.info(
+                "Profile of '%s' does not match template, reprojecting...",
+                geotiff,
+            )
+            geo_box = odc.geo.geobox.GeoBox(
+                shape=(height, width), affine=transform, crs=crs
+            )
+            tif = tif.astype("float32").odc.reproject(
+                how=geo_box, resampling=Resampling.nearest, dst_nodata=np.nan
+            )
+
+    if fillna is not None:
+        tif = tif.fillna(fillna)
 
     if band_index is not None:
         tif = tif.isel(band=band_index)

revrt/utilities/handlers.py

@@ -24,12 +24,11 @@
     revrtValueError,
 )
 from revrt.utilities.base import (
-    check_geotiff,
     delete_data_file,
     expand_dim_if_needed,
     transform_xy,
     save_data_array_to_geotiff,
-    TRANSFORM_ATOL,
+    load_data_using_layer_file_profile,
 )
 from revrt.utilities.monitoring import log_mem, log_runtime
 
@@ -529,21 +528,24 @@ def write_geotiff_to_file(
             geotiff,
         )
         with log_runtime(f"Writing GeoTIFF to {geotiff}"):
-            if check_tiff:
-                logger.debug("\t- Checking %s against %s", geotiff, self.fp)
-                check_geotiff(self.fp, geotiff, transform_atol=TRANSFORM_ATOL)
-
-            with rioxarray.open_rasterio(geotiff, chunks=tiff_chunks) as tif:
-                logger.debug(
-                    "\t- Writing data from %s to %s", geotiff, self.fp
-                )
+            tif = load_data_using_layer_file_profile(
+                self.fp,
+                geotiff,
+                tiff_chunks=tiff_chunks,
+                check_tiff=check_tiff,
+                fillna=None,
+            )
+            logger.debug("\t- Writing data from %s to %s", geotiff, self.fp)
+            try:
                 self.write_layer(
                     tif,
                     layer_name,
                     description=description,
                     overwrite=overwrite,
                     nodata=nodata,
                 )
+            finally:
+                tif.close()
 
     def layer_to_geotiff(
         self, layer, geotiff, ds_chunks="auto", lock=None, **profile_kwargs

tests/python/integration/test_utilties_handlers.py

@@ -260,10 +260,10 @@ def test_geotiff_to_layer_file(tif, tmp_path, test_utility_data_dir):
     lf.write_geotiff_to_file(in_tiff_fp, layer)
 
     with (
-        rioxarray.open_rasterio(in_tiff_fp) as truth_tif,
+        rioxarray.open_rasterio(in_tiff_fp, masked=True) as truth_tif,
         xr.open_dataset(test_fp, consolidated=False, engine="zarr") as ds,
     ):
-        assert np.allclose(ds[layer], truth_tif)
+        assert np.allclose(ds[layer], truth_tif, equal_nan=True)
         assert np.allclose(
             ds[layer].rio.transform(), truth_tif.rio.transform()
         )
@@ -293,12 +293,12 @@ def test_roundtrip(as_list, tmp_path, test_utility_data_dir):
     for layer, truth_fp in layer_names.items():
         test_tiff_fp = tmp_path / f"{layer}.tif"
         with (
-            rioxarray.open_rasterio(truth_fp) as truth_tif,
-            rioxarray.open_rasterio(test_tiff_fp) as test_tif,
+            rioxarray.open_rasterio(truth_fp, masked=True) as truth_tif,
+            rioxarray.open_rasterio(test_tiff_fp, masked=True) as test_tif,
             xr.open_dataset(test_fp, consolidated=False, engine="zarr") as ds,
         ):
-            assert np.allclose(ds[layer], truth_tif)
-            assert np.allclose(test_tif, truth_tif)
+            assert np.allclose(ds[layer], truth_tif, equal_nan=True)
+            assert np.allclose(test_tif, truth_tif, equal_nan=True)
             assert np.allclose(
                 ds[layer].rio.transform(), truth_tif.rio.transform()
             )
@@ -358,14 +358,14 @@ def test_roundtrip_cli(cli_runner, tmp_path, test_utility_data_dir, as_list):
     for layer, truth_fp in layer_names.items():
         test_tiff_fp = tmp_path / f"{layer}.tif"
         with (
-            rioxarray.open_rasterio(truth_fp) as truth_tif,
-            rioxarray.open_rasterio(test_tiff_fp) as test_tif,
+            rioxarray.open_rasterio(truth_fp, masked=True) as truth_tif,
+            rioxarray.open_rasterio(test_tiff_fp, masked=True) as test_tif,
             xr.open_dataset(
                 out_file_fp, consolidated=False, engine="zarr"
             ) as ds,
         ):
-            assert np.allclose(ds[layer], truth_tif)
-            assert np.allclose(test_tif, truth_tif)
+            assert np.allclose(ds[layer], truth_tif, equal_nan=True)
+            assert np.allclose(test_tif, truth_tif, equal_nan=True)
             assert np.allclose(
                 ds[layer].rio.transform(), truth_tif.rio.transform()
             )