Generate GDAL and RichDEM ground truths for xDEM tests

.gitignore

@@ -0,0 +1,2 @@
+venv
+.idea

generate_ground_truth_gdal.py

@@ -0,0 +1,124 @@
+from osgeo import gdal
+import os
+import numpy as np
+import geoutils as gu
+
+
+def generate_gdal_truth(filepath: str, outdir: str) -> None:
+    """Run GDAL's DEMProcessing for all attibutes and save in .tif"""
+
+    # Define attributes and options
+    gdal_processing_attr_option = {
+        "slope_Horn": ("slope", "Horn"),
+        "aspect_Horn": ("aspect", "Horn"),
+        "hillshade_Horn": ("hillshade", "hillshade_Horn"),
+        "slope_Zevenberg": ("slope", "Zevenberg"),
+        "aspect_Zevenberg": ("aspect", "Zevenberg"),
+        "hillshade_Zevenberg": ("hillshade", "hillshade_Zevenberg"),
+        "tri_Riley": ("TRI", "Riley"),
+        "tri_Wilson": ("TRI", "Wilson"),
+        "tpi": ("TPI", None),
+        "roughness": ("Roughness", None),
+    }
+
+    # Execute GDAL for every attribute
+    for attr, (processing, option) in gdal_processing_attr_option.items():
+        print(f"Processing GDAL truth for: {attr}")
+
+        # Convert GDAL options
+        if option is None:
+            gdal_option = gdal.DEMProcessingOptions(options=None)
+        else:
+            gdal_option_conversion = {
+                "Riley": gdal.DEMProcessingOptions(alg="Riley"),
+                "Wilson": gdal.DEMProcessingOptions(alg="Wilson"),
+                "Zevenberg": gdal.DEMProcessingOptions(alg="ZevenbergenThorne"),
+                "Horn": gdal.DEMProcessingOptions(alg="Horn"),
+                "hillshade_Zevenberg": gdal.DEMProcessingOptions(azimuth=315, altitude=45, alg="ZevenbergenThorne"),
+                "hillshade_Horn": gdal.DEMProcessingOptions(azimuth=315, altitude=45, alg="Horn"),
+            }
+            gdal_option = gdal_option_conversion[option]
+
+        # Create file for saving GDAL's results
+        output = os.path.join(outdir, f"{attr}.tif")
+
+        # Execute GDAL
+        gdal.DEMProcessing(
+            destName=output,
+            srcDS=filepath,
+            processing=processing,
+            options=gdal_option,
+        )
+
+
+def gdal_reproject_horizontal_shift_samecrs(filepath_example: str, xoff: float, yoff: float, outdir: str):
+    """
+    Reproject horizontal shift in same CRS with GDAL for testing purposes.
+
+    :param filepath_example: Path to raster file.
+    :param xoff: X shift in georeferenced unit.
+    :param yoff: Y shift in georeferenced unit.
+    :param outdir: output directory.
+
+    :return: Reprojected shift array in the same CRS.
+    """
+
+    from osgeo import gdal, gdalconst
+
+    # Open source raster from file
+    src = gdal.Open(filepath_example, gdalconst.GA_ReadOnly)
+
+    # Create output raster in memory
+    driver = "MEM"
+    method = gdal.GRA_Bilinear
+    drv = gdal.GetDriverByName(driver)
+    dest = drv.Create("", src.RasterXSize, src.RasterYSize, 1, gdal.GDT_Float32)
+    proj = src.GetProjection()
+    ndv = src.GetRasterBand(1).GetNoDataValue()
+    dest.SetProjection(proj)
+
+    # Shift the horizontally shifted geotransform
+    gt = src.GetGeoTransform()
+    gtl = list(gt)
+    gtl[0] += xoff
+    gtl[3] += yoff
+    dest.SetGeoTransform(tuple(gtl))
+
+    # Copy the raster metadata of the source to dest
+    dest.SetMetadata(src.GetMetadata())
+    dest.GetRasterBand(1).SetNoDataValue(ndv)
+    dest.GetRasterBand(1).Fill(ndv)
+
+    # Reproject with resampling
+    gdal.ReprojectImage(src, dest, proj, proj, method)
+
+    # Extract reprojected array
+    array = dest.GetRasterBand(1).ReadAsArray().astype("float32")
+    array[array == ndv] = np.nan
+
+    filename = os.path.join(outdir, f"shifted_reprojected_xoff{xoff}_yoff{yoff}.tif")
+    output_driver = gdal.GetDriverByName("GTiff")
+    output_raster = output_driver.Create(
+        filename, dest.RasterXSize, dest.RasterYSize, 1, gdal.GDT_Float32
+    )
+    output_raster.SetGeoTransform(dest.GetGeoTransform())
+    output_raster.SetProjection(dest.GetProjection())
+    output_band = output_raster.GetRasterBand(1)
+    output_band.WriteArray(array)
+    output_band.SetNoDataValue(ndv)
+    output_raster.FlushCache()
+
+
+if __name__ == "__main__":
+    filepath = "data/Longyearbyen/DEM_2009_ref.tif"
+    outdir = "test_data/gdal"
+
+    # Create output directory if needed
+    os.makedirs(outdir, exist_ok=True)
+
+    # generate_gdal_truth(filepath, output_dir)
+
+    dem = gu.Raster(filepath)
+    list_off = [(dem.res[0], dem.res[1]), (10 * dem.res[0], 10 * dem.res[1]), (-1.2 * dem.res[0], -1.2 * dem.res[1])]
+    for x_off, y_off in list_off:
+        gdal_reproject_horizontal_shift_samecrs(filepath, x_off, y_off, outdir)

generate_ground_truth_richdem.py

@@ -0,0 +1,167 @@
+import os
+from typing import Union, List, Any
+from geoutils.raster import RasterType
+
+import numpy as np
+import geoutils as gu
+import richdem as rd
+
+# xDEM attribute name : RichDEM attribute name, degrees
+attributes_richdem = {
+            "slope_Horn": ("slope", True),
+            "aspect_Horn": ("aspect", True),
+            "hillshade_Horn": "hillshade",
+            "curvature": "curvature",
+            "profile_curvature": "profile_curvature",
+            "planform_curvature": ("planform_curvature", True),
+        }
+
+def raster_to_rda(rst: RasterType) -> rd.rdarray:
+    """
+    Convert geoutils.Raster to richDEM rdarray.
+    """
+    arr = rst.data.filled(rst.nodata).squeeze()
+    rda = rd.rdarray(arr, no_data=rst.nodata)
+    rda.geotransform = rst.transform.to_gdal()
+    return rda
+
+
+def get_terrainattr_richdem(rst: RasterType, attribute: str = "slope_radians") -> np.ndarray[Any]:
+    """
+    Derive terrain attribute for DEM opened with geoutils.Raster using RichDEM.
+    """
+    rda = raster_to_rda(rst)
+    terrattr = rd.TerrainAttribute(rda, attrib=attribute)
+    terrattr[terrattr == terrattr.no_data] = np.nan
+    return np.array(terrattr)
+
+
+def get_terrain_attribute_richdem(
+    dem: RasterType,
+    attribute: Union[str, List[str]],
+    degrees: bool = True,
+    hillshade_altitude: float = 45.0,
+    hillshade_azimuth: float = 315.0,
+    hillshade_z_factor: float = 1.0,
+) -> Union[RasterType, List[RasterType]]:
+    """
+    Derive one or multiple terrain attributes from a DEM using RichDEM.
+    """
+    if isinstance(attribute, str):
+        attribute = [attribute]
+
+    if not isinstance(dem, gu.Raster):
+        raise ValueError("DEM must be a geoutils.Raster object.")
+
+    terrain_attributes = {}
+
+    # Check which products should be made to optimize the processing
+    make_aspect = any(attr in attribute for attr in ["aspect", "hillshade"])
+    make_slope = any(
+        attr in attribute
+        for attr in [
+            "slope",
+            "hillshade",
+            "planform_curvature",
+            "aspect",
+            "profile_curvature",
+            "maximum_curvature",
+        ]
+    )
+    make_hillshade = "hillshade" in attribute
+    make_curvature = "curvature" in attribute
+    make_planform_curvature = "planform_curvature" in attribute or "maximum_curvature" in attribute
+    make_profile_curvature = "profile_curvature" in attribute or "maximum_curvature" in attribute
+
+    if make_slope:
+        terrain_attributes["slope"] = get_terrainattr_richdem(dem, "slope_radians")
+
+    if make_aspect:
+        # The aspect of RichDEM is returned in degrees, we convert to radians to match the others
+        terrain_attributes["aspect"] = np.deg2rad(get_terrainattr_richdem(dem, "aspect"))
+        # For flat slopes, RichDEM returns a 90° aspect by default, while GDAL return a 180° aspect
+        # We stay consistent with GDAL
+        slope_tmp = get_terrainattr_richdem(dem, "slope_radians")
+        terrain_attributes["aspect"][slope_tmp == 0] = np.pi
+
+    if make_hillshade:
+        # If a different z-factor was given, slopemap with exaggerated gradients.
+        if hillshade_z_factor != 1.0:
+            slopemap = np.arctan(np.tan(terrain_attributes["slope"]) * hillshade_z_factor)
+        else:
+            slopemap = terrain_attributes["slope"]
+
+        azimuth_rad = np.deg2rad(360 - hillshade_azimuth)
+        altitude_rad = np.deg2rad(hillshade_altitude)
+
+        # The operation below yielded the closest hillshade to GDAL (multiplying by 255 did not work)
+        # As 0 is generally no data for this uint8, we add 1 and then 0.5 for the rounding to occur between
+        # 1 and 255
+        terrain_attributes["hillshade"] = np.clip(
+            1.5
+            + 254
+            * (
+                np.sin(altitude_rad) * np.cos(slopemap)
+                + np.cos(altitude_rad) * np.sin(slopemap) * np.sin(azimuth_rad - terrain_attributes["aspect"])
+            ),
+            0,
+            255,
+        ).astype("float32")
+
+    if make_curvature:
+        terrain_attributes["curvature"] = get_terrainattr_richdem(dem, "curvature")
+
+    if make_planform_curvature:
+        terrain_attributes["planform_curvature"] = get_terrainattr_richdem(dem, "planform_curvature")
+
+    if make_profile_curvature:
+        terrain_attributes["profile_curvature"] = get_terrainattr_richdem(dem, "profile_curvature")
+
+    # Convert the unit if wanted.
+    if degrees:
+        for attr in ["slope", "aspect"]:
+            if attr not in terrain_attributes:
+                continue
+            terrain_attributes[attr] = np.rad2deg(terrain_attributes[attr])
+
+    output_attributes = [terrain_attributes[key].reshape(dem.shape) for key in attribute]
+
+    if isinstance(dem, gu.Raster):
+        output_attributes = [
+            gu.Raster.from_array(attr, transform=dem.transform, crs=dem.crs, nodata=-99999)
+            for attr in output_attributes
+        ]
+
+    return output_attributes if len(output_attributes) > 1 else output_attributes[0]
+
+
+def generate_ground_truth_richdem(dem: gu.Raster,
+                                  attribute_dict: dict[str, str | tuple[str, bool]],
+                                  save_dir: str) -> None:
+    """
+    Generate and save ground truth terrain attributes using RichDEM.
+
+    :param dem: DEM as a geoutils.Raster object.
+    :param attribute_dict: dictionary of richdem terrain attributes.
+    :param save_dir: Directory where the ground truth data will be saved.
+    """
+    # Create output directory if needed
+    os.makedirs(save_dir, exist_ok=True)
+
+    for attribute, value in attribute_dict.items():
+        # Generate the terrain attribute using RichDEM
+        if isinstance(value, tuple):
+            terrain_attribute = get_terrain_attribute_richdem(dem, attribute=value[0], degrees=value[1])
+        else:
+            terrain_attribute = get_terrain_attribute_richdem(dem, attribute=value)
+
+        # Save the generated attribute as a .npy file
+        save_path = os.path.join(save_dir, f"{attribute}.tif")
+        terrain_attribute.save(save_path)
+
+
+if __name__ == "__main__":
+    filepath = "data/Longyearbyen/DEM_2009_ref.tif"
+    output_dir = "test_data/richdem"
+    dem = gu.Raster(filepath, load_data=True)
+    generate_ground_truth_richdem(dem=dem, attribute_dict=attributes_richdem, save_dir=output_dir)