I need help on coordinate transformation issue that I have been trying to fix for quite some days now. I want to calculate snow depth by differencing snow-on and Snow-off lidar DEM. When I did the differencing, I got about -18m as the snow depth. This value is not correct. The issue is that the DEMS are in different coordinate systems. When I checked the crs of the DEMs, the result I got is given below: snow-on: +proj=utm +zone=11 +datum=WGS84 +units=m +vunits=m +no_defs snow-off: +proj=utm +zone=11 +ellps=GRS80 +units=m +no_defs
The snow-on DEM acquired using a helipod is referenced to the WGS84 ellipsoid. The srs info above doesn't include the datum info for the snow-off DEM but according to QSI, it is in Geoid. The information is provided below: Coordinate System: UTM Zone 11 North Horizontal Datum: NAD83 (2011) Vertical Datum: NAVD88 (Geoid 12b)
Now the -18m value, though not right, makes sense because the separation between geoid and WGS 84 datum is about 17m for the UTM zone of my region of interest.
Now I know what to do in this case is perform coordinate transformation. I need to transform either of the DEMs to the same vertical datum as the other. I transformed the QSI data (snow free) to WGS 84 using gdalwarp. My first try is to specify the source and target coordinate in proj format and include the reference to the grid file for g2012a. The command that I used is given below:
!gdalwarp -overwrite ./data/Mores_Bare_Earth_DEMs_mosaic.tif ./data/Mores_Bare_Earth_DEMs_mosaic_WGS84_proj.tif -s_srs "+proj=utm +zone=11 +ellps=GRS80 +units=m +no_defs +geoidgrids=/usr/share/proj/g2012a_conus.gtx" -t_srs "+proj=utm +zone=11 +datum=WGS84 +units=m +no_defs"
Then I performed the differencing again. This time, I got about -2m which is still not right. My diagnosis is that I am using the wrong grid file.. so I tried a couple other gtx files such as egm96_15.gtx but still didn't get a good result. Then I specify the source and target coordinates in EPSG format. I used EPSG:26911+5703 to specify the horizontal and vertical crs of the source and 32611 as the target (WGS84).
!gdalwarp -overwrite ./data/Mores_Bare_Earth_DEMs_mosaic.tif ./data/Mores_Bare_Earth_DEMs_mosaic_WGS84_epsgv2.tif -s_srs EPSG:26911+5703 -t_srs EPSG:4979
The result I got here is similar to what I got when I differenced before transforming, suggesting that transformation didn't work.
I have tried several other options too. For instance, I transformed the helipod data to geoid but still no luck. At this time, I have exhausted all possible options I could think about.
