potential fix for resamplemethod in reproject #3542
Conversation
|
The fix works but I had to fall back to old (bad) behaviour in cases at the edge of valid target crs bounds. My method to compute target cellsize resulted in an error there. |
| val targetCellSizeInSrcCRS = | ||
| try{ | ||
| rasterExtent.reproject(dest,src).cellSize | ||
| }catch { | ||
| case e:Exception => //reprojection errors happen when going outside valid area of projection, need a better fix here | ||
| raster.rasterExtent.cellSize | ||
| } |
There was a problem hiding this comment.
are there any tests that cover this case?
There was a problem hiding this comment.
Not yet, but I now had some time to dig deeper.
The original stack trace was:
at geotrellis.raster.GridExtent.<init>(GridExtent.scala:46)
at geotrellis.raster.GridExtent.<init>(GridExtent.scala:58)
at geotrellis.raster.reproject.ReprojectRasterExtent$.apply(ReprojectRasterExtent.scala:71)
at geotrellis.raster.GridExtent$gridExtentMethods.$anonfun$reproject$2(GridExtent.scala:432)
at scala.Option.getOrElse(Option.scala:189)
at geotrellis.raster.GridExtent$gridExtentMethods.reproject(GridExtent.scala:432)
at geotrellis.raster.GridExtent$gridExtentMethods.reproject(GridExtent.scala:435)
at org.openeo.geotrellis.reproject.RasterRegionReproject$$anon$2.reprojectToBuffer(RasterRegionReproject.scala:268)
at org.openeo.geotrellis.reproject.RasterRegionReproject$$anon$2.regionReproject(RasterRegionReproject.scala:229)
at org.openeo.geotrellis.reproject.TileRDDReproject$.createCombiner$1(TileRDDReproject.scala:192)
at org.openeo.geotrellis.reproject.TileRDDReproject$.$anonfun$apply$6(TileRDDReproject.scala:206)
and got triggered by a more complicated test case here:
https://github.com/Open-EO/openeo-geotrellis-extensions/blob/106784e1cb83a87f1bacef3beb1b9df77c5ecb65/openeo-geotrellis/src/test/scala/org/openeo/geotrellis/OpenEOProcessesSpec.scala#L582
I could reproduce it in a simpler test, getting a similar but not exactly the same stack trace:
it("should (approximately) match a GDAL average interpolation on nlcd tile at edge of valid bounds") {
val tile = this.createTile(Array(1,4,1,4,1,4,1,4),4,2)
val raster = Raster(tile,Extent(-40.50000387499999, -89.99999, 40.499993875000015, -51.749994249999986))
val tempTiff = File.createTempFile("toReproject",".tif")
GeoTiff(raster, LatLng).write(tempTiff.getPath)
val tiffRe = RasterExtent(Extent(-4452779.631730943, -6.325724963354129E7, 0.0, -5.839130735403811E7), CellSize(15000.0,15000.0))
val alignment = TargetRegion(tiffRe)
val rs = GeoTiffReprojectRasterSource(tempTiff.getPath,WebMercator, alignment, Average)
val reprojectedFromRS = rs.read().get.tile.band(0)
val reprojected = raster.reproject(LatLng,WebMercator, Options(method = Average, errorThreshold = 0.0,targetCellSize = Some(CellSize(1000,1000))))
val refTile = this.createTile(Array(2.5,2.5),2226, 10363)
//assertEqual(refTile,reprojected.tile,0.1)
assertEqual(reprojectedFromRS,reprojected.tile,0.1)
}
is giving:
invalid rows: 0
geotrellis.raster.GeoAttrsError: invalid rows: 0
at geotrellis.raster.GridExtent.<init>(GridExtent.scala:46)
at geotrellis.raster.GridExtent.<init>(GridExtent.scala:58)
at geotrellis.raster.reproject.ReprojectRasterExtent$.apply(ReprojectRasterExtent.scala:71)
at geotrellis.raster.reproject.ReprojectRasterExtent$.apply(ReprojectRasterExtent.scala:76)
at geotrellis.raster.geotiff.GeoTiffReprojectRasterSource.closestTiffOverview$lzycompute(GeoTiffReprojectRasterSource.scala:86)
at geotrellis.raster.geotiff.GeoTiffReprojectRasterSource.closestTiffOverview(GeoTiffReprojectRasterSource.scala:81)
at geotrellis.raster.geotiff.GeoTiffReprojectRasterSource.readBounds(GeoTiffReprojectRasterSource.scala:110)
at geotrellis.raster.geotiff.GeoTiffReprojectRasterSource.read(GeoTiffReprojectRasterSource.scala:98)
at geotrellis.raster.geotiff.GeoTiffReprojectRasterSource.read(GeoTiffReprojectRasterSource.scala:94)
at geotrellis.raster.RasterSource.read(RasterSource.scala:134)
at geotrellis.raster.reproject.ReprojectSpec.$anonfun$new$5(ReprojectSpec.scala:92)
What's happening here is that we are going across valid bounds of webmercator. The root cause is then in these lines:
https://github.com/VitoTAP/geotrellis/blob/b071b333581ed7206ad14a06b1b70cfde6fb65a7/raster/src/main/scala/geotrellis/raster/reproject/ReprojectRasterExtent.scala#L63
The 'toLong' will round down to zero, which is not accepted later on in the code.
When that case happens, I think it also doesn't really matter any more what cellsize is used, because we warp into an area with cols or rows == 0.
There was a problem hiding this comment.
ic ic; the thing is that this try / catch is a bit sketchy 🤔 we definitely need to find a better way addressing it / or at least make exception more specific.
There was a problem hiding this comment.
Here's a case that really isolates and exactly reproduces the issue. Probably these extents are fairly problematic to begin with:
it("should (approximately) match a GDAL average interpolation on nlcd tile at edge of valid bounds") {
val tile = this.createNoData(288,272,ByteConstantNoDataCellType)
val raster = Raster(tile,Extent(-40.50000387499999, -89.99999, 40.499993875000015, -51.749994249999986))
val tiffRe = RasterExtent(Extent(-4452779.631730943, -6.325724963354129E7, 0.0, -5.839130735403811E7), CellSize(14744.303416327626,16112.391653984045))
val rr = implicitly[RasterRegionReproject[Tile]]
val reprojectedRasterRegion: Raster[Tile] = rr.regionReproject(
raster,
LatLng,
WebMercator,
tiffRe,
tiffRe.extent.toPolygon(),
Average,
0.2
)
//succesfull if it doesn't throw an exception
}There was a problem hiding this comment.
Hm ok I see what you need it for.
There was a problem hiding this comment.
I compared the behavior with the GDAL behavior. /tmp/raster-lat-lng.tif is GeoTiff(raster, LatLng).write("/tmp/raster-lat-lng.tif")
$ gdalwarp -t_srs EPSG:3857 -tr 14744.303416327626 16112.391653984045 -te -4452779.631730943 -6.325724963354129E7 0.0 -5.839130735403811E7 /tmp/raster-lat-lng.tif /tmp/raster-merc.tif`
$ gdalwarp -t_srs EPSG:4326 /tmp/raster-merc.tif /tmp/raster-merc-lat.tif
Creating output file that is 427P x 0L.
ERROR 1: Attempt to create 427x0 dataset is illegal,sizes must be larger than zero.This is exactly what happens in our reproject function: newCols: 427 newRows: 0 which throws geotrellis.raster.GeoAttrsError: invalid rows: 0
There was a problem hiding this comment.
I think the change is good, sadly I don't think we can merge this workaround with throws; it simply uses incorrect resolutions as an input for the functions.
Mb it is acceptable (it worked somehow up until now) but not sure that's what we really need to do here 🤔 The error thrown seems to be legit.
| case e:Exception => //reprojection errors happen when going outside valid area of projection, need a better fix here | ||
| raster.rasterExtent.cellSize | ||
| } | ||
| val resampler = Resample(resampleMethod, raster.tile, raster.extent, targetCellSizeInSrcCRS) |
There was a problem hiding this comment.
I see, the issue is that the targetCellSize was / is in the wrong CRS.
|
|
||
| val targetCellSizeInSrcCRS = | ||
| try{ | ||
| rasterExtent.reproject(dest,src).cellSize |
There was a problem hiding this comment.
I wonder if we have a broken function.
#3541
Overview
Proposed fix for a very important bug in reprojecting to lower resolution.
I think it's better (more correct) than previous behaviour, but not sure yet if it is 100% accurate.
It seems to me that the downsampling now happens in the original projection system, and then it is warping to the new CRS. This will work fine if the low resolution pixel is rectangular in both projection systems.
Checklist
docsguides update, if necessaryDemo
Optional. Screenshots/REPL
Notes
Optional. Ancillary topics, caveats, alternative strategies that didn't work out, anything else.
Closes #XXX