GlacioHack
Re-projecting with multi-step proj string pipeline
Dear xdem community,
Is there a way to use a multi-step proj pipeline string to reproject a raster?
I have two DEMs with the following metadata printed by GDAL.info:
DEM 1
Driver: GTiff/GeoTIFF
Files: /scratch-3/cogier/data/DEM_input/grid_bedrock_10m.tif
Size is 300, 149
Coordinate System is:
PROJCRS["CH1903+ / LV95",
BASEGEOGCRS["CH1903+",
DATUM["CH1903+",
ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4150]],
CONVERSION["Swiss Oblique Mercator 1995",
METHOD["Hotine Oblique Mercator (variant B)",
ID["EPSG",9815]],
PARAMETER["Latitude of projection centre",46.9524055555556,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8811]],
PARAMETER["Longitude of projection centre",7.43958333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8812]],
PARAMETER["Azimuth of initial line",90,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8813]],
PARAMETER["Angle from Rectified to Skew Grid",90,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8814]],
PARAMETER["Scale factor on initial line",1,
SCALEUNIT["unity",1],
ID["EPSG",8815]],
PARAMETER["Easting at projection centre",2600000,
LENGTHUNIT["metre",1],
ID["EPSG",8816]],
PARAMETER["Northing at projection centre",1200000,
LENGTHUNIT["metre",1],
ID["EPSG",8817]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping (large and medium scale)."],
AREA["Liechtenstein; Switzerland."],
BBOX[45.82,5.96,47.81,10.49]],
ID["EPSG",2056]]
Data axis to CRS axis mapping: 1,2
Origin = (2511507.080000000074506,976908.109999999986030)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 2511507.080, 976908.110) ( 6d19' 8.15"E, 44d56'23.21"N)
Lower Left ( 2511507.080, 975418.110) ( 6d19' 9.13"E, 44d55'34.97"N)
Upper Right ( 2514507.080, 976908.110) ( 6d21'24.91"E, 44d56'24.58"N)
Lower Right ( 2514507.080, 975418.110) ( 6d21'25.84"E, 44d55'36.34"N)
Center ( 2513007.080, 976163.110) ( 6d20'17.01"E, 44d55'59.78"N)
Band 1 Block=300x6 Type=Float32, ColorInterp=Gray
DEM 2
Driver: GTiff/GeoTIFF
Files: /scratch-3/cogier/data/DEM_input/Bedrock.tif
Size is 289, 155
Coordinate System is:
PROJCRS["RGF93 v1 / Lambert-93",
BASEGEOGCRS["RGF93 v1",
DATUM["Reseau Geodesique Francais 1993 v1",
ELLIPSOID["GRS 1980",6378137,298.257222101004,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4171]],
CONVERSION["Lambert Conic Conformal (2SP)",
METHOD["Lambert Conic Conformal (2SP)",
ID["EPSG",9802]],
PARAMETER["Latitude of false origin",46.5,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8821]],
PARAMETER["Longitude of false origin",3,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8822]],
PARAMETER["Latitude of 1st standard parallel",49,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8823]],
PARAMETER["Latitude of 2nd standard parallel",44,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8824]],
PARAMETER["Easting at false origin",700000,
LENGTHUNIT["metre",1],
ID["EPSG",8826]],
PARAMETER["Northing at false origin",6600000,
LENGTHUNIT["metre",1],
ID["EPSG",8827]]],
CS[Cartesian,2],
AXIS["easting",east,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["northing",north,
ORDER[2],
LENGTHUNIT["metre",1]],
ID["EPSG",2154]]
Data axis to CRS axis mapping: 1,2
Origin = (961677.866099768900312,6432278.289214908145368)
Pixel Size = (10.683080483055157,-10.683080483053720)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 961677.866, 6432278.289) ( 6d19' 5.69"E, 44d56'24.78"N)
Lower Left ( 961677.866, 6430622.412) ( 6d19' 2.52"E, 44d55'31.16"N)
Upper Right ( 964765.276, 6432278.289) ( 6d21'26.47"E, 44d56'20.56"N)
Lower Right ( 964765.276, 6430622.412) ( 6d21'23.25"E, 44d55'26.93"N)
Center ( 963221.571, 6431450.350) ( 6d20'14.48"E, 44d55'55.86"N)
Band 1 Block=128x128 Type=Float32, ColorInterp=Gray
Min=2402.215 Max=3221.648
Minimum=2402.215, Maximum=3221.648, Mean=2700.464, StdDev=175.148
NoData Value=-9999
Metadata:
RepresentationType=ATHEMATIC
STATISTICS_COVARIANCES=30676.71358956566
STATISTICS_MAXIMUM=3221.6481933594
STATISTICS_MEAN=2700.4638964504
STATISTICS_MINIMUM=2402.2153320313
STATISTICS_SKIPFACTORX=1
STATISTICS_SKIPFACTORY=1
STATISTICS_STDDEV=175.14769079142
I reprojected DEM 1 into the CRS of DEM 2 as follows:
import xdem
dem1 = xdem.DEM('grid_bedrock_10m.tif')
dem2 = xdem.DEM('Bedrock.tif')
dem1_reproj = dem1.reproject(dem2)
dem1_reproj.save('grid_bedrock_10m_reproj.tif')
Unfortunately, this transformation does not put the DEM in the correct place: it is slightly shifted.
However, when I load the DEM 1 and DEM 2 into QGIS, QGIS performs an on-the-fly re-projection that projects DEM 1 into the same CRS as DEM 2. QGIS tells me that it uses the following multi-step proj pipeline for the on-the-fly re-projection:
+proj=pipeline +step +inv +proj=somerc +lat_0=46.9524055555556 +lon_0=7.43958333333333 +k_0=1 +x_0=2600000 +y_0=1200000 +ellps=bessel +step +proj=push +v_3 +step +proj=cart +ellps=bessel +step +proj=helmert +x=674.374 +y=15.056 +z=405.346 +step +inv +proj=cart +ellps=GRS80 +step +proj=pop +v_3 +step +proj=lcc +lat_0=46.5 +lon_0=3 +lat_1=49 +lat_2=44 +x_0=700000 +y_0=6600000 +ellps=GRS80
How can I use this multi-step proj string pipeline in xdem?
Other ideas on how to repoject are also welcome
I can share the DEMs privately if needed.
Thanks a lot :)
Thanks @christopheogier for the issue.
The DEM is shifted vertically? If yes, the source of the problem is likely that Rasterio that doesn't support 3D reprojection, which is out of our hands: https://github.com/rasterio/rasterio/issues/2433.
The temporary solution via xDEM is to isolate your 2D horizontal CRSs and vertical CRSs (.crs and .vcrs) for each of your DEMs, and use the reprojection steps sequentially reproject(crs1, crs2) then set_vcrs(vcrs1) and to_vcrs(vcrs2), like shown here: https://xdem.readthedocs.io/en/stable/vertical_ref.html. It might not be a perfect match with QGIS/GDAL (not always equivalent to do the 2 steps in a row or having a single one), but should be much closer.
If you share your .crs/ .transform/.shape attributes for both DEMs (we don't need the real data), we could check if it works on placeholder data entirely in Python (by comparing directly with gdal.Reproject() like here: https://github.com/GlacioHack/xdem-data/blob/98004a09f84def4c78b253d41b212baca2b3cccb/generate_ground_truth_gdal.py#L93).
While ultimately I'd love for that rasterio issue to be closed, I learned relatively recently that you can encode the proj string in a GDAL VRT. If you do that, then opening the VRT with rasterio, GDAL will apply the pipeline on loading.
I have only done a few tests with some global DEMs, and the approach seems to work but has required some tinkering to get the bounds right (https://github.com/OSGeo/gdal/issues/11610)
We have some documentation on creating such a VRT here https://uw-cryo.github.io/3D_CRS_Transformation_Resources/cop30#encoding-transform-in-a-gdal-vrt
But since I'd like to reorganize that documentation in the near future (also contributions / suggestions for organization would be welcome) here is a link to what such a VRT looks like https://github.com/uw-cryo/3D_CRS_Transformation_Resources/blob/acf0787e7dd5a45c21eee33dba8f5b6d403285c4/globaldems/COP30_hh_7912.vrt#L34
Thanks for the quick answer, @rhugonnet and @scottyhq !
The DEM is not shifted vertically, only horizontally (both in x and y). So I guess the vcrs is not involved here. That said, and because we are talking about 3D reprojection, I remark that the ellipsoids are not similar (which could be a first hint on what's going on?):
DEM1
ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
LENGTHUNIT["metre",1]]]
DEM2
ELLIPSOID["GRS 1980",6378137,298.257222101004,
LENGTHUNIT["metre",1]]]
@rhugonnet: Here are the .crs / .transform and .shape for both DEMs (respectively)
DEM1
CRS.from_epsg(2056)
Affine(10.0, 0.0, 2511507.08, 0.0, -10.0, 976908.11)
(149, 300)
DEM2
CRS.from_epsg(2154)
Affine(10.683080483055157, 0.0, 961677.8660997689, 0.0, -10.68308048305372, 6432278.289214908)
(155, 289)
@scottyhq I will look at your link later, Thanks so much for your interest and the docs on this :)