batch processing script

Author: trchudley

README.md

@@ -44,6 +44,7 @@ Easy filtering of ocean/mélange, allowing for assessment of calving fronts and
    </a>
 </p>
 
+An example batch download and coregistration script is also included in the `batch` directory as a jumping-off point for large-scale projects.
 
 # Cite
 
@@ -409,6 +410,7 @@ If you have requested multiple variables, the output of the `.pdt.terrain()` is
 
 | Version | Date | Notes |
 | ------- | ---- | ----- |
+| 0.4 | November 2023 | Made batch processing script available as an example. |
 | 0.3 | September 2023 | Aligned search function with the new PGC-provided GeoParquet files |
 | 0.2 | August 2023 | Update `load.mosaic()` function to include the new ArcticDEM mosaic v4.1 |
 | 0.1 | May 2023 | Initial release |

batch/batch_download_and_coregister.py

@@ -0,0 +1,128 @@
+"""
+This script automatically downloads and coregisters a complete history of ArcticDEM or
+REMA strips within a given AOI. The strips will be created within a named directory in 
+the current working directory.
+
+The script will download the 2 m ArcticDEM or REMA mosaic for coregistration, and then
+loop through the strip record, downloading and coregistering against stable ground
+identified within the BedMachine mask.
+
+Note: this will only work in regions of Greenland and Antarctica where bare rock
+is exposed. If you require coregistration in other regions, please drop me an email - 
+I would be happy to help.
+
+Tom Chudley | [email protected]
+Durham University
+November 2023
+"""
+
+import os
+import rioxarray as rxr
+import numpy as np
+import pdemtools as pdt
+
+# ----------------------------------------------------------------------------------- #
+# EDIT THIS SECTION TO SPECIFY DEPENDENT FILES, THE REGION, AND BOUNDS OF INTEREST
+# ----------------------------------------------------------------------------------- #
+
+# Filepath to ArcticDEM parquet file -- download from https://data.pgc.umn.edu/elev/dem/setsm/ArcticDEM/indexes/ArcticDEM_Strip_Index_latest_gpqt.zip
+index_fpath = ".../ArcticDEM_Strip_Index_s2s041_gpqt.parquet"
+
+# Filepath to BedMachine v5 netcdf file -- download from  https://data.pgc.umn.edu/elev/dem/setsm/ArcticDEM/indexes/ArcticDEM_Strip_Index_latest_gpqt.zip
+bm_fpath = ".../data/bedmachine_5/BedMachineGreenland-v5.nc"
+
+# The name of your study area - this will be used to name output files and directories
+region = "isunguatasermia"
+
+# "arcticdem" or "rema"
+dataset = "arcticdem"
+
+# the bounds of your study area, in EPSG:3413 for ArcticDEM or EPSG:3031 for REMA
+xmin, ymin, xmax, ymax = -247000, -2500000, -212000, -2487000
+
+# Parameters with which to filter ArcticDEM dataset. Note that you may wish to add
+# more or less - feel free to modify the search function at line 73.
+dates = "20010101/20231231"
+baseline_max_hours = 24
+min_aoi_frac = 0.1
+
+# ----------------------------------------------------------------------------------- #
+# END OF EDITABLE PARAMETER SECTION
+# ----------------------------------------------------------------------------------- #
+
+# define AOI
+bounds = (xmin, ymin, xmax, ymax)
+
+print(f"Downloading data for {region}:")
+
+# Create output directory
+outdir = f"{region}_data"
+if not os.path.exists(outdir):
+    os.mkdir(outdir)
+
+# Get reference DEM (ArcticDEM mosaic) if it doesn't already exist
+reference_dem_fpath = os.path.join(outdir, f"{region}_arcticdem_mosaic_2m.tif")
+
+if not os.path.exists(reference_dem_fpath):
+    print("\nDownloading reference DEM...")
+
+    reference_dem = pdt.load.mosaic(
+        dataset=dataset,  # must be `arcticdem` or `rema`
+        resolution=2,  # must be 2, 10, or 32
+        bounds=bounds,  # (xmin, ymin, xmax, ymax) or shapely geometry
+        version="v4.1",  # optional: desired version (defaults to most recent)
+    )
+    reference_dem.rio.to_raster(
+        reference_dem_fpath, compress="ZSTD", predictor=3, zlevel=1
+    )
+
+else:
+    print("\nLoading reference DEM...")
+    reference_dem = pdt.load.from_fpath(
+        os.path.join(outdir, f"{region}_arcticdem_mosaic_2m.tif"), bounds=bounds
+    )
+
+reference_dem = reference_dem.squeeze()
+bedrock_mask = pdt.data.bedrock_mask_from_bedmachine(bm_fpath, reference_dem)
+
+# Search for DEM strips
+print("\nSearching for DEM strips...")
+gdf = pdt.search(
+    index_fpath,
+    bounds,
+    dates=dates,
+    # months=[6, 7, 8, 9],
+    # years=[2019],
+    baseline_max_hours=baseline_max_hours,
+    # sensors=["WV03", "WV02", "WV01"],
+    # accuracy=2,
+    min_aoi_frac=min_aoi_frac,
+)
+gdf = gdf.sort_values("acqdate1")
+
+n_strips = len(gdf)
+
+print(f"{n_strips} strips found")
+
+i = 0
+
+print("\nDownloading DEM strips...")
+for _, row in gdf.iterrows():
+    date = row.acqdate1.date()
+    date_str = date.strftime("%Y%m%d")
+    dem_id = row.dem_id
+
+    out_fpath = os.path.join(outdir, f"{date_str}_{dem_id}_coreg.tif")
+
+    if not os.path.exists(out_fpath):
+        print(f"\nDownloading {i}/{n_strips} {os.path.basename(out_fpath)}...")
+        dem = pdt.load.from_search(row, bounds=bounds, bitmask=True)
+        dem.compute()  # rioxarray uses lazy evaluation, so we can force the download using the `.compute()` function.
+        dem = dem.rio.pad_box(*bounds, constant_values=np.nan)
+        dem = dem.pdt.coregister(reference_dem, bedrock_mask, max_horiz_offset=50)
+        dem.rio.to_raster(out_fpath, compress="ZSTD", predictor=3, zlevel=1)
+        del dem
+
+    i += 1
+
+print("Finished")

notebooks/get_icebergs.ipynb

@@ -229,7 +229,7 @@
    },
    "outputs": [],
    "source": [
-    "bedmachine_fpath = '/Users/pgwc34/Library/CloudStorage/OneDrive-DurhamUniversity/data/bedmachine_5/BedMachineGreenland-v5.nc'\n",
+    "bedmachine_fpath = '.../BedMachineGreenland-v5.nc'\n",
     "\n",
     "geoid = pdt.data.geoid_from_bedmachine(bedmachine_fpath, dem)"
    ]

notebooks/mosaic_and_terrain.ipynb

@@ -82,10 +82,9 @@ def coregister(
         reference: DataArray,
         stable_mask: Optional[DataArray] = None,
         return_stats: Optional[bool] = False,
-        max_horiz_offset: float = 15,
+        max_horiz_offset: float = 50,
         rmse_step_thresh: float = -0.001,
         max_iterations: int = 5,
-
     ) -> DataArray:
         """
         Coregisters the scene against a reference DEM based on the Nuth and Kääb (2011)
@@ -111,9 +110,10 @@ def coregister(
             residuals
         """
 
-        if geospatial_match(self._obj, reference) == False:
+        check_match = geospatial_match(self._obj, reference, return_info=True)
+        if not check_match == True:
             raise ValueError(
-                "Input DEM and reference DEM do not share geospatial information"
+                f"Input DEM and reference DEM do not share geospatial information: {check_match}. Consider padding (`.rio.pad_box`) or reprojecting (`.rio.repoject_match`) DEMs."
             )
 
         if stable_mask is None:

notebooks/strip_search_and_dem_difference.ipynb

@@ -25,7 +25,7 @@ def coregisterdems(
     dem2,
     mask,
     res,
-    max_horiz_offset=15,
+    max_horiz_offset=50,
     rmse_step_thresh=-0.001,
     max_iterations=5,
 ):
@@ -73,7 +73,10 @@ def coregisterdems(
 
         # Break loop if conditions reached
         if np.any(np.abs(pn[1:]) > max_horiz_offset):
-            print(f"Maximum horizontal offset ({max_horiz_offset}) exceeded")
+            print(
+                f"Maximum horizontal offset ({max_horiz_offset}) exceeded."
+                "Consider raising the threshold if offsets are large."
+            )
             return_meddz = True
             break
 
@@ -122,7 +125,7 @@ def coregisterdems(
         if rmse_step > rmse_step_thresh or np.isnan(d0):
             print(
                 "RMSE step in this iteration ({:.5f}) is above threshold ({}), "
-                "stopping and returning values of prior iteration".format(
+                "stopping and returning values of prior iteration. ".format(
                     rmse_step, rmse_step_thresh
                 )
             )
@@ -278,7 +281,7 @@ def interp2_gdal(X, Y, Z, Xi, Yi, interp_str, extrapolate=False, oob_val=np.nan)
 def dtype_np2gdal(dtype_np):
     # TODO: Write docstring.
 
-    if dtype_np == bool: #np.bool:
+    if dtype_np == bool:  # np.bool:
         promote_dtype = np.uint8
     elif dtype_np == np.int8:
         promote_dtype = np.int16

src/pdemtools/_accessor.py

@@ -25,7 +25,9 @@ def get_resolution(xds: DataArray) -> float:
         return abs(resolutions[0])
 
 
-def geospatial_match(rxd_1: DataArray, rxd_2: DataArray) -> bool:
+def geospatial_match(
+    rxd_1: DataArray, rxd_2: DataArray, return_info: bool = False
+) -> bool:
     """Check whether two (rio)xarray DataArrays or Datasets have the same geospatial
     information.
 
@@ -38,17 +40,24 @@ def geospatial_match(rxd_1: DataArray, rxd_2: DataArray) -> bool:
     :rtype: bool
     """
 
-    matches = [
-        rxd_1.rio.shape == rxd_2.rio.shape,
-        rxd_1.rio.resolution() == rxd_2.rio.resolution(),
-        rxd_1.rio.bounds() == rxd_2.rio.bounds(),
-        rxd_1.rio.crs == rxd_2.rio.crs,
-    ]
+    failed = []
 
-    if all(matches):
+    if not rxd_1.rio.shape == rxd_2.rio.shape:
+        failed.append("shape")
+    if not rxd_1.rio.resolution() == rxd_2.rio.resolution():
+        failed.append("resolution")
+    if not rxd_1.rio.bounds() == rxd_2.rio.bounds():
+        failed.append("bounds")
+    if not rxd_1.rio.crs == rxd_2.rio.crs:
+        failed.append("crs")
+
+    if len(failed) == 0:
         return True
     else:
-        return False
+        if return_info = True:
+            return failed
+        else:
+            return False
 
 
 def clip(

src/pdemtools/_coreg.py

@@ -18,7 +18,7 @@
 from ._utils import clip
 
 # arctic and rema valid versions for STAC retrival
-VERSIONS = {'arcticdem': ['v3.0', 'v4.1'], 'rema': ['v2.0']}
+VERSIONS = {"arcticdem": ["v3.0", "v4.1"], "rema": ["v2.0"]}
 
 # filenames of mosaic indexes in ./src/pdemtools/mosaic_index directory
 ARCTICDEM_V3_INDEX_FNAME = "ArcticDEM_Mosaic_Index_v3_gpkg.gpkg"
@@ -263,37 +263,47 @@ def mosaic(
     # sanity check that datset and versioning is correct versioning is valid for selected dataset
     dataset = dataset.lower()
     if dataset not in VERSIONS.keys():
-        raise ValueError(f'Dataset must be one of {VERSIONS.keys}. Currently `{dataset}`.')
+        raise ValueError(
+            f"Dataset must be one of {VERSIONS.keys}. Currently `{dataset}`."
+        )
     if version == None:  # pick the most recent dataset
         version = VERSIONS[dataset][-1]
-    else: 
+    else:
         if version not in VERSIONS[dataset]:
-            raise ValueError(f'Version of {dataset} must be one of {VERSIONS[dataset]}. Currently `{version}`')
+            raise ValueError(
+                f"Version of {dataset} must be one of {VERSIONS[dataset]}. Currently `{version}`"
+            )
 
     # get resolution as str
     if type(resolution) == int:
         resolution = f"{resolution}m"
 
     # check if valid resolution
     if resolution not in VALID_MOSAIC_RES:
-        raise ValueError(f"Resolution must be one of {VALID_MOSAIC_RES}. Currently `{resolution}`")
+        raise ValueError(
+            f"Resolution must be one of {VALID_MOSAIC_RES}. Currently `{resolution}`"
+        )
 
     # Sanitise shapely geometry to bounds tuple
     if type(bounds) == Polygon:
         bounds = bounds.bounds
 
     # get dataset version
-    if dataset == "arcticdem" and version == 'v3.0':
+    if dataset == "arcticdem" and version == "v3.0":
         layer = ARCTICDEM_V3_INDEX_2M_LAYER_NAME
-    elif dataset == 'arcticdem' and version =='v4.1':
+    elif dataset == "arcticdem" and version == "v4.1":
         layer = ARCTICDEM_V4_INDEX_2M_LAYER_NAME
-    elif dataset == "rema" and version == 'v2.0':
+    elif dataset == "rema" and version == "v2.0":
         layer = REMA_V2_INDEX_2M_LAYER_NAME
     else:
-        raise ValueError("Cannot retrive internal index filepath for specified dataset and version.")
+        raise ValueError(
+            "Cannot retrive internal index filepath for specified dataset and version."
+        )
 
     # Load tiles that intersect with AOI
-    tiles = gpd.read_file(_get_index_fpath(dataset, version=version), layer=layer, bbox=bounds)
+    tiles = gpd.read_file(
+        _get_index_fpath(dataset, version=version), layer=layer, bbox=bounds
+    )
 
     if len(tiles) < 1:
         raise ValueError(
@@ -339,14 +349,16 @@ def _get_index_fpath(
     """
 
     # get dataset version
-    if dataset == "arcticdem" and version == 'v3.0':
+    if dataset == "arcticdem" and version == "v3.0":
         fname = ARCTICDEM_V3_INDEX_FNAME
-    elif dataset == 'arcticdem' and version =='v4.1':
+    elif dataset == "arcticdem" and version == "v4.1":
         fname = ARCTICDEM_V4_INDEX_FNAME
-    elif dataset == "rema" and version == 'v2.0':
+    elif dataset == "rema" and version == "v2.0":
         fname = REMA_V2_INDEX_FNAME
     else:
-        raise ValueError("Cannot retrive internal index filepath for specified dataset and version.")
+        raise ValueError(
+            "Cannot retrive internal index filepath for specified dataset and version."
+        )
 
     return resources.files("pdemtools.mosaic_index").joinpath(fname)
 
@@ -362,7 +374,7 @@ def _aws_link(
     `PREFIX`, construct the filepath of the relevant ArcticDEM or REMA mosaic tile.
     """
     # Construct appropriate suffix, considering ArcticDEM v3.0's alternate naming scheme
-    if dataset == "arcticdem" and version=="v3.0":
+    if dataset == "arcticdem" and version == "v3.0":
         suffix = f"_{resolution}_{version}_reg_dem.tif"
     else:
         suffix = f"_{resolution}_{version}_dem.tif"