Merge ae2338b into 4a292f6

CHANGELOG.rst

@@ -9,6 +9,7 @@ Changelog
 * don't close pool between zoom levels (#227)
 * ``_validate`` module renamed to ``validate`` (#230)
 * fix inverted hillshade & misleading tile reference (#229)
+* fix custom nodata values in overviews (#235)
 
 ----
 0.30

doc/source/apidoc/mapchete.rst

@@ -22,6 +22,7 @@ Submodules
    mapchete.index
    mapchete.log
    mapchete.tile
+   mapchete.validate
 
 Module contents
 ---------------

doc/source/apidoc/mapchete.validate.rst

@@ -0,0 +1,7 @@
+mapchete.validate module
+========================
+
+.. automodule:: mapchete.validate
+    :members:
+    :undoc-members:
+    :show-inheritance:

mapchete/_processing.py

@@ -133,7 +133,7 @@ def _interpolate_from_baselevel(self, baselevel=None):
                     in_affine=parent_tile.affine,
                     out_tile=self.tile,
                     resampling=self.config_baselevels["higher"],
-                    nodataval=self.output_reader.output_params["nodata"]
+                    nodata=self.output_reader.output_params["nodata"]
                 )
             # resample from children tiles
             elif baselevel == "lower":
@@ -150,16 +150,19 @@ def _interpolate_from_baselevel(self, baselevel=None):
                     lower_tiles = [
                         y for y in chain(*[x.get_children() for x in output_tiles])
                     ]
-                mosaic = raster.create_mosaic([
-                    (lower_tile, self.output_reader.read(lower_tile))
-                    for lower_tile in lower_tiles
-                ])
+                mosaic = raster.create_mosaic(
+                    [
+                        (lower_tile, self.output_reader.read(lower_tile))
+                        for lower_tile in lower_tiles
+                    ],
+                    nodata=self.output_reader.output_params["nodata"]
+                )
                 process_data = raster.resample_from_array(
                     in_raster=mosaic.data,
                     in_affine=mosaic.affine,
                     out_tile=self.tile,
                     resampling=self.config_baselevels["lower"],
-                    nodataval=self.output_reader.output_params["nodata"]
+                    nodata=self.output_reader.output_params["nodata"]
                 )
         logger.debug((self.tile.id, "generated from baselevel", str(t)))
         return process_data

mapchete/formats/default/gtiff.py

@@ -286,9 +286,12 @@ def empty(self, process_tile):
         profile = self.profile(process_tile)
         return ma.masked_array(
             data=np.full(
-                (profile["count"], ) + process_tile.shape, profile["nodata"],
-                dtype=profile["dtype"]),
-            mask=True
+                (profile["count"], ) + process_tile.shape,
+                profile["nodata"],
+                dtype=profile["dtype"]
+            ),
+            mask=True,
+            fill_value=profile["nodata"]
         )
 
     def profile(self, tile=None):

mapchete/io/raster.py

@@ -509,7 +509,8 @@ def resample_from_array(
     out_tile=None,
     in_crs=None,
     resampling="nearest",
-    nodataval=0
+    nodataval=None,
+    nodata=0
 ):
     """
     Extract and resample from array to target tile.
@@ -521,32 +522,37 @@ def resample_from_array(
     out_tile : ``BufferedTile``
     resampling : string
         one of rasterio's resampling methods (default: nearest)
-    nodataval : integer or float
+    nodata : integer or float
         raster nodata value (default: 0)
 
     Returns
     -------
     resampled array : array
     """
+    if nodataval is not None:
+        warnings.warn("'nodataval' is deprecated, please use 'nodata'")
+        nodata = nodata or nodataval
     # TODO rename function
     if isinstance(in_raster, ma.MaskedArray):
         pass
-    if isinstance(in_raster, np.ndarray):
-        in_raster = ma.MaskedArray(in_raster, mask=in_raster == nodataval)
+    elif isinstance(in_raster, np.ndarray):
+        in_raster = ma.MaskedArray(in_raster, mask=in_raster == nodata)
     elif isinstance(in_raster, ReferencedRaster):
         in_affine = in_raster.affine
         in_crs = in_raster.crs
         in_raster = in_raster.data
     elif isinstance(in_raster, tuple):
         in_raster = ma.MaskedArray(
             data=np.stack(in_raster),
-            mask=np.stack([
-                band.mask
-                if isinstance(band, ma.masked_array)
-                else np.where(band == nodataval, True, False)
-                for band in in_raster
-            ]),
-            fill_value=nodataval
+            mask=np.stack(
+                [
+                    band.mask
+                    if isinstance(band, ma.masked_array)
+                    else np.where(band == nodata, True, False)
+                    for band in in_raster
+                ]
+            ),
+            fill_value=nodata
         )
     else:
         raise TypeError("wrong input data type: %s" % type(in_raster))
@@ -556,29 +562,32 @@ def resample_from_array(
         pass
     else:
         raise TypeError("input array must have 2 or 3 dimensions")
-    if in_raster.fill_value != nodataval:
-        ma.set_fill_value(in_raster, nodataval)
-    out_shape = (in_raster.shape[0], ) + out_tile.shape
-    dst_data = np.empty(out_shape, in_raster.dtype)
-    in_raster = ma.masked_array(
-        data=in_raster.filled(), mask=in_raster.mask, fill_value=nodataval
+    if in_raster.fill_value != nodata:
+        ma.set_fill_value(in_raster, nodata)
+    dst_data = np.empty(
+        (in_raster.shape[0], ) + out_tile.shape,
+        in_raster.dtype
     )
     reproject(
-        in_raster,
+        in_raster.filled(),
         dst_data,
         src_transform=in_affine,
         src_crs=in_crs or out_tile.crs,
+        src_nodata=nodata,
         dst_transform=out_tile.affine,
         dst_crs=out_tile.crs,
+        dst_nodata=nodata,
         resampling=Resampling[resampling]
     )
-    return ma.MaskedArray(dst_data, mask=dst_data == nodataval)
+    return ma.MaskedArray(dst_data, mask=dst_data == nodata, fill_value=nodata)
 
 
 def create_mosaic(tiles, nodata=0):
     """
-    Create a mosaic from tiles. Tiles must be connected (also possible over Antimeridian),
-    otherwise strange things can happen!
+    Create a mosaic from tiles.
+
+    Tiles must be connected (also possible over Antimeridian), otherwise strange things
+    can happen!
 
     Parameters
     ----------

test/conftest.py

@@ -288,6 +288,13 @@ def baselevels_output_buffer():
     return ExampleConfig(path=path, dict=_dict_from_mapchete(path))
 
 
+@pytest.fixture
+def baselevels_custom_nodata():
+    """Fixture for baselevels_custom_nodata.mapchete."""
+    path = os.path.join(TESTDATA_DIR, "baselevels_custom_nodata.mapchete")
+    return ExampleConfig(path=path, dict=_dict_from_mapchete(path))
+
+
 @pytest.fixture
 def mapchete_input():
     """Fixture for mapchete_input.mapchete."""

test/test_io.py

@@ -342,9 +342,9 @@ def test_write_raster_window_errors():
 def test_extract_from_array():
     """Extract subdata from array."""
     in_tile = BufferedTilePyramid("geodetic", metatiling=4).tile(5, 5, 5)
-    shape = (in_tile.shape[0]//2, in_tile.shape[1])
+    shape = (in_tile.shape[0] // 2, in_tile.shape[1])
     data = ma.masked_array(
-        np.concatenate([np.ones(shape), np.ones(shape)*2])
+        np.concatenate([np.ones(shape), np.ones(shape) * 2])
     )
     # intersecting at top
     out_tile = BufferedTilePyramid("geodetic").tile(5, 20, 20)
@@ -386,6 +386,8 @@ def test_resample_from_array():
     in_data = (np.ones(in_tile.shape[1:]), )
     out_tile = BufferedTilePyramid("geodetic").tile(6, 10, 10)
     out_array = resample_from_array(in_data, in_tile.affine, out_tile)
+    # deprecated
+    resample_from_array(in_data, in_tile.affine, out_tile, nodataval=-9999)
     # errors
     with pytest.raises(TypeError):
         in_data = "invalid_type"

test/test_mapchete.py

@@ -182,6 +182,48 @@ def test_baselevels(mp_tmpdir, baselevels):
         ])
 
 
+def test_baselevels_custom_nodata(mp_tmpdir, baselevels_custom_nodata):
+    """Baselevel interpolation."""
+    fill_value = -32768.0
+    with mapchete.open(baselevels_custom_nodata.path, mode="continue") as mp:
+        # process data before getting baselevels
+        mp.batch_process()
+
+        # get tile from lower zoom level
+        for tile in mp.get_process_tiles(4):
+            lower_tile = mp.get_raw_output(tile)
+            assert not lower_tile.mask.all()
+            # assert fill_value is set and all data are not 0
+            assert lower_tile.fill_value == fill_value
+            assert lower_tile.data.all()
+            # write for next zoom level
+            mp.write(tile, lower_tile)
+            parent_tile = mp.get_raw_output(tile.get_parent())
+            assert not parent_tile.mask.all()
+            # assert fill_value is set and all data are not 0
+            assert parent_tile.fill_value == fill_value
+            assert parent_tile.data.all()
+
+        # get tile from higher zoom level
+        tile = next(mp.get_process_tiles(6))
+        # process and save
+        mp.write(tile, mp.get_raw_output(tile))
+        # read from baselevel
+        assert any([
+            not mp.get_raw_output(upper_tile).mask.all()
+            for upper_tile in tile.get_children()
+        ])
+        # assert fill_value is set and all data are not 0
+        assert all([
+            mp.get_raw_output(upper_tile).fill_value == fill_value
+            for upper_tile in tile.get_children()
+        ])
+        assert all([
+            mp.get_raw_output(upper_tile).data.all()
+            for upper_tile in tile.get_children()
+        ])
+
+
 def test_update_baselevels(mp_tmpdir, baselevels):
     """Baselevel interpolation."""
     conf = dict(baselevels.dict)

test/testdata/baselevels_custom_nodata.mapchete

@@ -0,0 +1,18 @@
+process: ../example_process.py
+zoom_levels:
+    min: 3
+    max: 7
+pyramid:
+    grid: geodetic
+input:
+    file1: dummy1.tif
+output:
+    dtype: float32
+    bands: 3
+    format: GTiff
+    path: tmp/baselevels
+    nodata: -32768.0
+baselevels:
+    min: 5
+    max: 6
+    lower: bilinear