Search code examples
transformgdalgeo

gdalwarp messing up GSD

I downloaded the DHM25 raster: https://www.swisstopo.admin.ch/en/geodata/height/dhm25.html

It comes with a dhm25_grid_raster.asc file that I would like to warp to EPSG:4326 (WGS84). That site lists the coordinate reference as "LV03 LN02" and the documentation that comes with the data lists it as "CH1903" which I believe is either EPSG:21781 or EPSG:5728.

Question 1: Does anyone know if it should be EPSG:21781 or EPSG:5728?

When I run the following code, it seems to work converting it to a .tif:

drv = gdal.GetDriverByName('GTiff')
ds_in = gdal.Open('dhm25_grid_raster.asc')
ds_out = drv.CreateCopy('dhm25_grid_raster.tif', ds_in)
srs = osr.SpatialReference()
srs.ImportFromEPSG(21781) #5728?
ds_out.SetProjection(srs.ExportToWkt())
ds_in = None
ds_out = None

I can print the GSD and it looks correct:

src = gdal.Open('dhm25_grid_raster.tif')
_, xres, _, _, _, yres  = src.GetGeoTransform()
print(xres, yres)
# 25.0 -25.0

Question 2: Is the yres supposed to be a negative value?

The negative value threw me so I tried using gdalwarp to convert the .asc straight to WGS84:

gdalwarp -s_srs EPSG:21781 -t_srs EPSG:4326 -dstnodata -32767.0 -r cubic -ot Int16 -of GTiff C:/somewhere/dhm25_grid_raster.asc C:/somewhere/warp25.tif

The resulting .tif looks to be correct. However when I load it and print the GSD, I get this:

src = gdal.Open('warp25.tif')
_, xres, _, _, _, yres  = src.GetGeoTransform()
print(xres, yres)
# 0.0003034045846788908 -0.0003034045846788908 !!!!

I tried setting the -tr 25.0 -25.0 flag but got the error:

ERROR 1: Attempt to create 0x0 dataset is illegal,sizes must be larger than zero.

I tried setting the -tr 25.0 25.0 flag but got the same error.

Question 3: How do I get gdalwarp to preserve the 25m GSD?

Solution

    1. I think LN02 is only the datum, related to the vertical dimension. It's needed (perhaps) to correctly interpret what the height values actually mean, how is 0.0 meters defined etc. https://en.wikipedia.org/wiki/Vertical_datum

      The geospatial projection seems to be EPSG:21781 indeed.

    2. The negative y-value is fine, it just means that the y-coordinate decreases when moving down ("~south") from the top. This is similar to EPSG:4326 for example, which starts with 90° at the North Pole, and decreases to -90°.

    3. You change the projection, which also changes the spatial unit from meters to degrees. So the value of ~0.0003 is the degree equivalent of 25 meters. If you don't set it, gdalwarp will automatically create this best approximation, but it's also possible to set it yourself.

      It's sometimes nice for example to have more rounded values for the bounding box and resolution. It can avoid errors when typing it manually, or avoid floating point precision issues etc. But that completely depends on the use case of course.

    If I assign EPSG:21781 in QGIS and compare it with Google Maps (which isn't necessarily perfectly located), at least the ballpark location looks reasonable to me.

    enter image description here