External DEM handling in SNAP: consolidated findings on geocoding, Terrain Flattening, pixel alignment and vertical datum

Microwave Toolbox
1

Over the past few years, several STEP forum discussions have reported differences between SNAP’s built-in DEM handling and the use of external DEMs. The information is now spread across multiple threads and SNAP versions, so I would like to summarize the current state here and connect some of the findings.

The main topics appear to be:

  • differences between built-in and external DEM sampling;
  • a possible half-pixel offset in FileElevationModel;
  • handling of projected external DEMs, e.g. UTM;
  • Terrain Flattening-specific DEM resampling and oversampling;
  • antimeridian handling;
  • EGM96 vs. EGM2008 for Copernicus DEM;
  • which of the older issues have actually been fixed.

Some of these points are confirmed from the SNAP source code, some have been demonstrated experimentally, and others still need verification.

Relevant previous discussions

The following threads seem particularly relevant:

1. Original observation: external DEMs can give substantially worse Terrain Flattening results

The original discussion compared Terrain Flattening results using:

  1. SNAP’s built-in Copernicus DEM;
  2. the same Copernicus DEM supplied externally in geographic coordinates;
  3. an externally prepared DEM reprojected to UTM.

The built-in Copernicus DEM gave the best result. The external geographic DEM was somewhat worse, while the UTM DEM showed substantially more residual terrain structure.

At the time, the main suspected causes were DEM resampling, grid alignment and Terrain Flattening oversampling.

Some of these issues were subsequently addressed. In particular, S1TBX 9.0.3 included:

See also the S1TBX 9.0.3 release.

The Terrain Flattening operator had indeed calculated the DEM oversampling factor differently for built-in and external DEMs. The fix for this discrepancy was discussed in the forum and included under SNAP-2548.

However, the external DEM discrepancy was still reproducible after SNAP 9.0.3. This suggests that the oversampling issue was only one part of the problem.

2. External DEMs appear to be sampled half a pixel off

A much more recent investigation provides a likely explanation for at least part of the long-standing discrepancy:

External DEM is sampled half a pixel off

The same Copernicus DEM data were accessed once through SNAP’s built-in DEM mechanism and once through externalDEMFile. The external DEM was apparently sampled approximately half a DEM pixel to the southeast.

This can also be explained directly from the current source code.

FileElevationModel.getElevation() first converts the requested geographic coordinate to a pixel position:

final PixelPos pix = tileGeocoding.getPixelPos(geoPos, null);

and then calls:

resampling.computeCornerBasedIndex(
    pix.x,
    pix.y,
    RASTER_WIDTH,
    RASTER_HEIGHT,
    newIndex
);

Current source: FileElevationModel.getElevation(), lines 92–100

However, computeCornerBasedIndex() explicitly adds 0.5 to both coordinates:

default void computeCornerBasedIndex(
        double x, double y, int width, int height, Index index) {
    computeIndex(x + 0.5, y + 0.5, width, height, index);
}

Current source: Resampling.computeCornerBasedIndex(), lines 88–90

If GeoCoding.getPixelPos() already returns the appropriate pixel position for the geographic coordinate, this introduces an additional +0.5,+0.5 shift before interpolation.

A likely fix would therefore be to use:

resampling.computeIndex(...)

instead of:

resampling.computeCornerBasedIndex(...)

inside FileElevationModel.getElevation().

As of the current snap-engine master revision, the computeCornerBasedIndex() call is still present.

Why this matters more for a coarse or reprojected DEM

For a native Copernicus GLO-30 DEM, half a pixel corresponds to roughly 10–15 m horizontally, depending on latitude.

For a DEM reprojected to 100 m UTM resolution, the same half-pixel offset corresponds to approximately:

50 m east
50 m south

or roughly 70 m diagonally.

This could explain why a UTM 100 m DEM showed much stronger residual terrain effects than the original geographic DEM. The difference may therefore not be caused mainly by the 30 m → 100 m resolution reduction itself.

3. Projected external DEMs are not treated as projected terrain grids in Terrain Flattening

One might expect that supplying a DEM in UTM would allow SNAP to operate directly on a regular metric terrain grid.

This is not what Terrain Flattening currently does.

For an external DEM, the operator derives a single angular DEM resolution and an oversampling factor from FileElevationModel:

FileElevationModel filedem = (FileElevationModel) dem;

demResolution = filedem.getPixelWidthInDegrees();

final double minSpacing = Math.min(rangeSpacing, azimuthSpacing);

overSamplingFactor =
    Math.ceil(filedem.getPixelWidthInMeters() / minSpacing)
    * oversamplingMultiple;

Current source: TerrainFlatteningOp external DEM resolution calculation, lines 570–580

The operator then constructs a new regular geographic DEM grid:

final int rows =
    (int) Math.round((latMax - latMin) / demResolution);
final int cols =
    (int) Math.round((lonMax - lonMin) / demResolution);

final double[][] height = new double[rows][cols];

for (int i = 0; i < rows; ++i) {
    final double lat = latMax - i * demResolution;

    for (int j = 0; j < cols; ++j) {
        final double lon = lonMin + j * demResolution;

        geoPos.setLocation(lat, lon);
        double alt = dem.getElevation(geoPos);
        height[i][j] = alt;
    }
}

Current source: TerrainFlatteningOp construction of the geographic working DEM, lines 582–609

Terrain Flattening therefore effectively does:

external UTM DEM
    ↓
transform geographic query coordinates to UTM DEM coordinates
    ↓
resample DEM
    ↓
new regular geographic DEM
    ↓
oversample geographic DEM
    ↓
terrain-facet simulation

In other words, reprojecting a DEM to UTM before giving it to SNAP does not avoid reprojection inside Terrain Flattening. It adds another resampling step before SNAP reconstructs the DEM on a geographic grid.

This is quite different from workflows such as GAMMA, where one can first prepare the DEM in the target map projection and use this projected DEM grid directly for forward geocoding.

4. Terrain Flattening assumes one geographic DEM resolution for both axes

There is another potential problem for projected DEMs.

FileElevationModel calculates its effective pixel width in degrees from the geographic coordinates of the left and right raster edges:

public double getTileWidthInDegrees() {
    final GeoPos gpUL = getGeoPos(new PixelPos(0, 0));
    final GeoPos gpUR =
        getGeoPos(new PixelPos(RASTER_WIDTH - 1, 0));

    return Math.abs(gpUL.lon - gpUR.lon);
}

public double getPixelWidthInDegrees() {
    return getTileWidthInDegrees() / RASTER_WIDTH;
}

It separately implements pixel height, but Terrain Flattening uses only getPixelWidthInDegrees() for its single demResolution.

Current source: FileElevationModel DEM dimension/resolution methods, lines 162–187

Terrain Flattening then applies this same demResolution to both:

delta longitude = demResolution
delta latitude  = demResolution

For a projected DEM with square metric pixels, this does not preserve square spacing after conversion to geographic coordinates.

For example, a 100 m UTM pixel at 60° latitude corresponds to approximately:

100 m east-west
100 m north-south

but the angular longitude interval corresponding to 100 m is roughly twice the angular latitude interval corresponding to 100 m.

Using the longitude-derived angular spacing for both axes therefore produces an internal geographic base grid corresponding approximately to:

100 m east-west
200 m north-south

at 60° latitude.

Approximate values are:

Latitude Original UTM pixel Effective SNAP TF base grid
30° 100 × 100 m ~100 × 115 m
45° 100 × 100 m ~100 × 141 m
50° 100 × 100 m ~100 × 156 m
60° 100 × 100 m ~100 × 200 m
70° 100 × 100 m ~100 × 292 m

This is particularly important at high latitudes.

5. getPixelWidthInMeters() is also based on the angular longitude width

For external DEM oversampling, Terrain Flattening uses:

filedem.getPixelWidthInMeters()

The implementation is:

public double getPixelWidthInMeters() {
    return getPixelWidthInDegrees()
        * Constants.DTOR
        * Constants.semiMinorAxis;
}

Current source: FileElevationModel.getPixelWidthInMeters(), lines 178–187

This converts degrees of longitude to metres using an Earth-radius conversion without a latitude-dependent cos(latitude) factor.

For a projected 100 m pixel transformed to geographic coordinates, this means that at 60°N the reported getPixelWidthInMeters() can be approximately 200 m even though the actual projected raster pixel is 100 m wide.

That value subsequently influences the Terrain Flattening oversampling factor.

6. Oversampling cannot recover terrain detail lost during the initial geographic sampling

The Terrain Flattening oversampling is applied only after the intermediate geographic DEM has already been generated:

final double delta =
    demResolution / overSamplingFactor;

final int nLat =
    (int) (overSamplingFactor * rows);

final int nLon =
    (int) (overSamplingFactor * cols);

Current source: TerrainFlatteningOp oversampling after working-DEM construction, lines 610–616

The processing sequence is therefore approximately:

external DEM
    ↓
sample to regular geographic base grid
    ↓
height[][]
    ↓
oversample/interpolate height[][]
    ↓
terrain facets

If the projected DEM has already been sampled onto an effectively coarser or anisotropic geographic grid, increasing oversamplingMultiple creates a denser interpolation grid but cannot restore terrain information that was lost during the first sampling step.

This could explain why increasing oversamplingMultiple did not resolve the differences reported in the older forum thread.

7. Why this matters strongly for Terrain Flattening

The DEM in Terrain Flattening is not used merely to obtain an elevation for geometric correction. Four terrain vertices are converted to WGS84 coordinates and their projected triangle areas are used to estimate the locally illuminated area:

final PosVector t00 =
    lg.t00geo2xyzWGS84.getXYZ(
        lg.t00Lon, lg.t00Height);

...

// project points onto the plane perpendicular to slant range

...

gamma0Area[j] =
    areaTriangle1 + areaTriangle2;

Current source: TerrainFlatteningOp.computeIlluminatedArea(), starting at line 1301

This means DEM resampling and spatial displacement directly alter the local facet geometry and therefore the radiometric normalization.

A reduction from 30 m to 100 m DEM resolution can certainly affect Terrain Flattening more strongly than ordinary terrain correction.

However, the combination of:

  • the external DEM half-pixel shift;
  • the extra UTM → geographic resampling;
  • the anisotropic geographic working grid;
  • and the original 30 → 100 m resampling

provides a more plausible explanation for the very large degradation reported for the UTM DEM than the resolution reduction alone.

8. Range-Doppler Terrain Correction handles the DEM differently

The projected-DEM working-grid issue described above is specific to Terrain Flattening.

Range-Doppler Terrain Correction creates its requested output grid first. For each target tile it creates a TileGeoreferencing based on that target product and then requests a local DEM:

final TileGeoreferencing tileGeoRef =
    new TileGeoreferencing(
        targetProduct,
        x0 - 1,
        y0 - 1,
        w + 2,
        h + 2);

final boolean valid = DEMFactory.getLocalDEM(
    dem,
    demNoDataValue,
    demResamplingMethod,
    tileGeoRef,
    x0,
    y0,
    w,
    h,
    sourceProduct,
    nodataValueAtSea,
    localDEM);

Current source: RangeDopplerGeocodingOp target-grid DEM sampling, lines 945–954

DEMFactory.getLocalDEM() then converts each target-grid pixel to geographic coordinates and calls the DEM:

tileGeoRef.getGeoPos(x, y, geoPos);

...

alt = dem.getElevation(geoPos);

Current source: DEMFactory.getLocalDEM(), lines 405–416

So even if:

target CRS = UTM
DEM CRS    = the same UTM CRS

the processing is conceptually:

target UTM pixel
    ↓
latitude / longitude
    ↓
FileElevationModel
    ↓
UTM DEM coordinates

rather than directly addressing the corresponding DEM pixel in the projected grid.

The geographic round trip itself should be numerically harmless if implemented correctly, but the FileElevationModel half-pixel problem affects this lookup as well.

9. Difference to GAMMA-style geocoding

This architectural difference may be worth discussing.

In GAMMA, a common workflow is approximately:

source DEM
    ↓
reproject to target map projection
    ↓
projected DEM, e.g. UTM
    ↓
forward geocoding / lookup-table generation
    ↓
SAR data mapped into the same projected geometry

The DEM working grid and the final map geometry are therefore directly related.

In SNAP Terrain Flattening the workflow is instead approximately:

external projected DEM
    ↓
geographic DEM queries
    ↓
new regular lat/lon terrain grid
    ↓
SAR geometry

and Range-Doppler Terrain Correction later performs:

requested output map grid
    ↓
lat/lon
    ↓
DEM lookup
    ↓
radar geometry

Therefore selecting UTM as the final Terrain Correction CRS does not make Terrain Flattening operate in UTM.

The two approaches are fundamentally different.

10. Antimeridian handling

One motivation for using a projected DEM is processing scenes crossing the antimeridian.

A single regular geographic GeoTIFF with longitudes represented conventionally in [-180°, 180°] is awkward for such scenes, while a suitable projected DEM can represent the area continuously.

However, using a projected DEM does not appear to make Terrain Flattening itself antimeridian-safe.

TerrainFlatteningOp.computeImageGeoBoundary() retrieves the four geographic corner positions and simply calculates the ordinary minimum and maximum longitude:

final double[] lons = {
    geoPosFirstNear.getLon(),
    geoPosFirstFar.getLon(),
    geoPosLastNear.getLon(),
    geoPosLastFar.getLon()
};

double lonMin = 180.0;
double lonMax = -180.0;

for (double lon : lons) {
    if (lon < lonMin) {
        lonMin = lon;
    }

    if (lon > lonMax) {
        lonMax = lon;
    }
}

Current source: TerrainFlatteningOp.computeImageGeoBoundary(), lines 713–753

A scene containing corners around:

+179.5°
-179.5°

can therefore be interpreted as spanning approximately 359° instead of approximately 1°.

There is a second potential antimeridian issue in FileElevationModel itself:

return Math.abs(gpUL.lon - gpUR.lon);

Current source: FileElevationModel.getTileWidthInDegrees(), lines 162–166

For a projected DEM crossing the antimeridian, the transformed left and right edges could have longitudes close to +180° and -180°, causing the raster to appear almost 360° wide.

This seems to indicate that antimeridian handling needs to be implemented explicitly inside Terrain Flattening rather than worked around by supplying a projected external DEM.

11. Vertical datum: Copernicus DEM uses EGM2008, SNAP applies EGM96

Another independent issue is discussed here:

Copernicus DEM built-in path applies EGM96 while the data is EGM2008

Copernicus DEM heights are referenced to EGM2008.

The current Copernicus DEM implementation nevertheless explicitly initializes EarthGravitationalModel96:

public CopernicusElevationTile(
        ElevationModel demModel,
        Product product) throws IOException {

    super(demModel, product);
    egm = EarthGravitationalModel96.instance();
}

Current source: CopernicusElevationTile constructor, lines 20–23

The EGM correction is subsequently added to the DEM heights here:

line[i] += egm.getEGM(
    geoPos.lat,
    geoPos.lon,
    v);

Current source: CopernicusElevationTile.addGravitationalModel(), lines 50–57

The same FileElevationModel mechanism can apply EGM96 to an external DEM when externalDEMApplyEGM=true.

In Terrain Flattening the external DEM is initialized here:

dem = new FileElevationModel(
    externalDEMFile,
    demResamplingMethod,
    externalDEMNoDataValue);

((FileElevationModel) dem)
    .applyEarthGravitionalModel(
        externalDEMApplyEGM);

Current source: TerrainFlatteningOp external elevation model initialization, lines 1030–1043

Using EGM96 is appropriate for DEMs referenced to EGM96, such as SRTM, but not for the original Copernicus DEM, which uses EGM2008.

As far as I can see, this is still unresolved.

A possible current workaround is to transform Copernicus heights to ellipsoidal heights outside SNAP using EGM2008 and use:

externalDEMApplyEGM = false

12. External DEM EGM GUI/default-value problems

There was additionally a GUI/default-value problem where the visible state of the “Apply Earth Gravitational Model” checkbox could disagree with the actual externalDEMApplyEGM parameter.

Relevant tickets are:

SNAP-3958 is listed among the fixes in Microwave Toolbox 12.0.0. Further affected operators remained in SNAP 12.0.0 and were subsequently addressed under SNAP-4021, which is listed as fixed in Microwave Toolbox 12.0.1.

So this GUI/default-value problem appears to be fixed in current versions.

It should not, however, be confused with the separate EGM96-vs-EGM2008 issue.

13. Copernicus DEM tile-boundary artifacts

Another related but independent issue concerns invalid or modified pixels at the borders of automatically downloaded Copernicus DEM tiles.

This was reported as:

The issue was reported as fixed and is included in the SNAP 11 release notes.

However, similar tile-border artifacts have subsequently been reported again, including with SNAP 13 and SNAP 14:

Artefacts after Terrain Flattening of Sentinel-1 with COP DEM

A SNAP 13.0.2 Terrain Flattening regression was tracked separately as:

and was fixed in Microwave Toolbox 13.0.3.

The later tile-border artifacts reported in SNAP 13.0.4 and SNAP 14 therefore appear not to be identical to SNAP-4148.

14. Summary of the apparent issues

Issue Effect Current status
External DEM sampled ~½ pixel southeast Wrong elevations used by TC, TF and other DEM-dependent operators Apparently still present in FileElevationModel
External TF oversampling calculated differently from built-in DEM Different / inferior TF results Addressed by SNAP-2548 / S1TBX 9.0.3
Projected external DEM converted to a geographic working grid in TF Extra resampling; projected grid not preserved Current design
TF uses external DEM pixel width as both lat/lon resolution Anisotropic effective resolution for projected DEMs Apparently still present
getPixelWidthInMeters() derives metres from longitude degrees without latitude correction Incorrect physical resolution estimate away from equator Apparently still present
Simple min/max longitude handling in TF Potential ~360° extent for antimeridian scenes Apparently unresolved
getTileWidthInDegrees() uses direct longitude difference Potentially incorrect width for projected DEM crossing ±180° Apparently unresolved
Copernicus EGM2008 treated using EGM96 Vertical height error Apparently unresolved
externalDEMApplyEGM GUI/default inconsistency EGM could be applied unexpectedly SNAP-3958 / SNAP-4021, fixed
Copernicus auto-download tile-boundary artifacts Lines / invalid pixels in TC/TF SNAP-3689 fixed one occurrence; similar artifacts later reported again
TF partial-scene/NaN regression in SNAP 13.0.2 Large invalid regions SNAP-4148, fixed in 13.0.3

15. Questions for the SNAP developers

It would be very helpful if the SNAP developers could comment on the following points.

1. FileElevationModel pixel convention

Is this sequence intentional?

final PixelPos pix =
    tileGeocoding.getPixelPos(geoPos, null);

resampling.computeCornerBasedIndex(
    pix.x,
    pix.y,
    RASTER_WIDTH,
    RASTER_HEIGHT,
    newIndex);

Given that computeCornerBasedIndex() explicitly adds 0.5, should this instead use computeIndex()?

The recent external DEM comparison strongly suggests that the current code produces a half-pixel shift.

2. Are projected external DEMs officially supported?

The operators can apparently read a projected GeoTIFF through its geocoding, but Terrain Flattening derives its working resolution in geographic degrees and reconstructs the DEM on a regular geographic grid.

Is this expected behavior, or are external DEMs formally expected to be geographic/WGS84 only?

3. Why does Terrain Flattening use only getPixelWidthInDegrees()?

For a projected DEM, using the transformed X pixel spacing as both latitude and longitude spacing can result in a strongly anisotropic terrain grid.

Would it make sense to use separate X/Y resolutions?

4. Could Terrain Flattening operate directly on a projected DEM grid?

The terrain correction factor is ultimately calculated from the ECEF coordinates of four terrain-facet vertices.

It therefore does not seem fundamentally necessary for the source facet grid itself to be regular in latitude/longitude.

Supporting a projected terrain grid directly could:

  • avoid an unnecessary DEM reprojection;
  • make projected DEMs first-class inputs;
  • retain square metric terrain cells;
  • simplify antimeridian handling.

5. What is the intended solution for antimeridian-crossing scenes?

Using a projected external DEM seems attractive, but Terrain Flattening itself converts the processing extent back to lonMin/lonMax.

Should longitude unwrapping or circular-longitude handling be added to the internal TF grid generation?

6. Should Copernicus DEM use EGM2008 instead of EGM96?

The original Copernicus DEM is EGM2008-referenced.

Is there already a ticket or planned change to distinguish the vertical reference system by DEM type rather than applying EGM96?

16. Possible controlled test

One experiment might help separate the different effects.

Starting from exactly the same Copernicus DEM data:

A: original geographic DEM at native resolution
B: UTM DEM at approximately the same physical resolution, e.g. 30 m

run Terrain Flattening with both.

Then repeat after changing only:

computeCornerBasedIndex(...)

to:

computeIndex(...)

inside FileElevationModel.getElevation().

If the interpretation above is correct, I would expect:

Current SNAP:
    geographic DEM and UTM DEM differ

After fixing the half-pixel lookup:
    the difference becomes significantly smaller

After additionally fixing the projected-DEM TF grid handling:
    the two should become very similar, apart from the unavoidable
    interpolation introduced during the original DEM reprojection

This would help separate:

  • the half-pixel external DEM issue;
  • projected-vs-geographic DEM handling;
  • actual loss caused by DEM reprojection/resampling.

I hope this summary helps consolidate the different observations. Several older threads seem to describe different manifestations of related DEM handling mechanisms, while some tickets that were considered fixed addressed only one particular component.

Corrections, existing Jira tickets that I have missed, or clarification of the intended handling of projected external DEMs would be very welcome.

Disclaimer: drafted with the help of ChatGPT. Please forgive me. I love SNAP and just want to help.

2 Likes
2

Hi, I am the person who reported Artifacts at DEM tile borders. Originally it was an issue exclusively northern then 70 degrees north and without any distinction between flat or mountainous areas - they appeared elsewhere (comparing to SNAP 8, where it was not an issue at all). After first published fix, flat areas became more or less adequate, while mountainous still had those artifacts. For whatever reason with newer and newer releases it became worse - artifacts now appear even to the south of 70 degrees, in every place which do not have smooth topography of plain or lowland.

1 Like