Descargar un DEM gratis: cinco DEM de 30 m comparados con el láser GEDI

// Accuracy of five free 30 m DEMs against the GEDI lidar in Colombia and Venezuela
// Ventuari · https://ventuari.dev/tutoriales/descargar-dem-gratis/
//
// Compares Copernicus GLO-30, FABDEM, NASADEM, SRTM and ALOS AW3D30 with the
// ground elevation measured by GEDI in six landscapes: the Andes, the Llanos
// and the Amazon in Venezuela, and the Andes, the Caribbean and the Llanos in
// Colombia. Error = DEM − GEDI, both as heights above the WGS84 ellipsoid.
// Positive: the DEM sits above the ground (tree canopies, roofs).
//
// The interactive tables can time out (5-minute limit). The CSV export always
// finishes: click RUN on the "dem_vs_gedi_colombia_venezuela" task.

// ── 1. Geoids ──────────────────────────────────────────────────────────────
// GEDI measures heights above the ellipsoid; the DEMs, above a geoid.
// Copernicus and FABDEM use EGM2008; NASADEM, SRTM and ALOS use EGM96.
// Each DEM goes to the ellipsoid with ITS geoid: h = H + N.
var EGM96 = ee.Image('projects/earthengine-legacy/assets/users/ejvalero/EGM96-GEOIDE').resample('bilinear');
var EGM2008 = ee.Image('projects/earthengine-legacy/assets/users/ejvalero/EGM2008-GEOIDE').resample('bilinear');

// ── 2. Study areas (about 44 × 44 km each) ─────────────────────────────────
var ZONAS = [
  {pais: 'Venezuela', nombre: 'VE · Andes (Mérida)', rect: [-71.35, 8.40, -70.95, 8.80]},
  {pais: 'Venezuela', nombre: 'VE · Llanos (Apure)', rect: [-67.90, 7.55, -67.50, 7.95]},
  {pais: 'Venezuela', nombre: 'VE · Amazonía (Amazonas)', rect: [-66.40, 3.30, -66.00, 3.70]},
  {pais: 'Colombia', nombre: 'CO · Andes (Caldas)', rect: [-75.65, 4.90, -75.25, 5.30]},
  {pais: 'Colombia', nombre: 'CO · Caribe (Córdoba)', rect: [-76.05, 8.55, -75.65, 8.95]},
  {pais: 'Colombia', nombre: 'CO · Llanos (Meta)', rect: [-72.95, 3.90, -72.55, 4.30]}
];
var region = ee.FeatureCollection(ZONAS.map(function (z) {
  return ee.Feature(ee.Geometry.Rectangle(z.rect), {zona: z.nombre});
}));

// ── 3. Reference: ground measured by GEDI (L2A, 2019–2025) ─────────────────
var gediCol = ee.ImageCollection('LARSE/GEDI/GEDI02_A_002_MONTHLY')
  .filterBounds(region);
var gedi = gediCol.map(function (im) {
  var bueno = im.select('quality_flag').eq(1)
    .and(im.select('degrade_flag').eq(0))
    .and(im.select('elevation_bias_flag').eq(0))
    .and(im.select('sensitivity').gte(0.95))   // sees the ground under dense canopies
    .and(im.select('solar_elevation').lt(0));   // night shots only: less noise
  return im.select('elev_lowestmode').updateMask(bueno);
}).mosaic()
  .setDefaultProjection(gediCol.first().projection())   // sample on the GEDI grid (25 m)
  .rename('gedi');

// ── 4. The five DEMs, as heights above the ellipsoid ───────────────────────
// Bilinear interpolation on each image BEFORE mosaicking them.
function mosaico(id, banda) {
  return ee.ImageCollection(id).filterBounds(region)
    .map(function (im) { return im.select([banda]).resample('bilinear'); })
    .mosaic();
}
var DEMS = ['copernicus', 'fabdem', 'nasadem', 'srtm', 'alos'];
var dems = ee.Image.cat([
  mosaico('COPERNICUS/DEM/GLO30_2024_1', 'DEM').add(EGM2008),
  mosaico('projects/sat-io/open-datasets/FABDEM', 0).add(EGM2008),
  ee.Image('NASA/NASADEM_HGT/001').select('elevation').resample('bilinear').add(EGM96),
  ee.Image('USGS/SRTMGL1_003').select('elevation').resample('bilinear').add(EGM96),
  mosaico('JAXA/ALOS/AW3D30/V4_1', 'DSM').add(EGM96)
]).rename(DEMS);

// ── 5. Slope, land cover and change ────────────────────────────────────────
var cop = ee.ImageCollection('COPERNICUS/DEM/GLO30_2024_1').filterBounds(region).select('DEM');
var pendiente = ee.Terrain.slope(cop.mosaic().setDefaultProjection(cop.first().projection()))
  .rename('pendiente');
var cobertura = ee.ImageCollection('ESA/WorldCover/v200').first().select('Map').rename('clase');
var CLASES = [10, 20, 30, 40, 50, 60];   // no water, snow, wetlands or mangroves
var NOMBRES = ee.Dictionary({
  '10': 'Bosque', '20': 'Arbustal', '30': 'Herbazal',
  '40': 'Cultivo', '50': 'Construido', '60': 'Suelo desnudo'
});
// No forest loss 2001–2025: the terrain did not change between the DEM and GEDI.
// In Earth Engine, lossyear is masked where there was no loss: unmask(0) first,
// or eq(0) would drop exactly the pixels we want to keep.
var sinPerdida = ee.Image('UMD/hansen/global_forest_change_2025_v1_13')
  .select('lossyear').unmask(0).eq(0);

var pila = gedi.addBands(dems).addBands(pendiente).addBands(cobertura)
  .updateMask(sinPerdida)
  .updateMask(cobertura.remap(CLASES, [1, 1, 1, 1, 1, 1], 0));

// ── 6. Sample: up to 400 shots per land cover class in each area ───────────
var COLS = DEMS.map(function (d) { return 'e_' + d; });
var CUAD = DEMS.map(function (d) { return 'e2_' + d; });

var puntos = ee.FeatureCollection(ZONAS.map(function (z) {
  return pila.stratifiedSample({
    numPoints: 400,
    classBand: 'clase',
    region: ee.Geometry.Rectangle(z.rect),
    seed: 1,
    geometries: true,
    tileScale: 4
  }).map(function (p) { return p.set({zona: z.nombre, pais: z.pais}); });
})).flatten().map(function (p) {
  var g = ee.Number(p.get('gedi'));
  var errores = DEMS.map(function (d) { return ee.Number(p.get(d)).subtract(g); });
  var s = ee.Number(p.get('pendiente'));
  var tramo = ee.Algorithms.If(s.lt(5), '0–5°',
    ee.Algorithms.If(s.lt(15), '5–15°',
      ee.Algorithms.If(s.lt(30), '15–30°', 'más de 30°')));
  return p.set(ee.Dictionary.fromLists(COLS, errores))
    .set(ee.Dictionary.fromLists(CUAD, errores.map(function (e) { return e.pow(2); })))
    .set({
      cobertura: NOMBRES.get(ee.Number(p.get('clase')).int().format('%d')),
      tramo_pendiente: tramo,
      // If the median of the five errors exceeds 50 m, GEDI did not see the ground
      mediana_errores: ee.List(errores).reduce(ee.Reducer.median())
    });
}).filter(ee.Filter.lt('mediana_errores', 50))
  .filter(ee.Filter.gt('mediana_errores', -50));

// ── 7. Bias, median, NMAD and RMSE of each DEM ─────────────────────────────
var MEDIANA = ee.Reducer.median(null, null, 100000);   // exact up to 100 000 values

function resumen(sub, grupo) {
  var medias = ee.List(sub.reduceColumns(ee.Reducer.mean().repeat(10), COLS.concat(CUAD)).get('mean'));
  var medianas = ee.List(sub.reduceColumns(MEDIANA.repeat(5), COLS).get('median'));
  var DESV = DEMS.map(function (d) { return 'd_' + d; });
  var desvios = sub.map(function (p) {
    return p.set(ee.Dictionary.fromLists(DESV, COLS.map(function (c, i) {
      return ee.Number(p.get(c)).subtract(medianas.get(i)).abs();
    })));
  });
  var mad = ee.List(desvios.reduceColumns(MEDIANA.repeat(5), DESV).get('median'));
  var r = function (x) { return ee.Number(x).multiply(100).round().divide(100); };
  return ee.FeatureCollection(DEMS.map(function (d, i) {
    return ee.Feature(null, {
      fila: ee.String(grupo).cat(' · ').cat(d),
      n: sub.size(),
      sesgo_m: r(medias.get(i)),
      mediana_m: r(medianas.get(i)),
      nmad_m: r(ee.Number(mad.get(i)).multiply(1.4826)),
      rmse_m: r(ee.Number(medias.get(i + 5)).sqrt())
    });
  }));
}

function tabla(campo) {
  var grupos = puntos.aggregate_array(campo).distinct().sort();
  return ee.FeatureCollection(grupos.map(function (g) {
    return resumen(puntos.filter(ee.Filter.eq(campo, g)), g);
  })).flatten();
}

function mostrar(fc, titulo) {
  print(ui.Chart.feature.byFeature({
    features: fc, xProperty: 'fila',
    yProperties: ['n', 'sesgo_m', 'mediana_m', 'nmad_m', 'rmse_m']
  }).setChartType('Table').setOptions({title: titulo}));
}

mostrar(tabla('pais'), 'Por país');
mostrar(tabla('zona'), 'Por zona');
mostrar(tabla('cobertura'), 'Por cobertura (ESA WorldCover)');
mostrar(tabla('tramo_pendiente'), 'Por pendiente');

// ── 8. Points to check outside Earth Engine ────────────────────────────────
Export.table.toDrive({
  collection: puntos,
  description: 'dem_vs_gedi_colombia_venezuela',
  fileFormat: 'CSV',
  selectors: ['pais', 'zona', 'cobertura', 'tramo_pendiente', 'pendiente', 'gedi']
    .concat(DEMS).concat(COLS).concat(['.geo'])
});

Map.addLayer(region.style({color: 'ff5a1f', fillColor: '00000000', width: 2}), {}, 'Zonas');
Map.centerObject(region, 5);

"""
Accuracy of five free 30 m DEMs against the GEDI lidar in Colombia and Venezuela
Ventuari · https://ventuari.dev/tutoriales/descargar-dem-gratis/

Compares Copernicus GLO-30, FABDEM, NASADEM, SRTM and ALOS AW3D30 with the
ground elevation measured by GEDI in six landscapes: the Andes, the Llanos
and the Amazon in Venezuela, and the Andes, the Caribbean and the Llanos in
Colombia. Error = DEM − GEDI, both as heights above the WGS84 ellipsoid.
Positive: the DEM sits above the ground (tree canopies, roofs).

The summary tables can time out (5-minute limit for interactive requests).
The CSV export always finishes and holds every sampled shot.
"""

import ee
import geemap

ee.Authenticate()
ee.Initialize(project='tu-proyecto-gcp')


# 1. Geoids
# GEDI measures heights above the ellipsoid; the DEMs, above a geoid.
# Copernicus and FABDEM use EGM2008; NASADEM, SRTM and ALOS use EGM96.
# Each DEM goes to the ellipsoid with ITS geoid: h = H + N.
EGM96 = ee.Image('projects/earthengine-legacy/assets/users/ejvalero/EGM96-GEOIDE').resample('bilinear')
EGM2008 = ee.Image('projects/earthengine-legacy/assets/users/ejvalero/EGM2008-GEOIDE').resample('bilinear')


# 2. Study areas (about 44 × 44 km each)
ZONAS = [
    {'pais': 'Venezuela', 'nombre': 'VE · Andes (Mérida)', 'rect': [-71.35, 8.40, -70.95, 8.80]},
    {'pais': 'Venezuela', 'nombre': 'VE · Llanos (Apure)', 'rect': [-67.90, 7.55, -67.50, 7.95]},
    {'pais': 'Venezuela', 'nombre': 'VE · Amazonía (Amazonas)', 'rect': [-66.40, 3.30, -66.00, 3.70]},
    {'pais': 'Colombia', 'nombre': 'CO · Andes (Caldas)', 'rect': [-75.65, 4.90, -75.25, 5.30]},
    {'pais': 'Colombia', 'nombre': 'CO · Caribe (Córdoba)', 'rect': [-76.05, 8.55, -75.65, 8.95]},
    {'pais': 'Colombia', 'nombre': 'CO · Llanos (Meta)', 'rect': [-72.95, 3.90, -72.55, 4.30]},
]
region = ee.FeatureCollection([
    ee.Feature(ee.Geometry.Rectangle(z['rect']), {'zona': z['nombre']}) for z in ZONAS
])


# 3. Reference: ground measured by GEDI (L2A, 2019–2025)
gedi_col = ee.ImageCollection('LARSE/GEDI/GEDI02_A_002_MONTHLY').filterBounds(region)


def calidad_gedi(im):
    bueno = (
        im.select('quality_flag').eq(1)
        .And(im.select('degrade_flag').eq(0))
        .And(im.select('elevation_bias_flag').eq(0))
        .And(im.select('sensitivity').gte(0.95))   # sees the ground under dense canopies
        .And(im.select('solar_elevation').lt(0))   # night shots only: less noise
    )
    return im.select('elev_lowestmode').updateMask(bueno)


gedi = (
    gedi_col.map(calidad_gedi).mosaic()
    .setDefaultProjection(gedi_col.first().projection())   # sample on the GEDI grid (25 m)
    .rename('gedi')
)


# 4. The five DEMs, as heights above the ellipsoid
# Bilinear interpolation on each image BEFORE mosaicking them.
def mosaico(id_, banda):
    return (
        ee.ImageCollection(id_).filterBounds(region)
        .map(lambda im: im.select([banda]).resample('bilinear'))
        .mosaic()
    )


DEMS = ['copernicus', 'fabdem', 'nasadem', 'srtm', 'alos']
dems = ee.Image.cat([
    mosaico('COPERNICUS/DEM/GLO30_2024_1', 'DEM').add(EGM2008),
    mosaico('projects/sat-io/open-datasets/FABDEM', 0).add(EGM2008),
    ee.Image('NASA/NASADEM_HGT/001').select('elevation').resample('bilinear').add(EGM96),
    ee.Image('USGS/SRTMGL1_003').select('elevation').resample('bilinear').add(EGM96),
    mosaico('JAXA/ALOS/AW3D30/V4_1', 'DSM').add(EGM96),
]).rename(DEMS)


# 5. Slope, land cover and change
cop = ee.ImageCollection('COPERNICUS/DEM/GLO30_2024_1').filterBounds(region).select('DEM')
pendiente = ee.Terrain.slope(
    cop.mosaic().setDefaultProjection(cop.first().projection())
).rename('pendiente')
cobertura = ee.ImageCollection('ESA/WorldCover/v200').first().select('Map').rename('clase')
CLASES = [10, 20, 30, 40, 50, 60]   # no water, snow, wetlands or mangroves
NOMBRES = ee.Dictionary({
    '10': 'Bosque', '20': 'Arbustal', '30': 'Herbazal',
    '40': 'Cultivo', '50': 'Construido', '60': 'Suelo desnudo',
})
# No forest loss 2001–2025: the terrain did not change between the DEM and GEDI.
# In Earth Engine, lossyear is masked where there was no loss: unmask(0) first,
# or eq(0) would drop exactly the pixels we want to keep.
sin_perdida = (
    ee.Image('UMD/hansen/global_forest_change_2025_v1_13')
    .select('lossyear').unmask(0).eq(0)
)

pila = (
    gedi.addBands(dems).addBands(pendiente).addBands(cobertura)
    .updateMask(sin_perdida)
    .updateMask(cobertura.remap(CLASES, [1] * len(CLASES), 0))
)


# 6. Sample: up to 400 shots per land cover class in each area
COLS = ['e_' + d for d in DEMS]
CUAD = ['e2_' + d for d in DEMS]


def muestra_zona(z):
    return pila.stratifiedSample(
        numPoints=400,
        classBand='clase',
        region=ee.Geometry.Rectangle(z['rect']),
        seed=1,
        geometries=True,
        tileScale=4,
    ).map(lambda p: p.set({'zona': z['nombre'], 'pais': z['pais']}))


def errores_punto(p):
    g = ee.Number(p.get('gedi'))
    errores = [ee.Number(p.get(d)).subtract(g) for d in DEMS]
    s = ee.Number(p.get('pendiente'))
    tramo = ee.Algorithms.If(s.lt(5), '0–5°',
            ee.Algorithms.If(s.lt(15), '5–15°',
            ee.Algorithms.If(s.lt(30), '15–30°', 'más de 30°')))
    return (
        p.set(ee.Dictionary.fromLists(COLS, errores))
        .set(ee.Dictionary.fromLists(CUAD, [e.pow(2) for e in errores]))
        .set({
            'cobertura': NOMBRES.get(ee.Number(p.get('clase')).int().format('%d')),
            'tramo_pendiente': tramo,
            # If the median of the five errors exceeds 50 m, GEDI did not see the ground
            'mediana_errores': ee.List(errores).reduce(ee.Reducer.median()),
        })
    )


puntos = (
    ee.FeatureCollection([muestra_zona(z) for z in ZONAS]).flatten()
    .map(errores_punto)
    .filter(ee.Filter.lt('mediana_errores', 50))
    .filter(ee.Filter.gt('mediana_errores', -50))
)


# 7. Bias, median, NMAD and RMSE of each DEM
MEDIANA = ee.Reducer.median(None, None, 100000)   # exact up to 100 000 values
DESV = ['d_' + d for d in DEMS]


def resumen(sub, grupo):
    medias = ee.List(sub.reduceColumns(ee.Reducer.mean().repeat(10), COLS + CUAD).get('mean'))
    medianas = ee.List(sub.reduceColumns(MEDIANA.repeat(5), COLS).get('median'))
    desvios = sub.map(lambda p: p.set(ee.Dictionary.fromLists(DESV, [
        ee.Number(p.get(c)).subtract(medianas.get(i)).abs() for i, c in enumerate(COLS)
    ])))
    mad = ee.List(desvios.reduceColumns(MEDIANA.repeat(5), DESV).get('median'))

    def r(x):
        return ee.Number(x).multiply(100).round().divide(100)

    return ee.FeatureCollection([
        ee.Feature(None, {
            'fila': ee.String(grupo).cat(' · ').cat(d),
            'n': sub.size(),
            'sesgo_m': r(medias.get(i)),
            'mediana_m': r(medianas.get(i)),
            'nmad_m': r(ee.Number(mad.get(i)).multiply(1.4826)),
            'rmse_m': r(ee.Number(medias.get(i + 5)).sqrt()),
        }) for i, d in enumerate(DEMS)
    ])


def tabla(campo):
    grupos = puntos.aggregate_array(campo).distinct().sort()
    return ee.FeatureCollection(grupos.map(
        lambda g: resumen(puntos.filter(ee.Filter.eq(campo, g)), g)
    )).flatten()


def mostrar(fc, titulo):
    print('\n' + titulo)
    print(f"{'fila':36} {'n':>6} {'sesgo_m':>8} {'mediana_m':>10} {'nmad_m':>7} {'rmse_m':>7}")
    for f in fc.getInfo()['features']:
        p = f['properties']
        print(f"{p['fila']:36} {p['n']:>6} {p['sesgo_m']:>8} {p['mediana_m']:>10} {p['nmad_m']:>7} {p['rmse_m']:>7}")


# 8. Points to check outside Earth Engine (start it first: it always finishes)
tarea = ee.batch.Export.table.toDrive(
    collection=puntos,
    description='dem_vs_gedi_colombia_venezuela',
    fileFormat='CSV',
    selectors=['pais', 'zona', 'cobertura', 'tramo_pendiente', 'pendiente', 'gedi']
    + DEMS + COLS + ['.geo'],
)
tarea.start()
print('Export started:', tarea.status()['state'])

mostrar(tabla('pais'), 'Por país')
mostrar(tabla('zona'), 'Por zona')
mostrar(tabla('cobertura'), 'Por cobertura (ESA WorldCover)')
mostrar(tabla('tramo_pendiente'), 'Por pendiente')

m = geemap.Map()
m.add_basemap('SATELLITE')
m.addLayer(region.style(color='ff5a1f', fillColor='00000000', width=2), {}, 'Zonas')
m.centerObject(region, 5)
m