diff --git a/telluric/georaster.py b/telluric/georaster.py index 0a65e1c..337308c 100644 --- a/telluric/georaster.py +++ b/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 diff --git a/tests/data/raster/dem_grayscale_raster.tif b/tests/data/raster/dem_grayscale_raster.tif new file mode 100644 index 0000000..9b39aee Binary files /dev/null and b/tests/data/raster/dem_grayscale_raster.tif differ diff --git a/tests/test_georaster.py b/tests/test_georaster.py index b48999a..61082ac 100644 --- a/tests/test_georaster.py +++ b/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