glaciohack / xdem Goto Github PK
View Code? Open in Web Editor NEWAnalysis of digital elevation models (DEMs)
Home Page: https://xdem.readthedocs.io/
License: MIT License
Analysis of digital elevation models (DEMs)
Home Page: https://xdem.readthedocs.io/
License: MIT License
I think this is good to do as soon as possible, before this package becomes more popular. Running pip install git+github.com/GlacioHack/DemUtils.git
only to import xdem
doesn't make much sense!
Right now it's still a mix. I'm up for homogenizing for in-function annotations everywhere (Erik's good at it, but I'm slowly learning ;D)
After discussion with @erikmannerfelt, we agreed to include the hillshade function of earthpy to xdem.
This seem to show different results than gdaldem hillshade by default, but can be solved by a zFactor of ~10-50.
Without scaling factor
https://files.slack.com/files-pri/T01D6EMSU0J-F01V0U7FWAH/figure_1.png
With scaling factor 50
https://files.slack.com/files-pri/T01D6EMSU0J-F01V47N24JX/figure_1.png
In setup.py
, richdem
is set as an extra_requires, but it is needed to import xdem
. In addition, scikit-gstat
is needed. geoutils
obviously a requirement, but since it is not yet on PyPi, it won't be downloaded by pip install xdem
. I don't know exactly how to fix that (apart from uploading geoutils on PyPi, of course)...
Right now, the "correct" way to install xdem without conda is quite complex:
pip install scikit-gstat
pip install richdem
pip install git+github.com/GlacioHack/GeoUtils.git
pip install git+github.com/GlacioHack/xdem.git # or pip install xdem :)
Skipping any of the above steps leads in an ImportError
when trying to import xdem
.
If we would add scikit-gstat
and richdem
as requirements, this would be simpler. Question is, however, if we want them as mandatory requirements? As far as I'm aware, richdem
only provides "nice to have" functions right now, unless @rhugonnet uses them more intricately somewhere? scikit-gstat
is mandatory for variograms, right? Could geoutils
be replaced somehow with git+git+github.com/GlacioHack/GeoUtils.git
so that pip
automatically fetches the latest master when trying to install xdem
?
I know that worrying about "ease of install" in such an early stage of development is kind of dumb. I feel, however, that it will be easier to keep up with these issues instead of delaying them to a "super-fix" later!
This includes:
Add option to sample_multirange_empirical_variogram():
Pick random point in raster + sample in rings of varying sizes to allow for a good distribution of samples for all distances and increase speed
Right now, the hypsometric approach (PR #36) only allows for equally spaced bins. An alternative to equal spacing is having bins with equal pixel counts.
Correct me if I'm wrong, but the advantage would be on glaciers with a combination of flat and steep surfaces. If a nonlinear elevation vs. elevation change relationship is assumed, flat surfaces would be under-represented in bins with equal spacing. Bins with equal pixel counts, on the other hand, would make sure that each part of the glacier is represented better.
Well, you can calculate a percentile bin, then use it in your percentile functions:
pbins = np.arange(min, max, bin_size)
bins = np.percentile(mean_dem, pbins)
Originally posted by @adehecq in #36 (comment)
Should we use the newer, older or mean elevation when calculating the mass balance?
We could write a small jupyter notebook on a test case that implements all solutions and look at the differences.
Currently, the implemented hypsometric approaches are local and (absolute) regional interpolation. Estimating a regional gradient using the normalized elevation range (0-1) for each glacier specifically has been shown by @rhugonnet to perform better. This should be included as an interpolation type somewhere.
In #71, we discussed errors associated to DEM gridding, since ICP (and possibly more future approaches) need to convert the DEM into a point cloud. ICP currently does:
Step 3 currently uses scipy.interpolate.griddata
which is extremely slow! A faster and more reliable approach would be fantastic to have. Something with the approximate signature:
def apply_matrix(dem: np.ndarray, transform: rio.transform.Affine, matrix: np.ndarray) -> np.ndarray:
# do matrix magic
return transformed_dem
One idea I had was to describe the 3D rigid transform as a 2D similarity transform and a linear vertical correction:
A 3D rigid transform from above (from 2D) would include horizontal translation, horizontal rotation, horizontal scaling (due to rotation-induced perspective differences), and elevation (bias + rotation) differences. As far as I understand, the first 3 parts could be described as a 2D similarity transform, and the final as a deramp.
The advantage of such an approach is that the DEM never has to be converted to a point cloud and regridded; it just needs to be moved and twisted a bit. I guess the most difficult part of this, however, would be to mathematically do the conversion...
This general-purpose function is integral for #79 to implement.
I am also all-ears to other alternatives.
Fixed for deramping in #60
Need to be fixed for amaury_coregister_dem.
The .copy()
method, which is inherited by the Raster
class, does not work as expected on the DEM class:
In [1]: import xdem
In [2]: dem = xdem.dem.DEM("examples/Longyearbyen/data/DEM_2009_ref.tif")
No metadata could be read from filename.
In [3]: dem_copy = dem.copy()
In [4]: assert isinstance(dem_copy, xdem.dem.DEM)
---------------------------------------------------------------------------
AssertionError Traceback (most recent call last)
<ipython-input-4-54169263df41> in <module>
----> 1 assert isinstance(dem_copy, xdem.dem.DEM)
AssertionError:
In [5]: type(dem_copy)
Out[5]: geoutils.georaster.Raster
In [6]: assert isinstance(dem, xdem.dem.DEM)
This could be argued to be a flaw in geoutils, but could also be fixed here.
Existing work here: https://github.com/iamdonovan/pyddem/blob/master/pyddem/fit_tools.py#L2093
To scale GP application matricially, kernel learning automated methods, need to add methods implemented in: https://gpytorch.ai/
Paper (1200+ citations in 10 years, I think this method is a winner): https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=5432191
Python implementation: https://github.com/siavashk/pycpd
There is an examples folder in the base folder, and when running the tests or examples, a new folder is created under xdem. I think it's best to keep only one...
Drafted in #60 by @adehecq, inspired by: https://github.com/jewettaij/superpose3d.
I am currently drafting a Coreg
subclass version of it.
Minimal example:
import geoutils as gu
import xdem
ref = xdem.dem.DEM(xdem.examples.FILEPATHS["longyearbyen_ref_dem"])
tba = xdem.dem.DEM(xdem.examples.FILEPATHS["longyearbyen_tba_dem"])
glacier_outlines = gu.geovector.Vector(xdem.examples.FILEPATHS["longyearbyen_glacier_outlines"])
tba_original = tba.copy()
xdem.coreg.coregister(ref, tba, mask=glacier_outlines)
assert tba_original == tba # Raises AssertionError
This will fail, because the tba
will be masked with the glacier_outlines
. This can obviously lead to some unexpected side-effects, but we should still discuss it before I do something about it:
If the DEMs in question are huge, I can't simply do tba.copy()
inside the coregister()
function due to memory limits. Maybe an optional keyword should exist, to choose between modifying inplace (default would be to copy) or not? What do you think, @rhugonnet and @adehecq?
Copying the code suggested by @rhugonnet, this (or something very similar) should be possible:
step1 = filters.NMAD()
step2 = coreg.BiasCorr()
step3. = filters.MedianKernel()
step4 = coreg.ICP()
step5 = biascorr.AlongTrack()
pipeline = step1 + step2 + step3 + step4 + step5
It either suggests a more general Pipeline
class, or just that the (not yet implemented) Filter
class just assumes it will be used in a CoregPipeline
when calling Filter.__add__()
Right now, all coreg
functions work on the entire DEM. This may or may not be necessary, and may consequently be a performance bottleneck on large DEMs. Perhaps coregistration can work almost equally well on a randomly subsampled set of pixels.
First of all, there are two idential functions:
xdem.coreg.calculate_nmad()
xdem.spatial_tools.nmad()
This is largely my fault, since I added my own in November when spatial_tools
existed but wasn't in my field of vision.
Second, both trigger odd warnings when passing masked_array
s to them:
In [1]: import xdem
In [2]: import numpy
In [3]: masked = np.ma.masked_array([1, 2, 3, 3], mask=[True, True, False, False])
In [4]: equivalent_array = np.array([3, 3])
In [5]: masked.mean() == equivalent_array.mean()
Out [5]: True
In [6]: masked_nmad = xdem.spatial_tools.nmad(masked)
/home/erik/.local/share/conda/xdem/lib/python3.8/site-packages/numpy/core/fromnumeric.py:753: UserWarning: Warning: 'partition' will ignore the 'mask' of the MaskedArray.
a.partition(kth, axis=axis, kind=kind, order=order)
In [7]: nmad = xdem.spatial_tools.nmad(equivalent_array)
In [8]: masked_nmad == nmad
Out [8]: False
In [9]: nmad
Out [9]: 0.0
In [10]: masked_nmad
Out [10]: 1.4826
For both masked
and equivalent_array
, I expected the NMAD to be 0. The output of the masked_nmad
suggests to me that it just returns the nfact
which is exactly 1.4826
.
I think the function has to check first if the array is a masked_array or not, then act accordingly.
Also, in the coreg.py
revision (#24), I should remove the duplicated calculate_nmad()
.
We need:
All available here: https://github.com/luc-girod/MMASTER-workflows/blob/master/pymmaster/mmaster_tools.py
See discussion in PR #85. For looping the local hypsometric interpolation, should we:
A local resource for ICP point-cloud alignment: https://github.com/ethz-asl/libpointmatcher
As suggested by @erikmannerfelt in #37
volume.linear_interpolation does not work as expected.
It is because of this line
Line 188 in 4e5ac64
Generalize the code in this to stack raster time series with low-memory writing on disk: https://github.com/iamdonovan/pyddem/blob/master/pyddem/stack_tools.py#L500
Get familiar with rioxarray to avoid duplicating functionalites developed in the past year for CRS parsing between NetCDF4 and rasters (https://github.com/iamdonovan/pyddem/blob/master/pyddem/stack_tools.py#L156)
Right now, the Longyearbyen test DEMs are void-free. Therefore, all the tests that pass are only guaranteed to pass with void-free examples. I've had large issues with the coreg functions and I think it is because of nans, so this should be added to further find out what's wrong with them.
I think this would be a good first step toward fixing #61.
PDAL is quite sensitive to the input data that it is given. The currently implemented approach can be improved dramatically by filtering the inputs from NaNs beforehand. This requires that the DEM is first converted to a point cloud (who could have thought PDAL worked better with those!?).
Right now, xdem.volume.hypsometric_binning()
allows for all kinds of binning alternatives.
This is currently not available in the DEMCollection.interpolate()
and dDEM.interpolate()
convenience methods (it will only use the default bins=50, kind="fixed"
arguments)
See discussion in PR #85. Some methods should be implemented to extrapolate hypsometric bins when they are missing at the bottom/top part of the glacier. @rhugonnet suggested the following:
0/ Existing: fill with the mean regional signal at a given elevation (Fanny's paper, and in Bob's paper, I don't like this much and I quote Bob's results here: the global hypsometric methods are furthest from the true values, which is perhaps not surprising in a region with a high variability of glacier elevation changes)
1/ Easy: loop through all local hypsometric pd.DataFrame, normalize those by their hypsometry and apply the area-weighted mean dhdt at the same normalized hypsometry to gap fill: at least a bit more representative of the regional signal independently of the elevation location at the scale of each glacier, but can still be off due to high inter-regional variability
2/ Best: loop through all local hypsometric pd.DataFrame, normalize by both hypsometry and elevation change and derive a mean normalized signal for the region; then, for glacier with local hypsometric (maybe with a criteria covering at least X% of the area), de-normalize by an estimated range of elevation change rate. If the rate at high elevation is not available (accumulation area full of voids), you can use the signal at half of the normalized elevation range, or maybe the best would be to use the mean elevation change rate at highest elevation captured in the region (essentially, you scale the regional hypsometric signal only by the dhdt in the ablation area of each glacier, keeping the accumulation dhdt similar than where measured for others in the region).
For 2/, the best would be a Bayesian approach, but you don't have time for this right now.
@adehecq I think that the detailed steps should be:
a) Normalize all local hypsometric estimation with more than X% coverage in accumulation area by both hypsometry and elevation change rate.
b) Derive area-weighted, mean hypsometric regional normalized signal from everything in a).
c) Derive area-weighted, mean regional elevation change rate at highest elevation from everything in a).
d) For glaciers with more than X% coverage in ablation area, de-normalize the regional signal of b) with the maximum dhdt of the glacier and the dhdt at highest elevation of c).
e) For glaciers with less than X% coverage in ablation area (glacier with pretty much no coverage), use with the mean elevation change rate in the region? (or for a more beautiful visual rendering, fill with the mean hypsometric signal without normalization)
To better fit the names of the script files, these may fit better in biascorr.py
.
I strongly suggest that they still inherit from the Coreg
class, however, to be fully compatible with each other. I suggest the same for e.g. biascorr.AlongTrack()
.
Raster mosaicking handling in GDAL just allows for "first" or "last" operations. Mean and median operations would make for much better mosaicking, and it would be nice to have it integrated with GeoUtils Rasters.
Minimal example:
In [1]: import xdem
In [2]: dem = xdem.DEM(xdem.examples.FILEPATHS["longyearbyen_ref_dem"], silent=True)
No metadata could be read from filename.
Should the silent
flag not silence the "No metadata could be read from filename."
text?
Running 'python DemUtils/coreg.py' fails with error:
File "DemUtils/coreg.py", line 196, in icp_coregistration
assert reference_dem.shape == dem_to_be_aligned.shape
AssertionError
This is because of several issues. Be aware that rasterio ds shapes will be (nbands, height, width), so even for single band raster, shape will be of length 3. Hence the following issues:
Also, at https://github.com/GlacioHack/DemUtils/blob/95cf00be26499b09830e300f273d25ceef4a3b57/DemUtils/coreg.py#L629 (and after), I would suggest to:
Should xdem
allow DEMs with the shape (1, height, width)
, (height, width)
or either of the two? Right now, there's no real consensus on which one works (or at least the first might not always work, but will in some cases).
(height, width)
, all functions should have:assert len(dem.shape) == 2, f"Given DEM has shape '{dem.shape}'. Expected 2D array"
(height, width)
:demc = dem.squeeze() # Create a squeezed copy and work on the copy
assert len(dem.shape) == 2, f"Given DEM has shape '{dem.shape}'. Expected 2D array or 3D array with shape[0] == 1"
(1, height, width)
:if len(dem.shape) == 3 and dem.shape[0] == 1:
demc = dem.copy()
elif len(dem.shape) == 2:
demc = dem.reshape((1,) + dem.shape)
else:
raise AssertionError(f"Given DEM has shape '{dem.shape}'. Expected 2D array or 3D array with shape[0] == 1")
I vote for either 1 or 2. What do you guys think?
It could be implemented in a def _check_dem(dem) -> None:
to avoid repeating code everywhere.
There is no consensus on the optimal coregistration pipeline, so adding modularity would help explore this.
It could be done in scikit-learn's fashion for regression models:
from sklearn.linear_model import LinearRegressor, RANSACRegressor
import numpy as np
X = np.linspace(0, 1, num=50).reshape(-1, 1)
y = np.linspace(0, 2, num=50)
for regressor in LinearRegressor, RANSACRegressor:
model = regressor()
model.fit(X, y)
print(model.score(X, y))
The point of the above example is to show how interchangeable the different regression models are.
The function sklearn.pipeline.make_pipeline()
allows for custom regression pipelines to be created. If for example the samples should be normalized before, this could easily be done:
from sklearn.linear_model import LinearRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import make_pipeline
model = make_pipeline([StandardScaler(), LinearRegressor(fit_intercept=True)])
model.fit(X,y)
# etc...
A similar approach could be done for coreg.py
:
from xdem import coreg
reference_dem = ....
dem_to_be_aligned = ....
other_dem = ....
model = coreg.AmauryCoreg()
model.fit(reference_dem, dem_to_be_aligned)
print(model.nmad)
model.transform(other_dem)
... or something similar.
If a custom pipeline is sought:
model = coreg.make_pipeline([coreg.BiasCorr(), coreg.OutlierRemoval(), coreg.ICPCoreg()])
model.fit(reference_dem, dem_to_be_aligned)
Basically, all coreg approaches should be subclasses of a BaseCoreg
class or similar, and all should have the methods .fit()
and .transform()
(or similar terminology). The make_pipeline()
could be adapted from the equivalent in scikit-learn.
If parts or all steps of a pipeline can be described as matrices, the current apply loop is inefficient and might induce resampling-related errors:
Line 1531 in 82b87eb
A better approach would be to merge transformation matrices (when possible) and therefore do fewer resampling iterations, potentially lowering the error and the processing time.
A conceptual example is:
Currently, the pipeline would modify the DEM four times. By combining the matrices of A, B and C, only two modifications are needed:
This could be done by first implementing a property of Coreg
, for example called is_rigid
which is True if it can be represented as a rigid transform (4x4 matrix). Then, the following pseudocode could work:
dem_mod = dem.copy()
matrices: list[np.ndarray] = []
for coreg in self.pipeline:
if coreg.is_rigid:
# Append the matrix to the list of matrices to be applied at the same time.
matrices.append(coreg.to_matrix())
else: # This marks when matrices can no longer be added, and need to be applied before the non-rigid apply is done
# Use matrix-magic to combine the matrices
combined_matrix = combine_matrix(matrices)
# Use more matrix-magic to transform the DEM
apply_matrix(dem_mod, combined_matrix)
# Reset the matrices list
matrices.clear()
# Apply the non-rigid transform
coreg.apply(dem_mod, transform=transform)
if len(matrices) > 0:
# do the same as above to apply all the matrices
where the combine_matrix()
and apply_matrix()
functions need to be created.
I give up on my attempt to do this for now... I can package xdem and use it on a fresh conda env, but scikit-gstat
and geoutils
are not on conda/conda-forge, so they can't be included. I don't know what to do at this point. I'll dump my progress here and then move on to something else for now..
meta.yaml
:
{% set setup_data = load_setup_py_data(setup_file='setup.py', from_recipe_dir=True) %}
package:
name: xdem
version: "{{ setup_data.get('version') }}"
build:
noarch: python
number: 3
requirements:
run:
- python
- numpy
- scipy
- rasterio
- geopandas
- pyproj
- tqdm
about:
summary: Tools for DEM comparisons and registration
I've read that conda-forge packages can be published automatically from the PyPi version. That would be great!
Right now, my (probably inefficient) way of building and testing is:
conda build .
anaconda upload "$(conda build . --output)" # Upload the result to the logged in channel (eriksh for me)
# Steps for setting up an environment
conda create -n clean_env python
conda activate clean_env
conda install -c eriksh xdem
pip install scikit-gstat git+https://github.com/GlacioHack/GeoUtils.git
# Now the command below should not result in an ImportError
python -c "import xdem"
For now, deramping and nuth_kaab algorithm do not return a transformation matrix. Returning a trasnformation matrix would make it easier to:
So far, the best solution I found is this library: https://github.com/jewettaij/superpose3d
There is a lot of unused code there, so it could be worth just rewriting some of the code.
This is more of @rhugonnet's area of expertise; I am just reiterating what we had on the hypsometric approach TODO list. Maybe someone can elaborate here?
@adehecq and I have discussed at multiple points in time that this would be a nice addition. For example, if one would want to use the results of the coregistration in another tool (ASP / Metashape etc.), 4x4 transformation matrices are the way to go.
For bias corrections, it would be as simple as:
[
[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, bias],
[0, 0, 0, 1]
]
for translations:
[
[1, 0, 0, xcorr],
[0, 1, 0, ycorr],
[0, 0, 1, bias],
[0, 0, 0, 1]
]
We also discussed that 1st degree deramps can be described as a bias, rotation and scale in the horizontal direction. 2nd degree and higher deramps are however impossible to describe as a 4x4 matrix, and could return a NotImplementedError
or similar.
PDAL and opencv already provide matrices from their ICP approaches, so this would be simple to return as a np.ndarray.
This could be part of the coreg.py
revision initiative (#24).
I could submit a PR, but my aggressive formatting would mess up your code, @rhugonnet ;)
If the user has a Raster
and wants to convert this to a DEM
, there is currently no easy way to do it. One hacky way is:
dem = xdem.dem.DEM.from_array(raster.data, raster.transform, raster.crs, raster.nodata)
A quite simple conversion could be done in a .from_raster()
static method. I recently learned this trick, and it's great because no copying is done in memory whatsoever!:
class DEM(...):
# all the old stuff
def from_raster(raster: gu.georaster.Raster, vref_name=None, vref_grid=None, silent=False):
dem = DEM.__new__(DEM)
dem.__dict__ = raster.__dict__
# user input
dem.vref = vref_name
dem.vref_grid = vref_grid
dem.ccrs = None
# trying to get vref from product name (priority to user input)
# dem.__parse_vref_from_fn(silent=silent)
return dem
This kind of works for me, but I just now realised that it won't work perfectly because it's supposed to be initialised through a SatelliteImage
, so it doesn't have all the right attributes.
Anyway, this would be a really nice addition!
There should be a quick way to figure out how well the coregistration went. NMAD is one alternative, but it could hide potential biases. I like RMSE's, but they often look dauntingly big!
The .error()
(or similar) method could alternatively take an argument to choose (metric="nmad"
).
Opinions?
Romain's presentation just made me think that we should make sure that by design all our coreg function allow passing weights. This way, we can later implement weighted coregistration.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
๐ Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. ๐๐๐
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google โค๏ธ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.