[GIS] Generate hillshading with GDAL, and respect nodata

Question

I'm still trying to nut out this problem of merging overlapping geotiffs.

From the original dataset (here, a 20m resolution digital elevation model of the state of Victoria), I have reprojected it and converted it to tiff. Rendered in TileMill, there is a clean, irregular boundary:

for f in dtm20m dtm10m_nw; 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 

The TileMill style is:

.dem {
  raster-colorizer-default-mode: linear;
  raster-colorizer-default-color: yellow;//hsl(60,50%,80%);//transparent;
  raster-colorizer-stops:
  stop(0,hsl(60,50%,80%))
  stop(392,hsl(110,80%,20%))
  stop(785,hsl(120,70%,20%))
  stop(1100,hsl(100,0%,50%))
  stop(1370,white);
  raster-scaling:near;
}

Next I generate hillshading:

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

Unfortunately now data seems to be generated even where there was none before.

The gdaldem documentation claims to respect nodata:

For all algorithms, except color-relief, a nodata value in the target
dataset will be emitted if at least one pixel set to the nodata value
is found in the 3×3 window centered around each source pixel. The
consequence is that there will be a 1-pixel border around each image
set with nodata value. From GDAL 1.8.0, if -compute_edges is
specified, gdaldem will compute values at image edges or if a nodata
value is found in the 3×3 window, by interpolating missing values.

I'm using GDAL 1.9.

I'm stuck at this point – I'm not sure how to even directly verify whether there is a 'nodata' value in the reprojected tif, or by what process TileMill knows how to leave transparent pixels. Any pointers would be greatly appreciated.

(For context, if needed, the reason all this matters is my current approach to the merging problem is to layer multiple hillshaded tifs on top of each other, which won't work if the "nodata" bits collide.)

Answer 1

I have managed to find a workaround, although it relies on there being an "extent" shapefile that came with the original DEM. I'd still like to know what I'm doing wrong in my original method, if anything.

gdalwarp -co "BIGTIFF=YES" -dstalpha -cutline dtm20m_ext_vg94.shp dtm20m-3785-hs.tif dtm20m-3785-hs-cut.tif

This produces the right kind of TIF:

...so I can composite two different resolutions of terrain:

(If it's not obvious, the seam runs roughly NW/SE through the middle. 20m resolution on the left, 90m on the right.)