Stitching DEM’s spanning CRS projections in QGIS

I didn’t think this was going to be so troublesome! but wanted to note down the steps to resolve.

The ‘Build Virtual Raster‘ process in QGIS will not bring together multiple rasters if they span across 2 CRS projections

Following my initial post of the issue on GIS stackexchange:
https://gis.stackexchange.com/questions/342830/cant-join-multiple-dems-using-qgis-build-virtual-raster-if-they-span-2-crss
A users tip re unifying the projections for loaded data sent me in the right direction.

Adding test selection of DEM’s from the following ALOS PALSAR dataset into QGIS:

As you can see we have only 1 DEM tile on 20S CRS, so we’ll see if we can simply re-project that one to match the others with 21S

First switch QGIS project CRS to 21S (bottom right globe icon)

Create shapefile for square extent in this case 250000 x 250000:

Extent: (229085.941613, 4146747.332439) – (479085.941613, 4396747.332439)

Run the Warp process on the single 20S DEM tile:

Or if you need to carry out the warp process on multiple dem files, just use the ‘Run as Batch Process‘ button bottom left of screen shot above. And select all the assets you want re-projected:

And fill out the appropriate fields:

Or you could use Python to carry out the same GDAL task:

import os
path=r"F:/GIS/Falklands/dem/"
for file in os.listdir(path):
    if file.endswith(".tif"):
        print path+"/"+file
        os.system('"gdalwarp.exe"'/
            +' -dstnodata -multi -s_srs EPSG:32720 '/
            +'-t_srs EPSG:32721 -r cubic '+path+"/"+file/
            +' F:/GIS/Falklands/dem/'+file+' -wo NUM_THREADS=8')

Layers over perfectly:

Now run Build virtual Raster selecting only the assets that are EPSG:32721

For some reason QGIS can’t render this new vrt very well:

But we will push forward and use the gdalwarp command to cut the new virtual raster above to the square extent created earlier:

Not forgetting to switch the CRS to 32721 which is what the virtual raster is set to:

gdalwarp -t_srs EPSG:32721 -wo SOURCE_EXTRA=1000 -tr 12.5 12.5 -srcnodata "-9999" -r cubic -ot Float32 -of GTiff -te 229085.941613 4146747.332439 479085.941613 4396747.332439 "F:/GIS/Falklands/virtual_final.vrt" F:/GIS/Falklands/dem_final.tif

All working nicely on UTM 21S CRS:

I’m not too happy having to process a DEM up to three times to get to a clean unified extract of the full dataset, so I’ll keep looking for an alternative, potentially doing it all via GDAL commands. But degradation of the data seems to be kept to a minimum, likley helped by using cubic convolution resampling.