[Git][debian-gis-team/rasterio][upstream] New upstream version 1.0.6
Bas Couwenberg
gitlab at salsa.debian.org
Tue Sep 25 06:21:41 BST 2018
Bas Couwenberg pushed to branch upstream at Debian GIS Project / rasterio
Commits:
0db27e8c by Bas Couwenberg at 2018-09-25T05:02:36Z
New upstream version 1.0.6
- - - - -
9 changed files:
- CHANGES.txt
- docs/topics/plotting.rst
- rasterio/__init__.py
- rasterio/_io.pyx
- rasterio/_warp.pyx
- rasterio/errors.py
- tests/test_boundless_read.py
- tests/test_overviews.py
- tests/test_warp_transform.py
Changes:
=====================================
CHANGES.txt
=====================================
@@ -1,6 +1,20 @@
Changes
=======
+1.0.6 (2018-09-24)
+------------------
+
+Bug fixes:
+
+- If the build_overviews method of a dataset is passed a list of factors that
+ specify more than one 1x1 pixel overview (#1333), Rasterio raises an
+ exception.
+- Calling calculate_default_transform for large extents should no longer result
+ in the out of memory error reported in #1131. The rio-warp command should
+ also now run more quickly and with a smaller memory footprint.
+- We have a more general fix for the problem of filling the empty regions of
+ boundless reads (#1471).
+
1.0.5 (2018-09-19)
------------------
@@ -8,7 +22,7 @@ Bug fixes:
- The fill value for boundless reads was ignored in Rasterio versions 1-1.0.4
but now applies (#1471).
-- An invalid shortcut has been eliminated a Rasterio now prroduces a proper
+- An invalid shortcut has been eliminated and Rasterio now produces a proper
mask in the boundless masked read case (#1449).
- Loss of a row or column in geometry_window() and mask() has been fixed
(#1472).
=====================================
docs/topics/plotting.rst
=====================================
@@ -12,7 +12,6 @@ two dimensional data can be accomplished directly with ``pyplot``.
>>> src = rasterio.open("tests/data/RGB.byte.tif")
>>> pyplot.imshow(src.read(1), cmap='pink')
<matplotlib.image.AxesImage object at 0x...>
- >>> pyplot.show = lambda : None # prevents showing during doctests
>>> pyplot.show()
=====================================
rasterio/__init__.py
=====================================
@@ -43,7 +43,7 @@ import rasterio.path
__all__ = ['band', 'open', 'pad', 'Env']
-__version__ = "1.0.5"
+__version__ = "1.0.6"
__gdal_version__ = gdal_version()
# Rasterio attaches NullHandler to the 'rasterio' logger and its
=====================================
rasterio/_io.pyx
=====================================
@@ -6,6 +6,7 @@ from __future__ import absolute_import
include "directives.pxi"
include "gdal.pxi"
+from collections import Counter
import logging
import sys
import uuid
@@ -24,7 +25,7 @@ from rasterio.enums import ColorInterp, MaskFlags, Resampling
from rasterio.errors import (
CRSError, DriverRegistrationError, RasterioIOError,
NotGeoreferencedWarning, NodataShadowWarning, WindowError,
- UnsupportedOperation
+ UnsupportedOperation, OverviewCreationError
)
from rasterio.sample import sample_gen
from rasterio.transform import Affine
@@ -235,14 +236,7 @@ cdef class DatasetReaderBase(DatasetBase):
log.debug("Output nodata value read from file: %r", ndv)
- # Change given nodatavals to the closest value that
- # can be represented by this band's data type to
- # match GDAL's strategy.
- if fill_value:
- ndv = fill_value
- log.debug("Output nodata value set from fill value")
-
- elif ndv is not None:
+ if ndv is not None:
if np.dtype(dtype).kind in ('i', 'u'):
info = np.iinfo(dtype)
dt_min, dt_max = info.min, info.max
@@ -321,12 +315,6 @@ cdef class DatasetReaderBase(DatasetBase):
else:
out = np.empty(out_shape, dtype=dtype)
- for i, (ndv, arr) in enumerate(zip(
- nodatavals, out if len(out.shape) == 3 else [out])):
-
- if ndv is not None:
- arr.fill(ndv)
-
# Masking
# -------
#
@@ -337,11 +325,10 @@ cdef class DatasetReaderBase(DatasetBase):
# read_masks(), invert them and use them in constructing masked
# arrays.
- if masked:
- enums = self.mask_flag_enums
- all_valid = all([MaskFlags.all_valid in flags for flags in enums])
- log.debug("all_valid: %s", all_valid)
- log.debug("mask_flags: %r", enums)
+ enums = self.mask_flag_enums
+ all_valid = all([MaskFlags.all_valid in flags for flags in enums])
+ log.debug("all_valid: %s", all_valid)
+ log.debug("mask_flags: %r", enums)
# We can jump straight to _read() in some cases. We can ignore
# the boundless flag if there's no given window.
@@ -353,7 +340,7 @@ cdef class DatasetReaderBase(DatasetBase):
out = self._read(indexes, out, window, dtype,
resampling=resampling)
- if masked:
+ if masked or fill_value is not None:
if all_valid:
mask = np.ma.nomask
else:
@@ -365,17 +352,23 @@ cdef class DatasetReaderBase(DatasetBase):
kwds = {'mask': mask}
# Set a fill value only if the read bands share a
# single nodata value.
- if len(set(nodatavals)) == 1:
+ if fill_value is not None:
+ kwds['fill_value'] = fill_value
+ elif len(set(nodatavals)) == 1:
if nodatavals[0] is not None:
kwds['fill_value'] = nodatavals[0]
+
out = np.ma.array(out, **kwds)
+ if not masked:
+ out = out.filled(fill_value)
+
# If this is a boundless read we will create an in-memory VRT
# in order to use GDAL's windowing and compositing logic.
else:
vrt_doc = _boundless_vrt_doc(
- self, nodata=ndv, hidenodata=bool(fill_value), width=max(self.width, window.width) + 1,
+ self, nodata=ndv, width=max(self.width, window.width) + 1,
height=max(self.height, window.height) + 1,
transform=self.window_transform(window)).decode('ascii')
@@ -389,7 +382,7 @@ cdef class DatasetReaderBase(DatasetBase):
indexes, out, Window(0, 0, window.width, window.height),
None, resampling=resampling)
- if masked:
+ if masked or fill_value is not None:
mask = np.zeros(out.shape, 'uint8')
mask = ~vrt._read(
indexes, mask, Window(0, 0, window.width, window.height), None, masks=True).astype('bool')
@@ -397,11 +390,17 @@ cdef class DatasetReaderBase(DatasetBase):
kwds = {'mask': mask}
# Set a fill value only if the read bands share a
# single nodata value.
- if len(set(nodatavals)) == 1:
+ if fill_value is not None:
+ kwds['fill_value'] = fill_value
+ elif len(set(nodatavals)) == 1:
if nodatavals[0] is not None:
kwds['fill_value'] = nodatavals[0]
+
out = np.ma.array(out, **kwds)
+ if not masked:
+ out = out.filled(fill_value)
+
if return2d:
out.shape = out.shape[1:]
@@ -1551,6 +1550,11 @@ cdef class DatasetWriterBase(DatasetReaderBase):
['Resampling.{0}'.format(Resampling(k).name) for k in
resampling_map.keys()])))
+ # Check factors
+ ovr_shapes = Counter([(int((self.height + f - 1) / f), int((self.width + f - 1) / f)) for f in factors])
+ if ovr_shapes[(1, 1)] > 1:
+ raise OverviewCreationError("Too many overviews levels of 1x1 dimension were requested")
+
# Allocate arrays.
if factors:
factors_c = <int *>CPLMalloc(len(factors)*sizeof(int))
=====================================
rasterio/_warp.pyx
=====================================
@@ -6,10 +6,12 @@ include "gdal.pxi"
import logging
import uuid
import warnings
+import xml.etree.ElementTree as ET
from affine import identity
import numpy as np
+import rasterio
from rasterio._base import gdal_version
from rasterio._err import (
CPLE_IllegalArgError, CPLE_NotSupportedError,
@@ -20,6 +22,7 @@ from rasterio.enums import Resampling, MaskFlags, ColorInterp
from rasterio.env import GDALVersion
from rasterio.crs import CRS
from rasterio.errors import (
+ GDALOptionNotImplementedError,
DriverRegistrationError, CRSError, RasterioIOError,
RasterioDeprecationWarning, WarpOptionsError)
from rasterio.transform import Affine, from_bounds, guard_transform, tastes_like_gdal
@@ -31,7 +34,7 @@ from rasterio._err cimport exc_wrap_pointer, exc_wrap_int
from rasterio._io cimport (
DatasetReaderBase, InMemoryRaster, in_dtype_range, io_auto)
from rasterio._features cimport GeomBuilder, OGRGeomBuilder
-from rasterio._shim cimport delete_nodata_value
+from rasterio._shim cimport delete_nodata_value, open_dataset
log = logging.getLogger(__name__)
@@ -520,7 +523,7 @@ def _calculate_default_transform(src_crs, dst_crs, width, height,
cdef double geotransform[6]
cdef OGRSpatialReferenceH osr = NULL
cdef char *wkt = NULL
- cdef InMemoryRaster temp = None
+ cdef GDALDatasetH hds = NULL
extent[:] = [0.0, 0.0, 0.0, 0.0]
geotransform[:] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
@@ -540,36 +543,46 @@ def _calculate_default_transform(src_crs, dst_crs, width, height,
OSRExportToWkt(osr, &wkt)
_safe_osr_release(osr)
- with InMemoryRaster(width=width, height=height, transform=transform,
- gcps=gcps, crs=src_crs) as temp:
+ if isinstance(src_crs, str):
+ src_crs = CRS.from_string(src_crs)
+ elif isinstance(src_crs, dict):
+ src_crs = CRS(**src_crs)
+
+ vrt_doc = _suggested_proxy_vrt_doc(width, height, transform=transform, crs=src_crs, gcps=gcps).decode('ascii')
+
+ try:
try:
- hTransformArg = exc_wrap_pointer(
- GDALCreateGenImgProjTransformer(
- temp._hds, NULL, NULL, wkt, 1, 1000.0,0))
- exc_wrap_int(
- GDALSuggestedWarpOutput2(
- temp._hds, GDALGenImgProjTransform, hTransformArg,
- geotransform, &npixels, &nlines, extent, 0))
-
- log.debug("Created transformer and warp output.")
-
- except CPLE_NotSupportedError as err:
- raise CRSError(err.errmsg)
-
- except CPLE_AppDefinedError as err:
- if "Reprojection failed" in str(err):
- # This "exception" should be treated as a debug msg, not error
- # "Reprojection failed, err = -14, further errors will be
- # suppressed on the transform object."
- log.info("Encountered points outside of valid dst crs region")
- pass
- else:
- raise err
- finally:
- if wkt != NULL:
- CPLFree(wkt)
- if hTransformArg != NULL:
- GDALDestroyGenImgProjTransformer(hTransformArg)
+ hds = open_dataset(vrt_doc, 0x00 | 0x02 | 0x04, ['VRT'], {}, None)
+ except GDALOptionNotImplementedError:
+ hds = open_dataset(vrt_doc, 0x00 | 0x02 | 0x04, None, None, None)
+
+ hTransformArg = exc_wrap_pointer(
+ GDALCreateGenImgProjTransformer(
+ hds, NULL, NULL, wkt, 1, 1000.0,0))
+ exc_wrap_int(
+ GDALSuggestedWarpOutput2(
+ hds, GDALGenImgProjTransform, hTransformArg,
+ geotransform, &npixels, &nlines, extent, 0))
+
+ log.debug("Created transformer and warp output.")
+
+ except CPLE_NotSupportedError as err:
+ raise CRSError(err.errmsg)
+
+ except CPLE_AppDefinedError as err:
+ if "Reprojection failed" in str(err):
+ # This "exception" should be treated as a debug msg, not error
+ # "Reprojection failed, err = -14, further errors will be
+ # suppressed on the transform object."
+ log.info("Encountered points outside of valid dst crs region")
+ pass
+ else:
+ raise err
+ finally:
+ if wkt != NULL:
+ CPLFree(wkt)
+ if hTransformArg != NULL:
+ GDALDestroyGenImgProjTransformer(hTransformArg)
# Convert those modified arguments to Python values.
dst_affine = Affine.from_gdal(*[geotransform[i] for i in range(6)])
@@ -922,3 +935,32 @@ cdef class WarpedVRTReaderBase(DatasetReaderBase):
raise ValueError("WarpedVRT does not permit boundless reads")
else:
return super(WarpedVRTReaderBase, self).read_masks(indexes=indexes, out=out, window=window, out_shape=out_shape, resampling=resampling)
+
+
+def _suggested_proxy_vrt_doc(width, height, transform=None, crs=None, gcps=None):
+ """Make a VRT XML document to serve _calculate_default_transform."""
+ vrtdataset = ET.Element('VRTDataset')
+ vrtdataset.attrib['rasterYSize'] = str(height)
+ vrtdataset.attrib['rasterXSize'] = str(width)
+ vrtrasterband = ET.SubElement(vrtdataset, 'VRTRasterBand')
+
+ srs = ET.SubElement(vrtdataset, 'SRS')
+ srs.text = crs.wkt if crs else ""
+
+ if gcps:
+ gcplist = ET.SubElement(vrtdataset, 'GCPList')
+ gcplist.attrib['Projection'] = crs.wkt if crs else ""
+ for point in gcps:
+ gcp = ET.SubElement(gcplist, 'GCP')
+ gcp.attrib['Id'] = str(point.id)
+ gcp.attrib['Info'] = str(point.info)
+ gcp.attrib['Pixel'] = str(point.col)
+ gcp.attrib['Line'] = str(point.row)
+ gcp.attrib['X'] = str(point.x)
+ gcp.attrib['Y'] = str(point.y)
+ gcp.attrib['Z'] = str(point.z)
+ else:
+ geotransform = ET.SubElement(vrtdataset, 'GeoTransform')
+ geotransform.text = ','.join([str(v) for v in transform.to_gdal()])
+
+ return ET.tostring(vrtdataset)
=====================================
rasterio/errors.py
=====================================
@@ -98,3 +98,7 @@ class WarpOptionsError(RasterioError):
class UnsupportedOperation(RasterioError):
"""Raised when reading from a file opened in 'w' mode"""
+
+
+class OverviewCreationError(RasterioError):
+ """Raised when creation of an overview fails"""
=====================================
tests/test_boundless_read.py
=====================================
@@ -80,3 +80,18 @@ def test_boundless_fill_value():
assert (filled[:, -1] == 5).all()
assert (filled[0, :] == 5).all()
assert (filled[-1, :] == 5).all()
+
+
+def test_boundless_fill_value_overview_masks():
+ """Confirm a more general resolution to issue #1471"""
+ with rasterio.open("tests/data/cogeo.tif") as src:
+ data = src.read(1, boundless=True, window=Window(-300, -335, 1000, 1000), fill_value=5, out_shape=(512, 512))
+ assert (data[:, 0] == 5).all()
+
+
+def test_boundless_masked_fill_value_overview_masks():
+ """Confirm a more general resolution to issue #1471"""
+ with rasterio.open("tests/data/cogeo.tif") as src:
+ data = src.read(1, masked=True, boundless=True, window=Window(-300, -335, 1000, 1000), fill_value=5, out_shape=(512, 512))
+ assert data.fill_value == 5
+ assert data.mask[:, 0].all()
=====================================
tests/test_overviews.py
=====================================
@@ -9,6 +9,8 @@ import pytest
import rasterio
from rasterio.enums import Resampling
from rasterio.env import GDALVersion
+from rasterio.errors import OverviewCreationError
+
gdal_version = GDALVersion()
@@ -79,3 +81,12 @@ def test_test_unsupported_algo(data):
with rasterio.open(inputfile, 'r+') as src:
overview_factors = [2, 4]
src.build_overviews(overview_factors, resampling=Resampling.q1)
+
+
+def test_issue1333(data):
+ """Fail if asked to create more than one 1x1 overview"""
+ inputfile = str(data.join('RGB.byte.tif'))
+ with pytest.raises(OverviewCreationError):
+ with rasterio.open(inputfile, 'r+') as src:
+ overview_factors = [1024, 2048]
+ src.build_overviews(overview_factors, resampling=Resampling.average)
=====================================
tests/test_warp_transform.py
=====================================
@@ -4,6 +4,7 @@ import pytest
from rasterio._warp import _calculate_default_transform
from rasterio.control import GroundControlPoint
+from rasterio.crs import CRS
from rasterio.errors import CRSError
from rasterio.transform import from_bounds
from rasterio.warp import calculate_default_transform, transform_bounds
@@ -60,7 +61,7 @@ def test_identity():
assert res_width == width
assert res_height == height
for res, exp in zip(res_transform, transform):
- assert round(res, 7) == round(exp, 7)
+ assert round(res, 3) == round(exp, 3)
def test_identity_gcps():
@@ -88,7 +89,7 @@ def test_identity_gcps():
assert res_width == width
assert res_height == height
for res, exp in zip(res_transform, transform):
- assert round(res, 7) == round(exp, 7)
+ assert round(res, 3) == round(exp, 3)
def test_transform_bounds():
@@ -173,3 +174,9 @@ def test_transform_bounds_identity():
"""Confirm fix of #1411"""
bounds = (12978395.906596646, 146759.09430753812, 12983287.876406897, 151651.06411778927)
assert transform_bounds("+init=epsg:3857", "+init=epsg:3857", *bounds) == bounds
+
+
+def test_issue1131():
+ """Confirm that we don't run out of memory"""
+ transform, w, h = calculate_default_transform(CRS.from_epsg(4326), CRS.from_epsg(3857), 455880, 454450, 13.0460235139, 42.6925552354, 13.2511695428, 42.8970561511)
+ assert (w, h) == (381595, 518398)
View it on GitLab: https://salsa.debian.org/debian-gis-team/rasterio/commit/0db27e8c3357e7b64133650c78091f873f240bbb
--
View it on GitLab: https://salsa.debian.org/debian-gis-team/rasterio/commit/0db27e8c3357e7b64133650c78091f873f240bbb
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/20180925/962b1ccc/attachment-0001.html>
More information about the Pkg-grass-devel
mailing list