Improve rpcs support

telluric/georaster.py

@@ -1326,7 +1326,7 @@ def copy_with(self, mutable=None, **kwargs):
         """Get a copy of this GeoRaster with some attributes changed. NOTE: image is shallow-copied!"""
         if mutable is None:
             mutable = isinstance(self, MutableGeoRaster)
-        init_args = {'affine': self.affine, 'crs': self.crs, 'band_names': self.band_names,
+        init_args = {'affine': self.affine, 'crs': self.crs, 'rpcs': self.rpcs, 'band_names': self.band_names,
                      'nodata': self.nodata_value}
         init_args.update(kwargs)
 
@@ -1408,7 +1408,14 @@ def _resize(self, ratio_x, ratio_y, resampling):
         """Return raster resized by ratio."""
         new_width = int(np.ceil(self.width * ratio_x))
         new_height = int(np.ceil(self.height * ratio_y))
-        dest_affine = self.affine * Affine.scale(1 / ratio_x, 1 / ratio_y)
+
+        if self.crs is None and self.rpcs is not None:
+            # resize raster with rpcs
+            dest_rpcs = self._resize_rpcs(ratio_x, ratio_y)
+            dest_affine = self.affine
+        else:
+            dest_affine = self.affine * Affine.scale(1 / ratio_x, 1 / ratio_y)
+            dest_rpcs = self.rpcs
 
         window = rasterio.windows.Window(0, 0, self.width, self.height)
         resized_raster = self.get_window(
@@ -1417,10 +1424,20 @@ def _resize(self, ratio_x, ratio_y, resampling):
             ysize=new_height,
             resampling=resampling,
             affine=dest_affine,
+            rpcs=dest_rpcs
         )
 
         return resized_raster
 
+    def _resize_rpcs(self, ratio_x, ratio_y):
+        """resize raster by using its rpcs"""
+        dest_rpcs = copy(self.rpcs)
+        dest_rpcs.line_off = dest_rpcs.line_off * ratio_y
+        dest_rpcs.samp_off = dest_rpcs.samp_off * ratio_x
+        dest_rpcs.line_scale = dest_rpcs.line_scale * ratio_y
+        dest_rpcs.samp_scale = dest_rpcs.samp_scale * ratio_x
+        return dest_rpcs
+
     def to_pillow_image(self, return_mask=False):
         """Return Pillow. Image, and optionally also mask."""
         img = np.rollaxis(np.rollaxis(self.image.data, 2), 2)
@@ -1707,20 +1724,35 @@ def image_corner(self, corner):
         y = 0 if corner[0] == 'u' else self.height
         return Point(x, y)
 
-    def corner(self, corner):
-        """Return footprint origin in world coordinates, as GeoVector."""
+    def corner(self, corner, **kwargs):
+        """Return footprint origin in world coordinates, as GeoVector.
+           param  kwargs: additional arguments passed to corners function: dem path and dem.crs for rpcs rasters."""
+        if self.crs is None and self.rpcs is not None:
+            new_raster = self.reproject(dst_crs=WGS84_CRS, rpcs=self.rpcs, **kwargs)
+            return new_raster.to_world(new_raster.image_corner(corner))
         return self.to_world(self.image_corner(corner))
 
-    def corners(self):
-        """Return footprint corners, as {corner_type -> GeoVector}."""
-        return {corner: self.corner(corner) for corner in self.corner_types()}
-
-    def origin(self):
-        """Return footprint origin in world coordinates, as GeoVector."""
-        return self.corner('ul')
-
-    def center(self):
-        """Return footprint center in world coordinates, as GeoVector."""
+    def corners(self, **kwargs):
+        """Return footprint corners, as {corner_type -> GeoVector}.
+           :param  kwargs: additional arguments passed to corners function: dem path and dem.crs for rpcs rasters."""
+        if self.crs is None and self.rpcs is not None:
+            new_raster = self.reproject(dst_crs=WGS84_CRS, rpcs=self.rpcs, **kwargs)
+            return {corner: new_raster.corner(corner) for corner in new_raster.corner_types()}
+        else:
+            return {corner: self.corner(corner) for corner in self.corner_types()}
+
+    def origin(self, **kwargs):
+        """Return footprint origin in world coordinates, as GeoVector.
+           :param  kwargs: additional arguments passed to origin function: dem path and dem.crs for rpcs rasters."""
+        return self.corner('ul', **kwargs)
+
+    def center(self, **kwargs):
+        """Return footprint center in world coordinates, as GeoVector.
+           :param  kwargs: additional arguments passed to corners function: dem path and dem.crs for rpcs rasters."""
+        if self.crs is None:
+            new_raster = self.reproject(dst_crs=WGS84_CRS, rpcs=self.rpcs, **kwargs)
+            image_center = Point(new_raster.width / 2, new_raster.height / 2)
+            return new_raster.to_world(image_center)
         image_center = Point(self.width / 2, self.height / 2)
         return self.to_world(image_center)
 
@@ -1729,8 +1761,11 @@ def bounds(self):
         corners = [self.image_corner(corner) for corner in self.corner_types()]
         return Polygon([[corner.x, corner.y] for corner in corners])
 
-    def _calc_footprint(self):
+    def _calc_footprint(self, **kwargs):
         """Return rectangle in world coordinates, as GeoVector."""
+        if self._crs is None and self.rpcs is not None:
+            return self._footprint_from_rpcs(**kwargs)
+
         corners = [self.corner(corner) for corner in self.corner_types()]
         coords = []
         for corner in corners:
@@ -1742,15 +1777,22 @@ def _calc_footprint(self):
         self._footprint = GeoVector(shp, self.crs)
         return self._footprint
 
-    def footprint(self):
+    def footprint(self, **kwargs):
+        """:param  kwargs: additional arguments passed to footprint function: dem path and dem.crs for rpcs rasters."""
         if self._footprint is not None:
             return self._footprint
-        return self._calc_footprint()
+        return self._calc_footprint(**kwargs)
+
+    def _footprint_from_rpcs(self, **kwargs):
+        """Return raster footprint from rpcs in the crs: EPSG:4326
+           :param  kwargs: additional arguments passed to footprint function: dem path and dem.crs for rpcs rasters."""
+        new_raster = self.reproject(dst_crs=WGS84_CRS, rpcs=self.rpcs, **kwargs)
+        return new_raster.footprint()
 
     def area(self):
         return self.footprint().area
 
-    #  geography:
+    # geography:
     def project(self, dst_crs, resampling):
         """Return reprojected raster."""
 
@@ -1765,7 +1807,7 @@ def to_world(self, shape, dst_crs=None):
         shp = transform(shape, self.crs, dst_crs, dst_affine=self.affine)
         return GeoVector(shp, dst_crs)
 
-    #  array:
+    # array:
     # array ops: bitness conversion, setitem/getitem slices, +-*/.. scalar
     def min(self):
         return self.reduce('min')
@@ -1972,7 +2014,7 @@ def _read_with_mask(raster, masked):
 
     def get_window(self, window, bands=None,
                    xsize=None, ysize=None,
-                   resampling=Resampling.cubic, masked=None, affine=None
+                   resampling=Resampling.cubic, masked=None, affine=None, rpcs=None
                    ):
         """Get window from raster.
 
@@ -1983,6 +2025,8 @@ def get_window(self, window, bands=None,
         :param resampling: which Resampling to use on reading, default Resampling.cubic
         :param masked: if True uses the maks, if False doesn't use the mask, if None looks to see if there is a mask,
                        if mask exists using it, the default None
+        :param affine: Set destination affine otherwise calculate from output window shape
+        :param rpcs: If not none set destination rpcs
         :return: GeoRaster2 of tile
         """
         bands = bands or list(range(1, self.num_bands + 1))
@@ -2004,7 +2048,7 @@ def get_window(self, window, bands=None,
                 array = raster.read(bands, **read_params)
             nodata = 0 if not np.ma.isMaskedArray(array) else None
             affine = affine or self._calculate_new_affine(window, out_shape[2], out_shape[1])
-            raster = self.copy_with(image=array, affine=affine, nodata=nodata)
+            raster = self.copy_with(image=array, affine=affine, nodata=nodata, rpcs=rpcs)
 
             return raster
 

tests/test_georaster.py

@@ -57,6 +57,11 @@
 some_raster_shrunk_mask = some_raster_multiband.copy_with(
     image=np.ma.array(some_raster_multiband.image.data, mask=np.ma.nomask))
 
+footprint_rpcs_raster = Polygon([[111.39116744, 34.92250008],
+                                 [111.44875176, 34.92250008],
+                                 [111.44875176, 34.90383334],
+                                 [111.39116744, 34.90383334]])
+
 
 def test_construction():
     # test image - different formats yield identical rasters:
@@ -248,6 +253,15 @@ def test_resize(raster):
     with pytest.raises(GeoRaster2Error):
         raster.resize(ratio_x=2)
 
+    # test rpcs resize
+    rpcs_raster = GeoRaster2.open("tests/data/raster/grayscale.tif")
+    test_raster = rpcs_raster.resize(ratio=0.5)
+    assert ((test_raster.width == 0.5 * rpcs_raster.width)
+            and (test_raster.rpcs.line_off == 0.5 * rpcs_raster.rpcs.line_off)
+            and (test_raster.rpcs.samp_off == 0.5 * rpcs_raster.rpcs.samp_off)
+            and (test_raster.rpcs.line_scale == 0.5 * rpcs_raster.rpcs.line_scale)
+            and (test_raster.rpcs.samp_scale == 0.5 * rpcs_raster.rpcs.samp_scale))
+
 
 def test_to_pillow_image():
     # without mask:
@@ -424,6 +438,29 @@ def test_corner():
         assert some_raster.corner(corner).almost_equals(GeoVector(expected_corners[corner], some_raster.crs))
         assert some_raster.image_corner(corner).equals_exact(expected_image_corners[corner], tolerance=5e-07)
 
+    coords = footprint_rpcs_raster.exterior.coords.xy
+    rpcs_raster = GeoRaster2.open("tests/data/raster/grayscale.tif")
+
+    expected_image_corners_rpcs = {
+        'ul': Point(0, 0),
+        'ur': Point(rpcs_raster.width, 0),
+        'bl': Point(0, rpcs_raster.height),
+        'br': Point(rpcs_raster.width, rpcs_raster.height)
+    }
+
+    expected_corners_rpcs = {
+        'ul': Point(coords[0][0], coords[1][0]),
+        'ur': Point(coords[0][1], coords[1][1]),
+        'bl': Point(coords[0][3], coords[1][3]),
+        'br': Point(coords[0][2], coords[1][2])
+    }
+
+    for corner in GeoRaster2.corner_types():
+        assert rpcs_raster.corner(corner).almost_equals(GeoVector(expected_corners_rpcs[corner],
+                                                                  WGS84_CRS), decimal=5)
+        assert rpcs_raster.image_corner(corner).equals_exact(expected_image_corners_rpcs[corner],
+                                                             tolerance=5e-07)
+
 
 def test_center():
     ul = some_raster.corner('ul').get_shape(some_raster.crs)
@@ -433,6 +470,16 @@ def test_center():
 
     assert expected_center_vector.equals_exact(some_raster.center(), tolerance=5e-07)
 
+    # test with rpcs
+    rpcs_raster = GeoRaster2.open("tests/data/raster/grayscale.tif")
+
+    ul = rpcs_raster.corner('ul').get_shape(WGS84_CRS)
+    br = rpcs_raster.corner('br').get_shape(WGS84_CRS)
+    expected_center = Point((ul.x + br.x) / 2, (ul.y + br.y) / 2)
+    expected_center_vector = GeoVector(expected_center, WGS84_CRS)
+
+    assert expected_center_vector.equals_exact(rpcs_raster.center(), tolerance=5e-07)
+
 
 def test_bounds():
     expected = Polygon([[0, 0],
@@ -452,6 +499,22 @@ def test_footprint():
 
     assert some_raster.footprint().almost_equals(expected)
 
+    rpcs_raster = GeoRaster2.open("tests/data/raster/grayscale.tif")
+
+    expected = GeoVector(footprint_rpcs_raster, WGS84_CRS)
+    footprint = rpcs_raster.footprint()
+
+    assert footprint.almost_equals(expected, decimal=5)
+
+    # test raster footprint by passing fake DEM
+    dem = GeoRaster2.open("tests/data/raster/dem_grayscale_raster.tif")
+
+    kwargs = {'RPC_DEM': dem.source_file,
+              'RPC_DEM_SRS': WGS84_CRS  # the dem shall be referred to the RPC EPSG
+              }
+
+    assert footprint != rpcs_raster.footprint(**kwargs)
+
 
 def test_area():
     scale = 2