[Git][debian-gis-team/rasterio][upstream] New upstream version 1.2.5

Bas Couwenberg (@sebastic) gitlab at salsa.debian.org
Tue Jun 22 05:52:07 BST 2021



Bas Couwenberg pushed to branch upstream at Debian GIS Project / rasterio


Commits:
d781a3fc by Bas Couwenberg at 2021-06-22T06:07:09+02:00
New upstream version 1.2.5
- - - - -


12 changed files:

- CHANGES.txt
- rasterio/__init__.py
- rasterio/_base.pyx
- rasterio/_io.pyx
- rasterio/merge.py
- rasterio/rio/warp.py
- tests/test_complex_dtypes.py
- + tests/test_int8.py
- tests/test_merge.py
- tests/test_plot.py
- tests/test_rio_info.py
- tests/test_rio_warp.py


Changes:

=====================================
CHANGES.txt
=====================================
@@ -1,6 +1,15 @@
 Changes
 =======
 
+1.2.5 (2021-06-21)
+------------------
+
+- Change rio-warp to unrotate imagery by default to match gdalwarp (#2125).
+- Internal to write() cast int8 arrays to uint8 (#2180).
+- Get correct nodata values for complex_int16 data (#2206).
+- Prevent merge failures due to window and slicing mismatches (#2204 and
+  #2202).
+
 1.2.4 (2021-05-31)
 ------------------
 


=====================================
rasterio/__init__.py
=====================================
@@ -40,7 +40,7 @@ import rasterio.enums
 import rasterio.path
 
 __all__ = ['band', 'open', 'pad', 'Env']
-__version__ = "1.2.4"
+__version__ = "1.2.5"
 __gdal_version__ = gdal_version()
 
 # Rasterio attaches NullHandler to the 'rasterio' logger and its


=====================================
rasterio/_base.pyx
=====================================
@@ -493,9 +493,12 @@ cdef class DatasetBase(object):
                     dtype = "float64"
                 elif dtype == "complex64":
                     dtype = "float32"
+                elif dtype == "complex_int16":
+                    dtype = "int16"
 
                 nodataval = GDALGetRasterNoDataValue(band, &success)
                 val = nodataval
+
                 # GDALGetRasterNoDataValue() has two ways of telling you that
                 # there's no nodata value. The success flag might come back
                 # 0 (FALSE). Even if it comes back 1 (TRUE), you still need
@@ -504,9 +507,7 @@ cdef class DatasetBase(object):
                 # there's no nodata value.
                 if dtype not in dtypes.dtype_ranges:
                     pass
-                elif (success == 0 or
-                        val < dtypes.dtype_ranges[dtype][0] or
-                        val > dtypes.dtype_ranges[dtype][1]):
+                elif (success == 0 or val < dtypes.dtype_ranges[dtype][0] or val > dtypes.dtype_ranges[dtype][1]):
                     val = None
                 log.debug(
                     "Nodata success: %d, Nodata value: %f", success, nodataval)


=====================================
rasterio/_io.pyx
=====================================
@@ -1042,6 +1042,7 @@ cdef class DatasetWriterBase(DatasetReaderBase):
         Returns
         -------
         dataset
+
         """
         cdef char **options = NULL
         cdef char *key_c = NULL
@@ -1137,17 +1138,18 @@ cdef class DatasetWriterBase(DatasetReaderBase):
             else:
                 gdal_dtype = dtypes.dtype_rev.get(self._init_dtype)
 
-            if self._init_dtype == 'int8':
+            if _getnpdtype(self._init_dtype) == _getnpdtype('int8'):
                 options = CSLSetNameValue(options, 'PIXELTYPE', 'SIGNEDBYTE')
 
             # Create a GDAL dataset handle.
             try:
-
                 self._hds = exc_wrap_pointer(
-                    GDALCreate(drv, fname, width, height,
-                               count, gdal_dtype, options))
+                    GDALCreate(drv, fname, width, height, count, gdal_dtype, options)
+                )
+
             except CPLE_BaseError as exc:
                 raise RasterioIOError(str(exc))
+
             finally:
                 if options != NULL:
                     CSLDestroy(options)
@@ -1376,6 +1378,9 @@ cdef class DatasetWriterBase(DatasetReaderBase):
 
         dtype = _getnpdtype(dtype)
 
+        if dtype.name == "int8":
+            arr = arr.astype("uint8")
+
         # Require C-continguous arrays (see #108).
         arr = np.require(arr, dtype=dtype, requirements='C')
 
@@ -1755,7 +1760,7 @@ cdef class InMemoryRaster:
                 "in a `with rasterio.Env()` or `with rasterio.open()` "
                 "block.")
 
-        if dtype == 'int8':
+        if _getnpdtype(dtype) == _getnpdtype("int8"):
             options = CSLSetNameValue(options, 'PIXELTYPE', 'SIGNEDBYTE')
 
         datasetname = str(uuid4()).encode('utf-8')
@@ -2019,7 +2024,7 @@ cdef class BufferedDatasetWriterBase(DatasetWriterBase):
             # We've mapped numpy scalar types to GDAL types so see
             # if we can crosswalk those.
             if hasattr(self._init_dtype, 'type'):
-                if self._init_dtype == 'int8':
+                if _getnpdtype(self._init_dtype) == _getnpdtype('int8'):
                     options = CSLSetNameValue(options, 'PIXELTYPE', 'SIGNEDBYTE')
 
                 tp = self._init_dtype.type


=====================================
rasterio/merge.py
=====================================
@@ -70,6 +70,7 @@ def merge(
     output_count=None,
     resampling=Resampling.nearest,
     method="first",
+    target_aligned_pixels=False,
     dst_path=None,
     dst_kwds=None,
 ):
@@ -140,6 +141,10 @@ def merge(
             coff: int
                 column offset in base array
 
+    target_aligned_pixels : bool, optional
+        Whether to adjust output image bounds so that pixel coordinates
+        are integer multiples of pixel size, matching the ``-tap``
+        options of GDAL utilities.  Default: False.
     dst_path : str or Pathlike, optional
         Path of output dataset
     dst_kwds : dict, optional
@@ -217,10 +222,6 @@ def merge(
             ys.extend([bottom, top])
         dst_w, dst_s, dst_e, dst_n = min(xs), min(ys), max(xs), max(ys)
 
-    logger.debug("Output bounds: %r", (dst_w, dst_s, dst_e, dst_n))
-    output_transform = Affine.translation(dst_w, dst_n)
-    logger.debug("Output transform, before scaling: %r", output_transform)
-
     # Resolution/pixel size
     if not res:
         res = first_res
@@ -228,18 +229,19 @@ def merge(
         res = (res, res)
     elif len(res) == 1:
         res = (res[0], res[0])
-    output_transform *= Affine.scale(res[0], -res[1])
-    logger.debug("Output transform, after scaling: %r", output_transform)
+
+    if target_aligned_pixels:
+        dst_w = math.floor(dst_w / res[0]) * res[0]
+        dst_e = math.ceil(dst_e / res[0]) * res[0]
+        dst_s = math.floor(dst_s / res[1]) * res[1]
+        dst_n = math.ceil(dst_n / res[1]) * res[1]
 
     # Compute output array shape. We guarantee it will cover the output
     # bounds completely
-    output_width = int(math.ceil((dst_e - dst_w) / res[0]))
-    output_height = int(math.ceil((dst_n - dst_s) / res[1]))
+    output_width = int(round((dst_e - dst_w) / res[0]))
+    output_height = int(round((dst_n - dst_s) / res[1]))
 
-    # Adjust bounds to fit
-    dst_e, dst_s = output_transform * (output_width, output_height)
-    logger.debug("Output width: %d, height: %d", output_width, output_height)
-    logger.debug("Adjusted bounds: %r", (dst_w, dst_s, dst_e, dst_n))
+    output_transform = Affine.translation(dst_w, dst_n) * Affine.scale(res[0], -res[1])
 
     if dtype is not None:
         dt = dtype
@@ -314,13 +316,17 @@ def merge(
             )
 
             # 4. Read data in source window into temp
-            src_window = src_window.round_shape(pixel_precision=0)
-            dst_window = dst_window.round_shape(pixel_precision=0)
-            trows, tcols = dst_window.height, dst_window.width
-            temp_shape = (src_count, trows, tcols)
-            temp = src.read(
+            src_window_rnd_shp = src_window.round_shape(pixel_precision=0)
+            dst_window_rnd_shp = dst_window.round_shape(pixel_precision=0)
+            dst_window_rnd_off = dst_window_rnd_shp.round_offsets(pixel_precision=0)
+            temp_height, temp_width = (
+                dst_window_rnd_off.height,
+                dst_window_rnd_off.width,
+            )
+            temp_shape = (src_count, temp_height, temp_width)
+            temp_src = src.read(
                 out_shape=temp_shape,
-                window=src_window,
+                window=src_window_rnd_shp,
                 boundless=False,
                 masked=True,
                 indexes=indexes,
@@ -328,9 +334,11 @@ def merge(
             )
 
         # 5. Copy elements of temp into dest
-        dst_window = dst_window.round_offsets(pixel_precision=0)
-        roff, coff = dst_window.row_off, dst_window.col_off
-        region = dest[:, roff:roff + trows, coff:coff + tcols]
+        roff, coff = (
+            max(0, dst_window_rnd_off.row_off),
+            max(0, dst_window_rnd_off.col_off),
+        )
+        region = dest[:, roff : roff + temp_height, coff : coff + temp_width]
 
         if math.isnan(nodataval):
             region_mask = np.isnan(region)
@@ -339,8 +347,9 @@ def merge(
         else:
             region_mask = region == nodataval
 
+        # Ensure common shape, resolving issue #2202.
+        temp = temp_src[:, : region.shape[1], : region.shape[2]]
         temp_mask = np.ma.getmask(temp)
-
         copyto(region, temp, region_mask, temp_mask, index=idx, roff=roff, coff=coff)
 
     if dst_path is None:


=====================================
rasterio/rio/warp.py
=====================================
@@ -1,7 +1,8 @@
 """$ rio warp"""
 
 import logging
-from math import ceil
+from math import ceil, floor
+import sys
 
 import click
 
@@ -148,7 +149,8 @@ def warp(ctx, files, output, driver, like, dst_crs, dimensions, src_bounds,
         setenv(CHECK_WITH_INVERT_PROJ=check_invert_proj)
 
         with rasterio.open(files[0]) as src:
-            l, b, r, t = src.bounds
+            left, bottom, right, top = src.bounds
+
             out_kwargs = src.profile
             out_kwargs.pop("driver", None)
             if driver:
@@ -232,7 +234,7 @@ def warp(ctx, files, output, driver, like, dst_crs, dimensions, src_bounds,
                 # Same projection, different dimensions, calculate resolution.
                 dst_crs = src.crs
                 dst_width, dst_height = dimensions
-                l, b, r, t = src_bounds or (l, b, r, t)
+                l, b, r, t = src_bounds or (left, bottom, right, top)
                 dst_transform = Affine(
                     (r - l) / float(dst_width),
                     0, l, 0,
@@ -249,21 +251,38 @@ def warp(ctx, files, output, driver, like, dst_crs, dimensions, src_bounds,
                 dst_crs = src.crs
                 xmin, ymin, xmax, ymax = (src_bounds or dst_bounds)
                 dst_transform = Affine(res[0], 0, xmin, 0, -res[1], ymax)
-                dst_width = max(int(ceil((xmax - xmin) / res[0])), 1)
-                dst_height = max(int(ceil((ymax - ymin) / res[1])), 1)
+                dst_width = max(int(round((xmax - xmin) / res[0])), 1)
+                dst_height = max(int(round((ymax - ymin) / res[1])), 1)
 
             elif res:
                 # Same projection, different resolution.
                 dst_crs = src.crs
-                dst_transform = Affine(res[0], 0, l, 0, -res[1], t)
-                dst_width = max(int(ceil((r - l) / res[0])), 1)
-                dst_height = max(int(ceil((t - b) / res[1])), 1)
+                dst_transform = Affine(res[0], 0, left, 0, -res[1], top)
+                dst_width = max(int(round((right - left) / res[0])), 1)
+                dst_height = max(int(round((top - bottom) / res[1])), 1)
 
             else:
                 dst_crs = src.crs
-                dst_transform = src.transform
-                dst_width = src.width
-                dst_height = src.height
+                inv_transform = ~src.transform
+                eps = sys.float_info.epsilon
+                c1, r1 = inv_transform * (left + eps, top + eps)
+                c2, r2 = inv_transform * (right + eps, top + eps)
+                c3, r3 = inv_transform * (right + eps, bottom + eps)
+                c4, r4 = inv_transform * (left + eps, bottom + eps)
+                col1 = min(c1, c2, c3, c4)
+                col2 = max(c1, c2, c3, c4)
+                row1 = min(r1, r2, r3, r4)
+                row2 = max(r1, r2, r3, r4)
+                col1 = floor(col1)
+                col2 = ceil(col2)
+                row1 = floor(row1)
+                row2 = ceil(row2)
+                px = (right - left) / (col2 - col1)
+                py = (top - bottom) / (row2 - row1)
+                res = max(px, py)
+                dst_width = max(int(round((right - left) / res)), 1)
+                dst_height = max(int(round((top - bottom) / res)), 1)
+                dst_transform = Affine.translation(left, top) * Affine.scale(res, -res)
 
             if target_aligned_pixels:
                 dst_transform, dst_width, dst_height = aligned_target(dst_transform, dst_width, dst_height, res)


=====================================
tests/test_complex_dtypes.py
=====================================
@@ -81,7 +81,7 @@ def test_complex_int16(tmpdir):
         driver="GTiff",
         height=Z1.shape[0],
         width=Z1.shape[1],
-        nodata=0,
+        nodata=None,
         count=1,
         dtype="complex_int16",
         crs="+proj=latlong",
@@ -94,5 +94,6 @@ def test_complex_int16(tmpdir):
     )
 
     with rasterio.open(tempfile) as dst:
+        assert dst.nodatavals == (None,)
         data = dst.read()
         assert data.dtype == np.complex64


=====================================
tests/test_int8.py
=====================================
@@ -0,0 +1,55 @@
+import affine
+import numpy
+import pytest
+
+import rasterio
+
+
+ at pytest.mark.xfail(reason="Likely upstream bug")
+ at pytest.mark.parametrize("nodata", [-1, -128])
+def test_write_int8_mem(nodata):
+    profile = {
+        "driver": "GTiff",
+        "width": 2,
+        "height": 1,
+        "count": 1,
+        "dtype": "int8",
+        "crs": "EPSG:3857",
+        "transform": affine.Affine(10, 0, 0, 0, -10, 0),
+        "nodata": nodata,
+    }
+
+    values = numpy.array([[nodata, nodata]], dtype="int8")
+
+    with rasterio.open("/vsimem/test.tif", "w", **profile) as src:
+        src.write(values, indexes=1)
+
+    with rasterio.open("/vsimem/test.tif") as src:
+        read = src.read(indexes=1)
+        assert read[0][0] == nodata
+        assert read[0][1] == nodata
+
+
+ at pytest.mark.parametrize("nodata", [None, -1, -128])
+def test_write_int8_fs(tmp_path, nodata):
+    filename = tmp_path.joinpath("test.tif")
+    profile = {
+        "driver": "GTiff",
+        "width": 2,
+        "height": 1,
+        "count": 1,
+        "dtype": "int8",
+        "crs": "EPSG:3857",
+        "transform": affine.Affine(10, 0, 0, 0, -10, 0),
+        "nodata": nodata,
+    }
+
+    values = numpy.array([[127, -128]], dtype="int8")
+
+    with rasterio.open(filename, "w", **profile) as src:
+        src.write(values, indexes=1)
+
+    with rasterio.open(filename) as src:
+        read = src.read(indexes=1)
+        assert read[0][0] == 127
+        assert read[0][1] == -128


=====================================
tests/test_merge.py
=====================================
@@ -1,5 +1,8 @@
 """Tests of rasterio.merge"""
 
+import boto3
+from hypothesis import given, settings
+from hypothesis.strategies import floats
 import numpy
 import pytest
 
@@ -58,7 +61,7 @@ def test_issue2163():
     with rasterio.open("tests/data/float_raster_with_nodata.tif") as src:
         data = src.read()
         result, transform = merge([src])
-        assert numpy.allclose(data, result)
+        assert numpy.allclose(data, result[:, : data.shape[1], : data.shape[2]])
 
 
 def test_unsafe_casting():
@@ -66,3 +69,38 @@ def test_unsafe_casting():
     with rasterio.open("tests/data/float_raster_with_nodata.tif") as src:
         result, transform = merge([src], dtype="uint8", nodata=0.0)
         assert not result.any()  # this is why it's called "unsafe".
+
+
+ at pytest.mark.skipif(
+    not (boto3.Session()._session.get_credentials()),
+    reason="S3 raster access requires credentials",
+)
+ at pytest.mark.network
+ at pytest.mark.slow
+ at settings(deadline=None, max_examples=5)
+ at given(
+    dx=floats(min_value=-0.05, max_value=0.05),
+    dy=floats(min_value=-0.05, max_value=0.05),
+)
+def test_issue2202(dx, dy):
+    import rasterio.merge
+    from shapely import wkt
+    from shapely.affinity import translate
+
+    aoi = wkt.loads(
+        r"POLYGON((11.09 47.94, 11.06 48.01, 11.12 48.11, 11.18 48.11, 11.18 47.94, 11.09 47.94))"
+    )
+    aoi = translate(aoi, dx, dy)
+
+    with rasterio.Env(AWS_NO_SIGN_REQUEST=True,):
+        ds = [
+            rasterio.open(i)
+            for i in [
+                "/vsis3/copernicus-dem-30m/Copernicus_DSM_COG_10_N47_00_E011_00_DEM/Copernicus_DSM_COG_10_N47_00_E011_00_DEM.tif",
+                "/vsis3/copernicus-dem-30m/Copernicus_DSM_COG_10_N48_00_E011_00_DEM/Copernicus_DSM_COG_10_N48_00_E011_00_DEM.tif",
+            ]
+        ]
+        aux_array, aux_transform = rasterio.merge.merge(datasets=ds, bounds=aoi.bounds)
+        from rasterio.plot import show
+
+        show(aux_array)


=====================================
tests/test_plot.py
=====================================
@@ -55,12 +55,15 @@ def test_show_cmyk_interp(tmpdir):
     """A CMYK TIFF has cyan, magenta, yellow, black bands."""
     matplotlib = pytest.importorskip('matplotlib')
     with rasterio.open('tests/data/RGB.byte.tif') as src:
-        meta = src.meta
-    meta['photometric'] = 'cmyk'
-    meta['count'] = 4
-    del meta["nodata"]
+        profile = src.profile
+
+    profile['photometric'] = "CMYK"
+    profile['count'] = 4
+    del profile["nodata"]
+
     tiffname = str(tmpdir.join('foo.tif'))
-    with rasterio.open(tiffname, 'w', **meta) as dst:
+
+    with rasterio.open(tiffname, 'w', **profile) as dst:
         assert dst.colorinterp == (
             ColorInterp.cyan,
             ColorInterp.magenta,


=====================================
tests/test_rio_info.py
=====================================
@@ -1,5 +1,6 @@
 import json
 
+import boto3
 import pytest
 
 import rasterio
@@ -416,7 +417,7 @@ def test_info_no_credentials(tmpdir, monkeypatch, runner):
         ['info', 'tests/data/RGB.byte.tif'])
     assert result.exit_code == 0
 
-
+ at pytest.mark.skipif(not(boto3.Session()._session.get_credentials()), reason="S3 raster access requires credentials")
 @requires_gdal23(reason="Unsigned S3 requests require GDAL ~= 2.3")
 @pytest.mark.network
 def test_info_aws_unsigned(runner):


=====================================
tests/test_rio_warp.py
=====================================
@@ -159,8 +159,8 @@ def test_warp_no_reproject_res(runner, tmpdir):
         with rasterio.open(outputname) as output:
             assert output.crs == src.crs
             assert np.allclose([30, 30], [output.transform.a, -output.transform.e])
-            assert output.width == 327
-            assert output.height == 327
+            assert output.width == 326
+            assert output.height == 326
 
 
 def test_warp_no_reproject_bounds(runner, tmpdir):
@@ -179,7 +179,7 @@ def test_warp_no_reproject_bounds(runner, tmpdir):
             assert np.allclose([src.transform.a, src.transform.e],
                                [output.transform.a, output.transform.e])
             assert output.width == 105
-            assert output.height == 210
+            assert output.height == 209
 
 
 def test_warp_no_reproject_bounds_res(runner, tmpdir):
@@ -196,7 +196,7 @@ def test_warp_no_reproject_bounds_res(runner, tmpdir):
             assert output.crs == src.crs
             assert np.allclose(output.bounds, out_bounds)
             assert np.allclose([30, 30], [output.transform.a, -output.transform.e])
-            assert output.width == 34
+            assert output.width == 33
             assert output.height == 67
 
     # dst-bounds should be an alias to bounds
@@ -609,3 +609,14 @@ def test_warp_resampling_not_yet_supported(
     assert result.exit_code == 2
     assert "Invalid value for" in result.output
     assert "--resampling" in result.output
+
+
+def test_unrotate(runner, tmp_path):
+    """rio-warp unrotates imagery by default, like gdalwarp"""
+    outputname = tmp_path.joinpath("test.tif").as_posix()
+    runner.invoke(main_group, ["warp", "tests/data/rotated.tif", outputname])
+
+    # There is no skew in the output.
+    with rasterio.open(outputname) as src:
+        assert src.transform.b == 0.0
+        assert src.transform.d == 0.0



View it on GitLab: https://salsa.debian.org/debian-gis-team/rasterio/-/commit/d781a3fc856316611c41fa0af7f970c2fe09b5fe

-- 
View it on GitLab: https://salsa.debian.org/debian-gis-team/rasterio/-/commit/d781a3fc856316611c41fa0af7f970c2fe09b5fe
You're receiving this email because of your account on salsa.debian.org.


-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/pkg-grass-devel/attachments/20210622/2a086878/attachment-0001.htm>


More information about the Pkg-grass-devel mailing list