Merge pull request #42 from airware/6.1.0

0.6.1

buzzard/__init__.py

@@ -3,7 +3,7 @@
 buzzard should always be imported the first time from the main thread
 """
 
-__version__ = "0.6.0"
+__version__ = "0.6.1"
 
 # Import osgeo before cv2
 import osgeo as _

buzzard/_a_source_raster.py

@@ -193,7 +193,7 @@ def __init__(self, channels_schema, dtype, fp_stored, **kwargs):
         if self.to_work is not None:
             fp = fp_stored.move(*self.to_work([
                 fp_stored.tl, fp_stored.tr, fp_stored.br
-            ]))
+            ]), round_coordinates=True)
         else:
             fp = fp_stored
 

buzzard/_dataset_back_conversions.py

@@ -136,7 +136,7 @@ def convert_footprint(self, fp, wkt):
         sr = osr.SpatialReference(wkt)
         _, to_virtual = self.get_transforms(sr, fp, 'work')
         if to_virtual:
-            fp = fp.move(*to_virtual([fp.tl, fp.tr, fp.br]))
+            fp = fp.move(*to_virtual([fp.tl, fp.tr, fp.br]), round_coordinates=True)
         return fp
 
     @staticmethod

buzzard/_footprint.py

@@ -24,10 +24,11 @@
 from buzzard._env import env
 from buzzard._footprint_tile import TileMixin
 from buzzard._footprint_intersection import IntersectionMixin
+from buzzard._footprint_move import MoveMixin
 
 LOGGER = logging.getLogger('buzzard')
 
-class Footprint(TileMixin, IntersectionMixin):
+class Footprint(TileMixin, IntersectionMixin, MoveMixin):
     """Immutable object representing the location and size of a spatially localized raster. All
     methods are thread-safe.
 
@@ -412,7 +413,7 @@ def _classify(obj):
             footprints, geoms, resolution, rotation, alignment
         )
 
-    def move(self, tl, tr=None, br=None):
+    def move(self, tl, tr=None, br=None, round_coordinates=False):
         """Create a copy of self moved by an Affine transformation by providing new points.
         `rsize` is always conserved
 
@@ -432,6 +433,19 @@ def move(self, tl, tr=None, br=None):
             New top right coordinates
         br: (nbr, nbr)
             New bottom right coordinates
+        round_coordinates: bool
+            Round the input coordinates with respect to `buzz.env.significant`, so that the output
+            `Footprint` is as much similar as possible as the input `Footprint` regarding those
+            properties:
+            - angle
+            - pxsize
+            - pxsizex / pxsizey
+
+            This option helps a lot if the input coordinates suffered from floating point
+            precision loss since it will cancel the noise in the resulting transformation matrix.
+
+            .. warning::
+                Only work when `tr` and `br` are both provided
 
         Returns
         -------
@@ -453,6 +467,12 @@ def move(self, tl, tr=None, br=None):
             if br is not None:
                 raise ValueError('If br present, bl should be present too') # pragma: no cover
 
+        if round_coordinates:
+            if br is None: # pragma: no cover
+                raise ValueError(
+                    'Can only round input coordinates when all parameters are provided'
+                )
+            tl, tr, br = self._snap_target_coordinates_before_move(tl, tr, br)
 
         if tr is None:
             angle = self.angle
@@ -487,10 +507,16 @@ def move(self, tl, tr=None, br=None):
             affine.Affine.rotation(angle) *
             affine.Affine.scale(*scale)
         )
-        return self.__class__(
-            gt=aff.to_gdal(),
-            rsize=self.rsize,
-        )
+        try:
+            return self.__class__(
+                gt=aff.to_gdal(),
+                rsize=self.rsize,
+            )
+        except ValueError as e:
+            if br is not None and round_coordinates is False and \
+               len(e.args) > 0 and 'north-up' in e.args[0]:
+                raise ValueError('Moving Footprint failed. Try using `round_coordinates=True`.')
+            raise
 
     # Export ************************************************************************************ **
     @property
@@ -1553,23 +1579,7 @@ def _build_neighbors_in_direction(mask, yx_vector):
             if isinstance(mline, shapely.geometry.LineString):
                 mline = sg.MultiLineString([mline])
 
-        # Step 8: Temporary check for badness, until further testing of this method ************* **
-        check = arr.copy()
-        check[:] = 0
-
-        for l in mline:
-            coords = np.asarray(l) - output_offset
-            coords = self.spatial_to_raster(coords)
-            x = coords[:, 0]
-            y = coords[:, 1]
-            check[y, x] = 1
-
-        # Burning size one components since they are not transformed to lines
-        for sly, slx in ndi.find_objects(ndi.label(arr)[0]):
-            if sly.stop - sly.start == 1 and slx.stop - slx.start == 1:
-                check[sly, slx] = 1
-
-        # Step 9: Return ************************************************************************ **
+        # Step 8: Return ************************************************************************ **
         if isinstance(mline, shapely.geometry.LineString):
             return [mline]
         if isinstance(mline, shapely.geometry.MultiLineString):

buzzard/_footprint_move.py

@@ -0,0 +1,108 @@
+""">>> help(TileMixin)"""
+
+import numpy as np
+from buzzard._env import env
+
+class MoveMixin(object):
+    """Private mixin for the Footprint class containing move subroutines"""
+
+    def _snap_target_coordinates_before_move(self, tl1, tr1, br1):
+        rw, rh = self.rsize
+
+        # Extract `source` values *********************************************** **
+        v0 = self.pxlrvec
+        w0 = self.pxtbvec
+        norm_v0, norm_w0 = self.pxsize
+        i0 = v0 / norm_v0
+        j0 = w0 / norm_w0
+
+        # Compute `transformed values` ****************************************** **
+        v1 = (tr1 - tl1) / rw
+        w1 = (br1 - tr1) / rh
+        norm_v1 = np.linalg.norm(v1)
+        norm_w1 = np.linalg.norm(w1)
+        i1 = v1 / norm_v1
+        j1 = w1 / norm_w1
+
+        # Prepare comparison function ******************************************* **
+        largest_coord = np.abs([tl1, tr1, br1]).max()
+        spatial_precision = largest_coord * 10 ** -env.significant
+        def _are_close_enough(p, q):
+            return (np.abs(p - q) < spatial_precision).all()
+
+        # Compute `transformed and snapped values` ****************************** **
+        tl2 = tl1
+
+        # Attempt rounding to preserve `angle` and `pxsize`
+        # e.g. The transformation is just a shift
+        # e.g. The transformation is a flip or a mirror
+        i2 = np.copysign(i0, i1)
+        j2 = np.copysign(j0, j1)
+        norm_v2 = norm_v0
+        norm_w2 = norm_w0
+
+        v2 = i2 * norm_v2
+        w2 = j2 * norm_w2
+        tr2 = tl2 + v2 * rw
+        br2 = tr2 + w2 * rh
+        if _are_close_enough(tr1, tr2) and _are_close_enough(br1, br2):
+            return tl2, tr2, br2
+
+        # Attempt rounding to preserve `angle` and `pxsizex / pxsizey`
+        # e.g. The transformation is a change of unit
+        i2 = np.copysign(i0, i1)
+        j2 = np.copysign(j0, j1)
+        norm_v2 = norm_v1
+        norm_w2 = norm_v1 / norm_v0 * norm_w0
+
+        v2 = i2 * norm_v2
+        w2 = j2 * norm_w2
+        tr2 = tl2 + v2 * rw
+        br2 = tr2 + w2 * rh
+        if _are_close_enough(tr1, tr2) and _are_close_enough(br1, br2):
+            return tl2, tr2, br2
+
+        # Attempt rounding to preserve `angle`
+        # e.g. The transformation is a change of unit different on each axis
+        i2 = np.copysign(i0, i1)
+        j2 = np.copysign(j0, j1)
+        norm_v2 = norm_v1
+        norm_w2 = norm_w1
+
+        v2 = i2 * norm_v2
+        w2 = j2 * norm_w2
+        tr2 = tl2 + v2 * rw
+        br2 = tr2 + w2 * rh
+        if _are_close_enough(tr1, tr2) and _are_close_enough(br1, br2):
+            return tl2, tr2, br2
+
+        # Attempt rounding to preserve `pxsize`
+        # e.g. The transformation is a rotation
+        i2 = i1
+        j2 = j1
+        norm_v2 = norm_v0
+        norm_w2 = norm_w0
+
+        v2 = i2 * norm_v2
+        w2 = j2 * norm_w2
+        tr2 = tl2 + v2 * rw
+        br2 = tr2 + w2 * rh
+        if _are_close_enough(tr1, tr2) and _are_close_enough(br1, br2):
+            return tl2, tr2, br2
+
+        # Attempt rounding to preserve `pxsizex / pxsizey`
+        # e.g. The transformation is a rotation and a change of unit
+        i2 = i1
+        j2 = j1
+        norm_v2 = norm_v1
+        norm_w2 = norm_v1 / norm_v0 * norm_w0
+
+        v2 = i2 * norm_v2
+        w2 = j2 * norm_w2
+        tr2 = tl2 + v2 * rw
+        br2 = tr2 + w2 * rh
+        if _are_close_enough(tr1, tr2) and _are_close_enough(br1, br2):
+            return tl2, tr2, br2
+
+        # Don't perform snapping
+        return tl1, tr1, br1

buzzard/_numpy_raster.py

@@ -73,11 +73,16 @@ def get_data(self, fp, channel_ids, dst_nodata, interpolation):
                 dst_nodata,
                 self.dtype
             )
-        key = list(samplefp.slice_in(self.fp)) + [self._best_indexers_of_channel_ids(channel_ids)]
+        chans_indexer = self._best_indexers_of_channel_ids(channel_ids)
+        key = list(samplefp.slice_in(self.fp)) + [chans_indexer]
         key = tuple(key)
         array = self._arr[key]
         if self._should_tranform:
-            array = array * self.channels_schema['scale'] + self.channels_schema['offset']
+            array = (
+                array *
+                np.asarray(self.channels_schema['scale'])[chans_indexer] +
+                np.asarray(self.channels_schema['offset'])[chans_indexer]
+            )
         array = self.remap(
             samplefp,
             fp,

buzzard/srs/__init__.py

@@ -19,15 +19,14 @@ def wkt_of_any(string):
         return out
     else:
         prj = None
-        with Env(_osgeo_use_exceptions=False):
-            gdal_ds = gdal.OpenEx(string, conv.of_of_str('raster'))
-            if gdal_ds is not None:
-                prj = gdal_ds.GetProjection()
-            gdal_ds = gdal.OpenEx(string, conv.of_of_str('vector'))
-            if gdal_ds is not None:
-                lyr = gdal_ds.GetLayerByIndex(0)
-                if lyr is not None:
-                    prj = lyr.GetSpatialRef()
+        gdal_ds = gdal.OpenEx(string, conv.of_of_str('raster'))
+        if gdal_ds is not None:
+            prj = gdal_ds.GetProjection()
+        gdal_ds = gdal.OpenEx(string, conv.of_of_str('vector'))
+        if gdal_ds is not None:
+            lyr = gdal_ds.GetLayerByIndex(0)
+            if lyr is not None:
+                prj = lyr.GetSpatialRef()
         if prj is not None:
             return prj.ExportToWkt()
     raise ValueError('Could not convert to wkt ({})'.format(str(gdal.GetLastErrorMsg()).strip('\n')))
@@ -41,25 +40,24 @@ def wkt_same(a, b):
     return bool(sra.IsSame(srb))
 
 def _details_of_file(path):
-    with Env(_osgeo_use_exceptions=False):
-        gdal_ds = gdal.OpenEx(path, conv.of_of_str('raster'))
-        if gdal_ds is not None:
-            aff = affine.Affine.from_gdal(*gdal_ds.GetGeoTransform())
-            w, h = gdal_ds.RasterXSize, gdal_ds.RasterYSize
-            cx, cy = aff * [w / 2, h / 2]
-            return gdal_ds.GetProjection(), (cx, cy)
-        gdal_ds = gdal.OpenEx(path, conv.of_of_str('vector'))
-        if gdal_ds is not None:
-            lyr = gdal_ds.GetLayerByIndex(0)
-            if lyr is None:
-                raise ValueError('Could not open file layer')
-            extent = lyr.GetExtent()
-            if extent is None:
-                raise ValueError('Could not compute extent')
-            minx, maxx, miny, maxy = extent
-            cx, cy = (maxx + minx) / 2, (maxy + miny) / 2
-            return lyr.GetSpatialRef().ExportToWkt(), (cx, cy)
-        raise ValueError('Could not open file')
+    gdal_ds = gdal.OpenEx(path, conv.of_of_str('raster'))
+    if gdal_ds is not None:
+        aff = affine.Affine.from_gdal(*gdal_ds.GetGeoTransform())
+        w, h = gdal_ds.RasterXSize, gdal_ds.RasterYSize
+        cx, cy = aff * [w / 2, h / 2]
+        return gdal_ds.GetProjection(), (cx, cy)
+    gdal_ds = gdal.OpenEx(path, conv.of_of_str('vector'))
+    if gdal_ds is not None:
+        lyr = gdal_ds.GetLayerByIndex(0)
+        if lyr is None:
+            raise ValueError('Could not open file layer')
+        extent = lyr.GetExtent()
+        if extent is None:
+            raise ValueError('Could not compute extent')
+        minx, maxx, miny, maxy = extent
+        cx, cy = (maxx + minx) / 2, (maxy + miny) / 2
+        return lyr.GetSpatialRef().ExportToWkt(), (cx, cy)
+    raise ValueError('Could not open file')
 
 def wkt_of_file(path, center=False, unit=None, implicit_unit='m'):
     """Retrieve projection of file as wkt.