Merge DEM tiles with different NoData values

How-to Last tested 24 October 2026

Keep each tile’s true NoData value in its band metadata, let gdalbuildvrt read those values, and use -vrtnodata to choose one output missing-value marker¹. Do not pass a common -srcnodata merely because the input sentinels differ: that override can make a missing value from one tile appear to be a valid elevation. The worked example below checks both real overlapping elevations and an explicitly constructed mixed-NoData copy.

Separate source meaning from output encoding

A NoData value is the marker a band uses for missing observations. Start by inspecting that value in each input with gdalinfo; do not infer it from a dark area in a map. Also check the CRS, geotransform, dimensions, band type and vertical-reference metadata. Agreement on a sentinel does not make incompatible grids suitable for mosaicking.

By default, gdalbuildvrt uses each source band’s intrinsic NoData. The generated XML records the source-specific exclusion inside its ComplexSource elements. -vrtnodata sets the VRT band’s output marker. These are different jobs: excluding an input cell versus representing an output location that has no contributing observation¹.

-srcnodata overrides the input-band value; a list supplies values for bands, not a different value for each file. If two single-band files use different correct sentinels, a single override cannot describe both. Retain their metadata, or deliberately normalise missing pixels and metadata together before using a common override. Changing a label without changing the relevant stored values is a different operation.

Use a pair with actual shared observations

The inputs are USGS_1M_13_x50y449_CO_EasternColorado_2018_A18.tif and its x51y449 neighbour, published on 26 May 2023. Their product IDs are 6480141fd34eac007b56ac82 and 64801405d34eac007b56ac31. Both official GeoTIFFs (western tile, eastern tile) and their metadata XML were downloaded; the acquisition record retains exact URLs, byte counts and SHA-256 hashes. Use the download guide for the acquisition workflow.

These rasters have 10012×10012 Float32 cells, one-metre spacing and NAD83/UTM zone 13N. Their XML describes NAVD88 elevations in metres². Both original bands use -999999. The western and eastern origins differ by nearly 10,000 metres, leaving a 12-column overlap. The script calculates the column offset from the stored transforms rather than treating the filenames as coordinates.

The source overlap contains 120,144 valid pairs. Of those, 872 differ; the maximum absolute difference is 0.0001220703125 metres. These are differences between delivered raster values, not an independent estimate of terrain accuracy. No tolerance was used to turn a small nonzero difference into equality.

Check order before introducing missing values

The first two VRTs use unmodified official files. east-last.vrt lists west then east; west-last.vrt reverses them. At all 872 differing pairs, the first result matches the eastern source and the second matches the western source. This gives observable evidence of source priority, instead of trying to identify the chosen source where both values are identical.

One recorded pair is 1456.52294921875 metres west and 1456.5228271484375 metres east. The two orderings return those respective source values. The mosaic is selecting an available value, not averaging the pair. Keep source order explicit when differences matter; the default behaviour documented for this GDAL release gives later sources priority while considering their NoData¹.

Construct the mixed-NoData case transparently

The script creates west-constructed.tif; it never edits either official input. It replaces the western tile’s 49,883,836 existing missing pixels with -32768 and assigns that band value as NoData. It also creates two artificial missing pixels: one where both original tiles have an observation, and one in the western-only footprint. Those holes are test inputs, not defects attributed to USGS.

The overlap hole is western pixel (10001, 1). Its original western value is 1530.009521484375; the corresponding eastern value is 1530.0096435546875. The western-only hole is pixel (9999, 5006), originally 1506.078369140625. Both coordinates are zero-based column/row offsets. The first exercises fallback to the other source; the second exercises an output with no remaining observation.

Set P10_WORK to a working directory containing the two verified official files and activate the recorded GDAL environment. This script creates a separate constructed raster, six VRTs, a JSON comparison report and the translated output. Keep those working rasters outside a published evidence package.

# SPDX-License-Identifier: MIT
set -euo pipefail
: "${P10_WORK:?Set P10_WORK to the directory containing the two verified USGS tiles}"
cd "$P10_WORK"
west=USGS_1M_13_x50y449_CO_EasternColorado_2018_A18.tif
east=USGS_1M_13_x51y449_CO_EasternColorado_2018_A18.tif
gdalinfo --version
python - <<'PY'
from osgeo import gdal
from pathlib import Path
import numpy as np, json, subprocess, xml.etree.ElementTree as ET
gdal.UseExceptions()
wname='USGS_1M_13_x50y449_CO_EasternColorado_2018_A18.tif'
ename=wname.replace('x50y449','x51y449')
w,e=[gdal.Open(n) for n in [wname,ename]]
wg,eg=w.GetGeoTransform(),e.GetGeoTransform()
assert w.GetProjection()==e.GetProjection() and wg[1:3]==eg[1:3] and wg[4:]==eg[4:]
off=round((eg[0]-wg[0])/wg[1]); width=w.RasterXSize-off
wa=w.ReadAsArray(off,0,width,w.RasterYSize); ea=e.ReadAsArray(0,0,width,e.RasterYSize)
wn,en=[d.GetRasterBand(1).GetNoDataValue() for d in [w,e]]
both=(wa!=wn)&(ea!=en); diff=both&(wa!=ea); assert both.any() and diff.sum()>1
locations=np.argwhere(diff); hy,hx=map(int,locations[0]); dy,dx=map(int,locations[-1])
cy,cx=w.RasterYSize//2,off-1; assert w.ReadAsArray(cx,cy,1,1)[0,0]!=wn
copy=gdal.Translate('west-constructed.tif',w,creationOptions=['TILED=YES','COMPRESS=DEFLATE'])
band=copy.GetRasterBand(1); rewritten=0
for y in range(0,w.RasterYSize,512):
    values=w.ReadAsArray(0,y,w.RasterXSize,min(512,w.RasterYSize-y)); missing=values==wn
    rewritten+=int(missing.sum()); values[missing]=-32768; band.WriteArray(values,0,y)
for x,y in [(hx+off,hy),(cx,cy)]:band.WriteArray(np.array([[-32768]],dtype='float32'),x,y)
band.SetNoDataValue(-32768); band.FlushCache(); copy=None
cases=[('east-last',[wname,ename],[]),('west-last',[ename,wname],[]),
       ('intrinsic',[ename,'west-constructed.tif'],[]),
       ('output',[ename,'west-constructed.tif'],['-vrtnodata','-9999']),
       ('wrong-source',[ename,'west-constructed.tif'],['-srcnodata',str(wn),'-vrtnodata','-9999']),
       ('ignore-source',[ename,'west-constructed.tif'],['-srcnodata','None','-vrtnodata','-9999'])]
points=[('valid-difference',dx+off,dy),('constructed-overlap-hole',hx+off,hy),
        ('constructed-west-only-hole',cx,cy)]
results={}
for name,inputs,options in cases:
    command=['gdalbuildvrt','-overwrite','-strict','-resolution','same',*options,name+'.vrt',*inputs]
    subprocess.run(command,check=True); ds=gdal.Open(name+'.vrt'); xml=ET.parse(name+'.vrt')
    strip=ds.ReadAsArray(off,0,width,w.RasterYSize)
    sources=[dict(file=s.findtext('SourceFilename'),nodata=s.findtext('NODATA'),mask=s.findtext('UseMaskBand'))
             for s in xml.findall('.//ComplexSource')+xml.findall('.//SimpleSource')]
    results[name]=dict(command=command,source_xml=sources,nodata=ds.GetRasterBand(1).GetNoDataValue(),
      samples={label:float(ds.ReadAsArray(x,y,1,1)[0,0]) for label,x,y in points},
      differing_pair_matches_west=int((strip[diff]==wa[diff]).sum()),
      differing_pair_matches_east=int((strip[diff]==ea[diff]).sum()))
assert results['east-last']['differing_pair_matches_east']==int(diff.sum())
assert results['west-last']['differing_pair_matches_west']==int(diff.sum())
assert results['intrinsic']['samples']['constructed-overlap-hole']==float(ea[hy,hx])
assert results['output']['samples']['constructed-west-only-hole']==-9999
assert results['wrong-source']['samples']['constructed-overlap-hole']==-32768
assert results['ignore-source']['samples']['constructed-overlap-hole']==-32768
report=dict(column_offset=off,overlap_columns=width,strip_cells=int(both.size),both_valid=int(both.sum()),
  valid_pair_differences=int(diff.sum()),max_abs_difference=float(np.max(np.abs(wa[both]-ea[both]))),
  mean_abs_difference=float(np.mean(np.abs(wa[both].astype('float64')-ea[both]))),
  original_nodata=[wn,en],rewritten_missing_cells=rewritten,constructed_holes=[[hx+off,hy],[cx,cy]],
  sample_inputs={'valid-difference':[float(wa[dy,dx]),float(ea[dy,dx])],
    'constructed-overlap-hole':[float(wa[hy,hx]),float(ea[hy,hx])],
    'constructed-west-only-hole':float(w.ReadAsArray(cx,cy,1,1)[0,0])},results=results)
Path('nodata-results.json').write_text(json.dumps(report,indent=2)+'\n'); print(json.dumps(report,indent=2))
PY
gdal_translate -of GTiff -co TILED=YES -co COMPRESS=DEFLATE output.vrt mosaic-e2.tif
gdalinfo -json mosaic-e2.tif > mosaic-e2-gdalinfo.json
gdalinfo mosaic-e2.tif > mosaic-e2-gdalinfo.txt
python - <<'PY'
from osgeo import gdal
import json
gdal.UseExceptions()
r=json.load(open('nodata-results.json')); d=gdal.Open('mosaic-e2.tif')
hx,hy=r['constructed_holes'][0]; cx,cy=r['constructed_holes'][1]
assert float(d.ReadAsArray(hx,hy,1,1)[0,0])==r['sample_inputs']['constructed-overlap-hole'][1]
assert float(d.ReadAsArray(cx,cy,1,1)[0,0])==-9999
print('Translated GeoTIFF: constructed overlap fallback and output NoData checks passed.')
PY

Read the XML and the resulting cells together

The following are actual values from the constructed-copy run. In every row, the input order is eastern original followed by western constructed copy. The output marker is separate from the two source markers recorded in the XML.

With intrinsic metadata, the constructed western hole is excluded and the eastern observation supplies the overlap. Changing only -vrtnodata retains that fallback and changes the unsupported western-only location to -9999. The incorrect source override instead treats -32768 as a contributing value; ignoring intrinsic source NoData does the same in these tested cells. Neither result is a sensible elevation at those locations.

The saved VRTs make the distinction inspectable: NODATA belongs to each source, while NoDataValue belongs to the output band³. The selected output VRT is translated to a GeoTIFF, then reopened. The final two assertions confirm that its overlap fallback and western-only missing value match the VRT results⁴.

GDAL 3.11.5; output: 20012 × 10012, Float32
CRS: NAD83 / UTM zone 13N; source XML: NAVD88, metres
Geotransform: (499993.9999999989, 1.0, 0.0, 4490006.000036587, 0.0, -1.0)
Output NoData: -9999
Translated overlap hole: 1530.0096435546875
Translated western-only hole: -9999

Download the recorded cell and VRT results for exact source-order counts and constructed-case values.

Keep the decision reproducible

For a production mosaic, retain the source hashes, original metadata, explicit input order, generated VRT and cell-check output. Check both a location where sources compete and one where a source is missing. A correct output NoData tag alone does not establish that the inputs were interpreted correctly. Use these separate checks to catch a leaked sentinel before it enters a downstream elevation calculation.

Not covered on this page

• Filling gaps or interpolating a replacement surface.
• Reconciling horizontal or vertical reference systems.
• Establishing survey accuracy or deciding which source is preferable.
• Alpha masks, multiband NoData lists or another GDAL release.