Support for 3D transformations and PROJ pipelines #2433
Comments
|
One thing @scottyhq mentioned in the previous discussion is that |
|
Yeah, we really need to prioritize #1990. It's on my mind for 1.3 and I may even get to do something about it next week. |
|
I think this is a duplicate of #2354 after all. Or not 😄 |
|
@scottyhq if I read https://github.com/OSGeo/gdal/blob/master/apps/gdalwarp_lib.cpp#L613-L614 correctly, vertical grid shift with GDAL 3.4 might be "just" this: My work on #1990 is blocked by missing features of the GDALWarp API, so I'm looking at other ways around the problem. |
|
Thanks for looking into this @sgillies! I wonder if that could be set as a default somehow to mirror gdalwarp behavior? "Starting with GDAL 2.2, if the SRS has an explicit vertical datum that points to a PROJ.4 geoidgrids, and the input dataset is a single band dataset, a vertical correction will be applied to the values of the dataset." I can confirm But... The result is not the same result as gdalwarp (which seems to more accurately reflect the shift grid values from https://cdn.proj.org/us_nga_egm08_25.tif). I think this is best illustrated with a figure:
|
|
@scottyhq are you sure you're using the same versions of GDAL, PROJ, its grids, on each side of the comparison? It's easy to get into a place where we have software and data discrepancies. |
As far as I can tell, yes. I saw #2349, and wasn't quite sure how to verify command line gdal and dependencies match the rasterio wheels, so I'm going off https://github.com/rasterio/rasterio-wheels/blob/master/env_vars.sh: I've gone ahead and documented things in more detail here using https://github.com/scottyhq/riodatum. Since this is all public data you can run the exact same environment and generate the above figure via One thing that jumps out at me is what looks like a vertical artifact at -108.5 longitude, and I notice the rasterio log outputs the following: |
After a bit more exploring, this artifact appears to be coming from the
|
With
|
|
I observe the same behavior documented above with rasterio=1.3a4, but it seems that brings GDAL3.4.2 and PROJ8.2.1 whereas the GDAL3.4.2 docker containers or GDAL from conda-forge bring PROJ 9.0.0. I also tried passing the pipeline directly to reproject in 1.3a4, but the result is the same as documented above: proj_pipeline="proj=pipeline step proj=axisswap order=2,1 step proj=unitconvert xy_in=deg xy_out=rad step proj=vgridshift grids=us_nga_egm08_25.tif multiplier=1 step proj=unitconvert xy_in=rad xy_out=deg step proj=axisswap order=2,1"
with rasterio.open(url) as src:
# original metadata unaware of vertical reference
src_crs3D = rasterio.crs.CRS.from_epsg('9518')
dst_crs = rasterio.crs.CRS.from_epsg('7661')
transform, width, height = calculate_default_transform(src_crs3D,
dst_crs,
src.width,
src.height,
*src.bounds,
coordinate_operation=proj_pipeline)
kwargs = src.meta.copy()
kwargs.update({
'crs': dst_crs,
'transform': transform,
'width': width,
'height': height
})
with rasterio.open('epsg7661_rio13.tif', 'w', **kwargs) as dst:
reproject(
rasterio.band(src, src.indexes),
rasterio.band(dst, dst.indexes),
resampling=Resampling.nearest,
apply_vertical_shift=True,
coordinate_operation=proj_pipeline,
warp_mem_limit=10,
) |
Need to install with |
@snowman2 - I'm fairly convinced there is something with Nevertheless, I tried the following with rasterio 1.3b1 (perhaps this deserves a separate issue or discussion thread): Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File "/usr/local/lib/python3.8/dist-packages/rasterio/__init__.py", line 13, in <module>
from rasterio.crs import CRS
File "rasterio/crs.pyx", line 1, in init rasterio.crs
File "rasterio/_base.pyx", line 31, in init rasterio._base
File "/usr/local/lib/python3.8/dist-packages/rasterio/transform.py", line 14, in <module>
from rasterio.env import env_ctx_if_needed
File "/usr/local/lib/python3.8/dist-packages/rasterio/env.py", line 16, in <module>
from rasterio._env import (
File "rasterio/_env.pyx", line 1, in init rasterio._env
ImportError: /usr/local/lib/python3.8/dist-packages/rasterio/_filepath.cpython-38-x86_64-linux-gnu.so: undefined symbol: VSIAllocFilesystemPluginCallbacksStruct |
I am seeing similar, slightly different, errors with a GDAL 3.5 test here: https://github.com/rasterio/rasterio/runs/6538351664 Do you see a similar issue with GDAL 3.4? |
|
@sgillies @snowman2 I reran this analysis with the latest matched GDAL and rasterio and am pretty confident it is a bug. Interestingly, I'm observing that reducing
So maybe breaking the image up by |
|
I wonder if there is any connection with this fix #2512? |
|
Wonder if this is related? conda-forge/rasterio-feedstock#239 |
|
This thread reminds me of our discussion with @snowman2 here: pyproj4/pyproj#761. Is it really useful to have this implemented directly in (btw, hi @scottyhq! I just started at UW) |
|
@snowman @sgillies apologies for letting this issue drift a bit with intermediate updates, but as documented above we've seen that PROJ pipelines are supported in rasterio>=1.3, and rasterio does make use shift grids as long as the However, for outputs requiring a vertical shift grid rasterio is giving incorrect outputs apparently due to different sampling of the shift grid compared to gdalwarp. I keep seeing this for all recent versions of rasterio no matter how it's installed. I can open a separate issue for this specific bug if that would help organize things. The discrepancy is best illustrated using a coarse grid (like EGM96 for SRTM data below):
code for above figure: https://gist.github.com/scottyhq/af8ee3c440cbe0a3a5ffaa538643b104 It's hard to debug this further because I'm not sure how inspect intermediate data structures or follow the difference in code paths between GDAL and rasterio. It was suggested earlier that switching _reproject to GDALWarp internally might fix this and the previous blocker seems fixed on the GDAL side to continue pursuing that (#2437 (comment))? I'd be happy to test if possible. The only differences standing out to me in logs at this point are:
|
Expected behavior and actual behavior.
Currently it seems
rasterio.warp.reprojectdoes not perform 3D reprojections, for example converting elevation rasters from geoid to ellipsoid heights (CompoundCRS -> CompoundCRS).Below is an example for a relatively straightforward transform (
projinfo -s EPSG:9518 -t EPSG:7661 -o PROJ --hide-ballpark --spatial-test intersects). Looking at DEBUG logs rasterio does invoke PROJ which retrieves the correct vertical shift grid, but the grid is not applied to raster values as expected.The PROJ CRS conversion guide shows how tricky these pipelines can get. But often for modern datasets such as Copernicus DEM, it's just a matter of applying a vertical shift grid.
Previous discussion with @snowman2 (corteva/rioxarray#494) suggests this could be fixed by #1990?
Steps to reproduce the problem.
DEBUG log
Click to show log
Operating system
Ubuntu 18.04.6 LTS
Rasterio version and provenance
rasterio 1.2.10 (conda-forge). also tested via pypi install
The text was updated successfully, but these errors were encountered: