Skip to content

clip_to_poly fails for multiband raster

Hi,

while using arosics, I ran into a small bug when trying to clip the resulting GeoArray to its footprint. Here is an example with a multiband raster:

from geoarray import GeoArray
import numpy as np
import rasterio

# just some random data
crs = rasterio.CRS.from_epsg(32632)
transform = rasterio.Affine(100.0, 0.0, 797601, 0.0, -100.0, 5827090)
nodata = 0

array_data = np.random.randn(100, 200, 3)
# force some no data into the first column (not sure how it behaves if there is nothing to clip...)
array_data[:,0,:] = nodata

geo_array = GeoArray(array_data, geotransform=transform.to_gdal(), projection=crs.to_wkt(), nodata=nodata)

# fails with AssertionError
geo_array.clip_to_footprint()
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
Cell In[52], line 1
----> 1 geo_array.clip_to_footprint()

File ~/Documents/misc/2023-02-09_bug_geoarray/./geoarray-main/geoarray/baseclasses.py:1601, in GeoArray.clip_to_footprint(self)
   1599 def clip_to_footprint(self):
   1600     """Clip the GeoArray instance to the outer bounds of the actual footprint."""
-> 1601     self.clip_to_poly(self.footprint_poly)

File ~/Documents/misc/2023-02-09_bug_geoarray/./geoarray-main/geoarray/baseclasses.py:1612, in GeoArray.clip_to_poly(self, poly)
   1609 self.arr, self.gt, self.projection = self.get_mapPos(mapBounds=poly.bounds)
   1610 self.mask_nodata.arr, self.mask_nodata.gt, self.mask_nodata.projection = \
   1611     self.mask_nodata.get_mapPos(mapBounds=poly.bounds, mapBounds_prj=self.prj)
-> 1612 assert self.shape == self.mask_nodata.shape
   1614 if self._mask_baddata is not None:
   1615     self.mask_baddata.arr, self.mask_baddata.gt, self.mask_baddata.projection = \
   1616         self.mask_baddata.get_mapPos(mapBounds=poly.bounds)

AssertionError:

I guess this could easily be solved by changing assert self.shape == self.mask_nodata.shape to assert self.shape[0:2] == self.mask_nodata.shape to limit the shape check to first to dimensions. Likewise for assert self.shape == self.mask_baddata.shape a couple lines later.

I just recently started to use geoarray, so I a not entirely sure if I missed something here. But as far as I can see, the mask is always 2D (height, width) and not (height, width, n_channels), so checking against the first two dimension should suffice? BTW, Everything works as expected once I drop the third dimension from the numpy array.

Best, geopanda1

Edited by geopanda1