Merge overlapping DEMs with irregular boundaries

Question

I have a 10m DEM covering part of a state, a 20m DEM covering the whole state, and a 90m DEM covering the whole country. How can I merge them together, getting the best quality available in each area?

The 10m and 20m are provided in ADF format, but I've converted them to .tif using ogr2ogr.

Using gdal_merge doesn't work, because the boundaries of each area are irregular. This is the result:

(The northeastern edge of the state is irregular, but there obviously shouldn't be the gap in the middle.)

UPDATE If it's useful, I also have a shapefile for each 10m DEM, marking its extent.

UPDATE 2

Ok, here's what the whole state looks like with the 20m DEM.

The process I'm going through so far is this:

for f in dtm20m dtm10m_nw dtm10m_e; do
export GDAL_CACHEMAX=1000
echo -n "Re-projecting: "
gdalwarp  -s_srs EPSG:3111 -t_srs EPSG:3785 -r bilinear ./vmelev_${f}/${f} ${f}-3785.tif 

echo -n "Generating hill shading: "
gdaldem hillshade -z 5 $f-3785.tif $f-3785-hs.tif
echo and overviews:
gdaladdo -r average $f-3785-hs.tif 2 4 8 16

done

Answer 1

I've done this using the workflow I explain below. I'll make the explanation for only two DEMs, with three you just have to repeat the process twice, the first time to merge the coarse and medium DEMs, and the second time to merge the resulting one with the finer DEM. For step 1 I use Global Mapper and for the following steps, Matlab, but the workflow would be the same using other software.

Let's say DEM "C" is the coarser and DEM "F" is the finer

1. Export both DEMs to geotiff using the same boundaries and the pixel size of the finer DEM. For that coarse DEMs will be interpolated. So both GeoTiff files will have the same size.
2. Create a DEM "D" corresponding to the difference C-F
3. Make void (Null, NaN, -9999 or whatever correspond to no data) all the pixels of D that are not in the edge of the patch of valid data. In this way D will only have a ring of valid values that correspond to the pixels on the edge of F.
4. Interpolate D to fill the void we just create (use Natural interpolation or some other weighted distance method, otherwise you might get sharp ridges that will mess the data).
5. Make F = D + F. In this way we have transform F so that it perfectly matches C on the edges but preserve the details of the relief of the original F.
6. Replace the void values of the adjusted F with the corresponding values of C.

That's it.

Might look complicated, but in Matlab is not many lines of code and the result is very good. If you are keen to use Matlab let me know and I can post some code.

This is a 3D visualization I've made of a DEM I patched this way and compared with SRTM v3. The annotations show some SRTM errors and the same areas in the fixed DEM.