Conversion of INGV geotif DEM to ISCE DEM

[falkamelung]

The INGV dems have a mask layer inside the *tif. A conversion using rasterio (in the clean_ingv_dem.py worked. 44 is the geoid undulation for Italy (we need +44 in the DEM file).

rm ele*dem* ellipsoidal.tif cleaned.tif
clean_ingv_dem.py w42510_s10.tif cleaned.tif
gdal_calc.py -A cleaned.tif --outfile=ellipsoidal.tif --calc="A+44" --type=Int16
gdalwarp -t_srs EPSG:4326 -tr 0.000277777777778 0.000277777777778 -r bilinear -ot Int16 -of ISCE ellipsoidal.tif elevation.dem
gdalbuildvrt elevation.dem.vrt elevation.dem
fixImageXml.py -i elevation.dem -f

Here the code for the *tif cleaning. It is committed as clean_ingv_dem.py:

#!/usr/bin/env python3

import sys
import numpy as np
import rasterio

if len(sys.argv) != 3:
    print("Usage: python clean_dem.py <input_file> <output_file>")
    sys.exit(1)

input_file = sys.argv[1]
output_file = sys.argv[2]

# Open input DEM
with rasterio.open(input_file) as src:
    profile = src.profile.copy()
    data = src.read(1)

    # Detect NoData value
    nodata = src.nodata
    if nodata is None:
        print("⚠️  Warning: NoData value not set in original file. Assuming -9999.")
        nodata = -9999

    # Replace NoData with 0
    data_clean = np.where(data == nodata, 0, data)

    # Update profile
    profile.update({
        "dtype": "int16",
        "nodata": None
    })

    # Write to output
    with rasterio.open(output_file, "w", **profile) as dst:
        dst.write(data_clean.astype(np.int16), 1)