Merge pull request #288 from satellogic/optcrop

Add option to disable crop in merge_all

setup.cfg

@@ -5,7 +5,10 @@ exclude = versioneer.py
 [tool:pytest]
 addopts = --cov telluric --cov-report=term-missing --verbose
 
-[mypy-telluric._version]
+[mypy]
+# empty section required in 0.800; https://github.com/python/mypy/issues/9940
+
+[mypy-*._version]
 ignore_errors = True
 
 # See the docstring in versioneer.py for instructions. Note that you must

setup.py

@@ -24,7 +24,7 @@
         'fiona>=1.8.4,<2.*',
         'pyproj<3; python_version>="3.5"',
         'shapely>=1.6.3,<2.*',
-        'rasterio>=1.0.21,<2.*',
+        'rasterio>=1.0.21,< 1.2.0',
         'pillow',
         'mercantile>=0.10.0',
         'boltons',

telluric/_version.py

@@ -6,7 +6,7 @@
 # that just contains the computed version number.
 
 # This file is released into the public domain. Generated by
-# versioneer-0.18 (https://github.com/warner/python-versioneer)
+# versioneer-0.19 (https://github.com/python-versioneer/python-versioneer)
 
 """Git implementation of _version.py."""
 
@@ -57,7 +57,7 @@ class NotThisMethod(Exception):
 
 
 def register_vcs_handler(vcs, method):  # decorator
-    """Decorator to mark a method as the handler for a particular VCS."""
+    """Create decorator to mark a method as the handler of a VCS."""
     def decorate(f):
         """Store f in HANDLERS[vcs][method]."""
         if vcs not in HANDLERS:
@@ -93,9 +93,7 @@ def run_command(commands, args, cwd=None, verbose=False, hide_stderr=False,
         if verbose:
             print("unable to find command, tried %s" % (commands,))
         return None, None
-    stdout = p.communicate()[0].strip()
-    if sys.version_info[0] >= 3:
-        stdout = stdout.decode()
+    stdout = p.communicate()[0].strip().decode()
     if p.returncode != 0:
         if verbose:
             print("unable to run %s (error)" % dispcmd)
@@ -165,6 +163,10 @@ def git_versions_from_keywords(keywords, tag_prefix, verbose):
         raise NotThisMethod("no keywords at all, weird")
     date = keywords.get("date")
     if date is not None:
+        # Use only the last line.  Previous lines may contain GPG signature
+        # information.
+        date = date.splitlines()[-1]
+
         # git-2.2.0 added "%cI", which expands to an ISO-8601 -compliant
         # datestamp. However we prefer "%ci" (which expands to an "ISO-8601
         # -like" string, which we must then edit to make compliant), because
@@ -300,6 +302,9 @@ def git_pieces_from_vcs(tag_prefix, root, verbose, run_command=run_command):
     # commit date: see ISO-8601 comment in git_versions_from_keywords()
     date = run_command(GITS, ["show", "-s", "--format=%ci", "HEAD"],
                        cwd=root)[0].strip()
+    # Use only the last line.  Previous lines may contain GPG signature
+    # information.
+    date = date.splitlines()[-1]
     pieces["date"] = date.strip().replace(" ", "T", 1).replace(" ", "", 1)
 
     return pieces
@@ -338,18 +343,18 @@ def render_pep440(pieces):
 
 
 def render_pep440_pre(pieces):
-    """TAG[.post.devDISTANCE] -- No -dirty.
+    """TAG[.post0.devDISTANCE] -- No -dirty.
 
     Exceptions:
-    1: no tags. 0.post.devDISTANCE
+    1: no tags. 0.post0.devDISTANCE
     """
     if pieces["closest-tag"]:
         rendered = pieces["closest-tag"]
         if pieces["distance"]:
-            rendered += ".post.dev%d" % pieces["distance"]
+            rendered += ".post0.dev%d" % pieces["distance"]
     else:
         # exception #1
-        rendered = "0.post.dev%d" % pieces["distance"]
+        rendered = "0.post0.dev%d" % pieces["distance"]
     return rendered
 
 
@@ -385,7 +390,7 @@ def render_pep440_old(pieces):
 
     The ".dev0" means dirty.
 
-    Eexceptions:
+    Exceptions:
     1: no tags. 0.postDISTANCE[.dev0]
     """
     if pieces["closest-tag"]:

telluric/georaster.py

@@ -128,11 +128,16 @@ def _dest_resolution(first_raster, crs):
 
 def merge_all(rasters, roi=None, dest_resolution=None, merge_strategy=MergeStrategy.UNION,
               shape=None, ul_corner=None, crs=None, pixel_strategy=PixelStrategy.FIRST,
-              resampling=Resampling.nearest):
-    """Merge a list of rasters, cropping by a region of interest.
+              resampling=Resampling.nearest, crop=True):
+    """Merge a list of rasters, cropping (optional) by a region of interest.
        There are cases that the roi is not precise enough for this cases one can use,
        the upper left corner the shape and crs to precisely define the roi.
-       When roi is provided the ul_corner, shape and crs are ignored
+       When roi is provided the ul_corner, shape and crs are ignored.
+
+       NB: Reading rotated rasters with GDAL (and rasterio) gives unpredictable result
+       and in order to overcome this you must use the warping algorithm to apply the rotation (it
+       might be acomplished by using gdalwarp utility). Hence we should have the possibility to
+       disable cropping, otherwise calling merge_all on rotated rasters may cause fails.
     """
 
     first_raster = rasters[0]
@@ -148,8 +153,15 @@ def merge_all(rasters, roi=None, dest_resolution=None, merge_strategy=MergeStrat
         dtype=first_raster.dtype, shape=shape, ul_corner=ul_corner, crs=crs)
 
     # Create a list of single band rasters
+    if not crop:
+        warnings.warn(
+            "The option to disable crop has been added to overcome rare issues that happen "
+            "while working with rotated rasters and it is not yet well tested.",
+            stacklevel=2
+        )
+
     all_band_names, projected_rasters = _prepare_rasters(rasters, merge_strategy, empty,
-                                                         resampling=resampling)
+                                                         resampling=resampling, crop=crop)
     assert len(projected_rasters) == len(rasters)
 
     prepared_rasters = _apply_pixel_strategy(projected_rasters, pixel_strategy)
@@ -221,8 +233,14 @@ def key(rs):
     return rasters_final
 
 
-def _prepare_rasters(rasters, merge_strategy, first, resampling=Resampling.nearest):
-    # type: (List[GeoRaster2], MergeStrategy, GeoRaster2, Resampling) -> Tuple[IndexedSet[str], List[Optional[_Raster]]]
+def _prepare_rasters(
+    rasters,  # type: List[GeoRaster2]
+    merge_strategy,  # type: MergeStrategy
+    first,  # type: GeoRaster2
+    resampling=Resampling.nearest,  # type: Resampling
+    crop=True,  # type: bool
+):
+    # type: (...) -> Tuple[IndexedSet[str], List[Optional[_Raster]]]
     """Prepares the rasters according to the baseline (first) raster and the merge strategy.
 
     The baseline (first) raster is used to crop and reproject the other rasters,
@@ -234,7 +252,7 @@ def _prepare_rasters(rasters, merge_strategy, first, resampling=Resampling.neare
     all_band_names = IndexedSet(first.band_names)
     projected_rasters = []
     for raster in rasters:
-        projected_raster = _prepare_other_raster(first, raster, resampling=resampling)
+        projected_raster = _prepare_other_raster(first, raster, resampling=resampling, crop=crop)
 
         # Modify the bands only if an intersecting raster was returned
         if projected_raster:
@@ -265,22 +283,23 @@ def _explode_raster(raster, band_names=[]):
     return [_Raster(image=raster.bands_data([band_name]), band_names=[band_name]) for band_name in band_names]
 
 
-def _prepare_other_raster(one, other, resampling=Resampling.nearest):
-    # type: (GeoRaster2, GeoRaster2, Resampling) -> Union[_Raster, None]
-    # Crop and reproject the second raster, if necessary
+def _prepare_other_raster(one, other, resampling=Resampling.nearest, crop=True):
+    # type: (GeoRaster2, GeoRaster2, Resampling, bool) -> Union[_Raster, None]
+    # Crop and reproject the second raster, if necessary.
     if not (one.crs == other.crs and one.affine.almost_equals(other.affine) and one.shape == other.shape):
         if one.footprint().intersects(other.footprint()):
-            if one.crs != other.crs:
-                src_bounds = one.footprint().get_bounds(other.crs)
-                src_vector = GeoVector(Polygon.from_bounds(*src_bounds), other.crs)
-                src_width, src_height = (
-                    src_bounds.right - src_bounds.left,
-                    src_bounds.top - src_bounds.bottom)
-                buffer_ratio = int(os.environ.get("TELLURIC_MERGE_CROP_BUFFER", 10))
-                buffer_size = max(src_width, src_height) * (buffer_ratio / 100)
-                other = other.crop(src_vector.buffer(buffer_size))
-            else:
-                other = other.crop(one.footprint(), resolution=one.resolution())
+            if crop:
+                if one.crs != other.crs:
+                    src_bounds = one.footprint().get_bounds(other.crs)
+                    src_vector = GeoVector(Polygon.from_bounds(*src_bounds), other.crs)
+                    src_width, src_height = (
+                        src_bounds.right - src_bounds.left,
+                        src_bounds.top - src_bounds.bottom)
+                    buffer_ratio = int(os.environ.get("TELLURIC_MERGE_CROP_BUFFER", 10))
+                    buffer_size = max(src_width, src_height) * (buffer_ratio / 100)
+                    other = other.crop(src_vector.buffer(buffer_size))
+                else:
+                    other = other.crop(one.footprint(), resolution=one.resolution())
 
             other = other._reproject(new_width=one.width, new_height=one.height,
                                      dest_affine=one.affine, dst_crs=one.crs,
@@ -462,13 +481,12 @@ def __init__(self, image=None, band_names=None, shape=None, nodata=None):
         self._band_names = None
         self._blockshapes = None
         self._shape = copy(shape)
+        self._dtype = None
         if band_names:
             self._set_bandnames(copy(band_names))
         if image is not None:
             self._set_image(image.copy(), nodata)
             self._dtype = np.dtype(image.dtype)
-        else:
-            self._dtype = None
 
     def _build_masked_array(self, image, nodata):
         return np.ma.masked_array(image, image == nodata)
@@ -484,7 +502,7 @@ def _set_image(self, image, nodata=None):
         # convert to masked array:
         if isinstance(image, np.ma.core.MaskedArray):
             masked = image
-        elif isinstance(image, np.core.ndarray):
+        elif isinstance(image, np.ndarray):
             masked = self._build_masked_array(image, nodata)
         else:
             raise GeoRaster2NotImplementedError('only ndarray or masked array supported, got %s' % type(image))
@@ -1101,7 +1119,7 @@ def type_min(dtype):
                 conversion_gain = (omax - omin) / (imax - imin)
 
             # temp conversion, to handle saturation
-            dst_array = conversion_gain * (self.image.astype(np.float) - imin) + omin
+            dst_array = conversion_gain * (self.image.astype(np.float_) - imin) + omin
             dst_array = np.clip(dst_array, omin, omax)
         else:
             dst_array = self.image

telluric/util/raster_utils.py

@@ -28,7 +28,7 @@ def _join_masks_from_masked_array(data):
     """Union of masks."""
     if not isinstance(data.mask, np.ndarray):
         # workaround to handle mask compressed to single value
-        mask = np.empty(data.data.shape, dtype=np.bool)
+        mask = np.empty(data.data.shape, dtype=np.bool_)
         mask.fill(data.mask)
         return mask
     mask = data.mask[0].copy()

tests/test_merge_all.py

@@ -482,15 +482,16 @@ def test_merge_all_non_overlapping_has_correct_metadata():
     assert metadata == expected_metadata
 
 
-def test_merge_all_different_crs():
+@pytest.mark.parametrize("crop", [True, False])
+def test_merge_all_different_crs(crop):
     roi = GeoVector(Polygon.from_bounds(-6321833, -3092272, -6319273, -3089712), WEB_MERCATOR_CRS)
     affine = Affine.translation(-57, -26) * Affine.scale(0.00083, -0.00083)
     expected_resolution = 10
     expected_crs = WEB_MERCATOR_CRS
 
     # from memory
     raster_0 = make_test_raster(1, [1], height=1200, width=1200, affine=affine, crs=WGS84_CRS)
-    result_0 = merge_all([raster_0], roi=roi, dest_resolution=expected_resolution, crs=expected_crs)
+    result_0 = merge_all([raster_0], roi=roi, dest_resolution=expected_resolution, crs=expected_crs, crop=crop)
     assert(result_0.resolution() == expected_resolution)
     assert(result_0.crs == expected_crs)
     assert(result_0.footprint().envelope.almost_equals(roi.envelope, decimal=3))
@@ -499,7 +500,7 @@ def test_merge_all_different_crs():
     path = "/vsimem/raster_for_test.tif"
     result_0.save(path)
     raster_1 = GeoRaster2.open(path)
-    result_1 = merge_all([raster_1], roi=roi, dest_resolution=expected_resolution, crs=expected_crs)
+    result_1 = merge_all([raster_1], roi=roi, dest_resolution=expected_resolution, crs=expected_crs, crop=crop)
 
     assert(result_1.resolution() == expected_resolution)
     assert(result_1.crs == expected_crs)
@@ -508,5 +509,5 @@ def test_merge_all_different_crs():
 
     # preserve the original resolution if dest_resolution is not provided
     raster_2 = make_test_raster(1, [1], height=1200, width=1200, affine=affine, crs=WGS84_CRS)
-    result_2 = merge_all([raster_2], roi=roi, crs=expected_crs)
+    result_2 = merge_all([raster_2], roi=roi, crs=expected_crs, crop=crop)
     assert pytest.approx(result_2.resolution()) == 97.9691