GlacioHack / geoutils

Analysis of georeferenced rasters and vectors
https://geoutils.readthedocs.io/
BSD 3-Clause "New" or "Revised" License
85 stars 19 forks source link

`Raster.reproject()` not working when only providing `dst_bounds` #147

Closed erikmannerfelt closed 3 years ago

erikmannerfelt commented 3 years ago

Minimal example:

In [1]: import geoutils as gu

In [2]: r = gu.Raster(gu.datasets.get_path("landsat_B4"))

In [3]: r.bounds
Out[3]: BoundingBox(left=478000.0, bottom=3088490.0, right=502000.0, top=3108140.0)

In [4]: import rasterio as rio

In [5]: new_bounds = rio.coords.BoundingBox(left=478000.0, bottom=3088490.0, right=512000.0, top=3108140.0)

In [6]: r.reproject(dst_bounds=new_bounds)
/home/erik/Projects/ETH/GlacioHack/GeoUtils/geoutils/georaster.py:717: UserWarning: Output projection, bounds and size are identical -> return self (not a copy!)
  warnings.warn("Output projection, bounds and size are identical -> return self (not a copy!)")
Out[6]: 
geoutils.georaster.Raster(BoundingBox(left=478000.0, bottom=3088490.0, right=502000.0, top=3108140.0), 1, EPSG:32645, <built-in method dataset_mask of DatasetReader object at 0x7fe3fbec6840>, GTiff, ('uint8',), 655, (1,), /home/erik/Projects/ETH/GlacioHack/GeoUtils/geoutils/datasets/LE71400412000304SGS00_B4_crop.TIF, None, (30.0, 30.0), (655, 800), | 30.00, 0.00, 478000.00|
| 0.00,-30.00, 3108140.00|
| 0.00, 0.00, 1.00|, 800)

See that the "copy" warning turns up and the returned bounds are not equal to the dst_bounds.

This is because dst_transform will only be created (or rather not None) if dst_res is provided:

In [7]: r.reproject(dst_bounds=new_bounds, dst_res=r.res)
Out[7]: 
geoutils.georaster.Raster(BoundingBox(left=478000.0, bottom=3088490.0, right=512020.0, top=3108140.0), 1, EPSG:32645, <built-in method dataset_mask of DatasetReader object at 0x7fe3fbec6b40>, GTiff, ('uint8',), 655, (1,), /vsimem/0306adbb-68df-4ba2-9e48-0219ea1c8d46/0306adbb-68df-4ba2-9e48-0219ea1c8d46.tif, None, (30.0, 30.0), (655, 1134), | 30.00, 0.00, 478000.00|
| 0.00,-30.00, 3108140.00|
| 0.00, 0.00, 1.00|, 1134)

Below is the line that messes this up:

https://github.com/GlacioHack/GeoUtils/blob/2c38c0be21c9526ce8cbf5ef2495bdacc91feb4b/geoutils/georaster.py#L666

erikmannerfelt commented 3 years ago

Wait, I already reported this in #140! I feel like I'm going crazy!

adehecq commented 3 years ago

Haha, yes I was wondering why a new issue!