From 766e122f981e895991fc16ac81ad1ee60481bfd7 Mon Sep 17 00:00:00 2001 From: Simone Date: Tue, 26 Apr 2022 12:31:02 +0200 Subject: [PATCH 1/2] Fix copy_with for rpcs, add rpcs for footprint, corners, corner, center, origin methods --- telluric/georaster.py | 61 +++++++++++++++------ tests/data/raster/dem_grayscale_raster.tif | Bin 0 -> 2348 bytes tests/test_georaster.py | 54 ++++++++++++++++++ 3 files changed, 97 insertions(+), 18 deletions(-) create mode 100644 tests/data/raster/dem_grayscale_raster.tif diff --git a/telluric/georaster.py b/telluric/georaster.py index 0a65e1c..1cde588 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) @@ -1707,20 +1707,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 +1744,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 +1760,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 +1790,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') diff --git a/tests/data/raster/dem_grayscale_raster.tif b/tests/data/raster/dem_grayscale_raster.tif new file mode 100644 index 0000000000000000000000000000000000000000..9b39aee444ae05d297e7fdbe7f32854a861f0788 GIT binary patch literal 2348 zcmebD)M7Zo#=s!Rz`)4Nz{tSB5DCPLP_{gf%>-pT0NKn?HY5gK zt$^$e?K})jKz0?7y|JB%fdk0y0u5uzu|-sZH( zkfTOPTDaR~(FTnpMdp0n#~*LvxYFjUI>lAx(@LMTNd4(P$(MFcnUrWe?VQ@@%r$+r zN(+3|gKzOXt8`ux_%isHs9v}AvcOfQT(v^8gI9;F+A`~C(RRPp*I#ery4rR&bxWwu zx7EIBH}kj0WMA7ob<)G)ZTIxPXRq!1n{*&nKl~2wv!BgJ5m;t{rvMTo@?#CORuDw{95am9=ZN{PVue1(j? z&7<$}J+D0cr0{3-KXLu;-Ombl+49v3&%XS+WY?Z~$BVb0d;R_QKECViXP1quv81WQQ= zq*4-B-L(!_SmQ1wA@$oLpc+D@By!!g3#cAecU=OqVRaW1u4getf@t6Pr literal 0 HcmV?d00001 diff --git a/tests/test_georaster.py b/tests/test_georaster.py index b48999a..a6a4fed 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: @@ -424,6 +429,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 +461,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 +490,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 From d3c232205a5e74cafbd6432c5e0bfd116ca9b567 Mon Sep 17 00:00:00 2001 From: Simone Date: Tue, 26 Apr 2022 12:40:44 +0200 Subject: [PATCH 2/2] Add resize method for rpc rasters and unit test --- telluric/georaster.py | 25 ++++++++++++++++++++++--- tests/test_georaster.py | 9 +++++++++ 2 files changed, 31 insertions(+), 3 deletions(-) diff --git a/telluric/georaster.py b/telluric/georaster.py index 1cde588..337308c 100644 --- a/telluric/georaster.py +++ b/telluric/georaster.py @@ -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) @@ -1997,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. @@ -2008,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)) @@ -2029,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/test_georaster.py b/tests/test_georaster.py index a6a4fed..61082ac 100644 --- a/tests/test_georaster.py +++ b/tests/test_georaster.py @@ -253,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: