Skip to content

Commit

Permalink
Incremental changes
Browse files Browse the repository at this point in the history
  • Loading branch information
rhugonnet committed Dec 6, 2024
1 parent d5cade7 commit 71ed3e5
Show file tree
Hide file tree
Showing 2 changed files with 27 additions and 5 deletions.
1 change: 1 addition & 0 deletions doc/source/conf.py
Original file line number Diff line number Diff line change
Expand Up @@ -79,6 +79,7 @@
"geopandas": ("https://geopandas.org/en/stable", None),
"xarray": ("https://docs.xarray.dev/en/stable/", None),
"rioxarray": ("https://corteva.github.io/rioxarray/stable/", None),
"sliderule": ("https://slideruleearth.io/web/rtd/", None),
}

sphinx_gallery_conf = {
Expand Down
31 changes: 26 additions & 5 deletions doc/source/sliderule.md
Original file line number Diff line number Diff line change
Expand Up @@ -14,15 +14,15 @@ kernelspec:
---
(cheatsheet)=

# Elevation reference with SlideRule
# Pair with reference data

Most analysis of **xDEM relies on independent, high-precision elevation data to use as reference**, whether for
coregistration, bias-corrections or uncertainty analysis.

[SlideRule](https://slideruleearth.io/) provides the ideal way to retrieve such data by accessing big data archives
of high-precision elevations, such as ICESat-2 or GEDI, efficiently in the cloud.

Below, a short example to coregister and perform uncertainty analysis of a DEM with ICESat-2 ATL06.
Below, a short example to coregister our example DEM with ICESat-2 ATL06 reference data.

```{code-cell} ipython3
:tags: [remove-cell]
Expand All @@ -34,6 +34,10 @@ pyplot.rcParams['savefig.dpi'] = 600
pyplot.rcParams['font.size'] = 9 # Default 10 is a bit too big for coregistration plots
```

## Retrieving reference data

We load our example DEM data, and define the region of interest from its footprint.

```{code-cell} ipython3
from sliderule import sliderule, icesat2
import geoutils as gu
Expand All @@ -52,7 +56,11 @@ inlier_mask = ~glacier_outlines.create_mask(dem)
# Define region of interest as DEM footprint
bounds = list(dem.get_bounds_projected(4326))
region = sliderule.toregion(bounds)["poly"]
```

We can then initialize the SlideRule client, and fetch the ICESat-2 ATL06 reference data using {func}`icesat2.atl06sp`.

```{code-cell} ipython3
# Initialize SlideRule client
sliderule.init("slideruleearth.io")
Expand All @@ -71,15 +79,22 @@ gdf = gdf[np.isfinite(gdf["h_li"])] # Keep valid data
gdf = gdf[gdf["atl06_quality_summary"]==0] # Keep very high-confidence data
```

## Coregistration

Running the coregistration is as simple as passing the geodataframe to {func}`~xdem.coreg.Coreg.fit_and_apply`,
specifying the name of the column to use.

```{code-cell} ipython3
# Run a translation coregistration
nk = xdem.coreg.NuthKaab()
aligned_dem = nk.fit_and_apply(reference_elev=gdf, to_be_aligned_elev=dem, inlier_mask=inlier_mask, z_name="h_li")
aligned_dem = nk.fit_and_apply(reference_elev=gdf, to_be_aligned_elev=dem, inlier_mask=inlier_mask, z_name="h_li", random_state=42)
# Print the estimated translation parameters
print([k+f': {nk.meta["outputs"]["affine"][k]:.2f} meters' for k in ["shift_x", "shift_y", "shift_z"]])
```

We can plot the improvement due to the coregistration.

```{code-cell} ipython3
:tags: [hide-input]
:mystnb:
Expand All @@ -96,9 +111,13 @@ pc_mask = inlier_mask.astype("uint8").interp_points((gdf.geometry.x.values, gdf.
vect.ds = vect.ds[pc_mask == 1]
# Get point differences before and after
z_pc = dem.interp_points((vect.ds.geometry.x.values, vect.ds.geometry.y.values))
vect.ds["dh_bef"] = vect.ds["h_li"] - z_pc
z_pc_aligned = aligned_dem.interp_points((vect.ds.geometry.x.values, vect.ds.geometry.y.values))
vect.ds["dh_aft"] = vect.ds["h_li"] - z_pc_aligned
import warnings
# GeoPandas raising unexpected warning
with warnings.catch_warnings():
warnings.simplefilter("ignore")
vect["dh_bef"] = vect["h_li"] - z_pc
vect["dh_aft"] = vect["h_li"] - z_pc_aligned
# Plot before and after
import matplotlib.pyplot as plt
Expand All @@ -113,6 +132,8 @@ _ = ax[1].set_yticklabels([])
plt.tight_layout()
```

And the full information around the coregistration using {func}`~xdem.coreg.Coreg.info`.

```{code-cell} ipython3
# Show full coregistration summary
nk.info()
Expand Down

0 comments on commit 71ed3e5

Please sign in to comment.