pyproj seems to sometimes not work properly from the conda-forge channel, and may need to be installed with pip instead.
conda remove -f pyproj
pip install pyproj
This code that you wrote is theoretically perfect! It's more or less what would be done in an element-wise assignment under the hood normally. An issue with this, however, is that python loops in practice are very slow. Maybe this works now, but if you have 8 million points, this would take days.
Lines 160 to 187 in 806d6e2
I expected this to be possible in many fewer lines, but I may be missing something. The DEM.value_at_coords()
method can take arrays of coordinates, which we can use in our favour to speed things up:
is2["dem_elevation"] = "index", dem.value_at_coords(is2["easting"], is2["northing"])
Note that I haven't tried this code, so it might be more difficult.
As far as I remember, the only xarray coordinate in the dataset is one called "index"
, right? There's a shorthand way of saying "assign these values with the following (already existing) coordinates", which is what I've done above. Basically, you assign a two-element tuple where the first is the name of the coordinate, and the second is the value array.
It might be that DEM.value_at_coords()
needs loaded values. In that case, you can change the above to one of these:
dem.value_at_coords(is2["easting"].values, is2["northing"].values)
# OR
dem.value_at_coords(is2["easting"].load(), is2["northing"].load())
When you have DEM elevations as a variable in the xarray dataset, you can easily make dh by simply subtracting the two:
is2["dh"] = is2["dem_elevation"] - is2["h_te_best_fit"]
An advantage to retaining the DEM elevation as a variable is that it should optimally be used for the hypsometric binning (#1) instead of IS2 elevations.
The hypsometric binning function is almost correct! xdem.volume.hypsometric_binning
takes two values of equal shape and bins them in a specified way. This means that the .shape
of the first argument must be equal to the second. Note that I mentioned in my email that this function won't work at all, but I take that back after looking at the code.
When you open the raster (L199), it will have a shape of, say, (3000, 2000)
, while your is2 data may have a shape of something like (10000, C)
(where C is the amount of columns). This will fail directly, because the hypsometric binning function has no code that transforms between different types of data (such as rasters to points).
Lines 195 to 222 in 806d6e2
What you need to do is to sample the DEM elevation and store it in the point data. That way, the shape will be equal. A quick (but perhaps not long-term great) way is to instead use the IS2 elevation:
def hypsometric_binning(is2_path):
data = xr.open_dataset(is2_path)
binned = xdem.volume.hypsometric_binning(ddem=data["dh"], ref_dem=data["h_te_best_fit"])
return binned
In the future, or now directly, we may want to use quantile binning (just add kind="quantile"
as an argument) as this will make comparisons between glaciers much easier. Also maybe ease down on the bin count (e.g. bins=10
).
Right now, we have a pretty good view of the dependencies that are required for the code to run. To make it easy to install on new computers, it's nice with an "environment file". This file is an instruction on how to recreate this environment using conda/mamba. It's relatively easy to set up:
name: SvalbardSurges
channels:
- conda-forge
dependencies:
- python >= 3.10
- geopandas
- xarray
- netcdf4
- matplotlib
- rasterio
- pip:
- variete
Note that the pyproj
hack that we made with pip is not present here. It could be added as a pip dependency, but the problem is (or at least should be) unique to your computer for some reason. Therefore it's best to not add this hack, but perhaps add that to the (not yet existing) README.md
:
## Installation
pyproj seems to sometimes not work properly from the conda-forge channel, and may need to be installed with pip instead.
conda remove -f pyproj
pip install pyproj
Identical to #9 apart from the type of data. See #9 for the rationale.
The same (long term) issue is present here with data having to be downloaded manually, but could be solved by just automatically downloading from the https://npolar.no website instead.
Lines 137 to 140 in 2fa3e53
Here's the link to the dataset repo:
https://data.npolar.no/dataset/f6afca5c-6c95-4345-9e52-cfe2f24c7078
And here's the direct link to the glacier outline dataset:
https://api.npolar.no/dataset/f6afca5c-6c95-4345-9e52-cfe2f24c7078/_file/3df9512e5a73841b1a23c38cf4e815e3
geopandas
can read shapefiles within zipfiles, so we don't have to worry about extracting them!
It's great that you've added the spread (std)! This will help us figure out a qualitative or quantitative signal to noise ratio on the plots.
SvalbardSurges/svalbardsurges/analysis.py
Lines 110 to 124 in 4430e91
Right now you're calculating a global spread (on all bins). Come to think of it, this might actually be really useful as a surging glacier should have a massive std while a nonsurging should have a small one. I think that an even more useful tool, however, would be a per-bin std:
stds = xdem.volume.hypsometric_binning(
ddem=data_subset["dh_corr"],
ref_dem=data_subset["dem_elevation"],
kind="custom",
bins=bins,
aggregation_function=np.nanstd
)
Then you'd have two hypsometric dataframes. We could experiment on different ways to combine them. I think it should just be as simple as first renaming the median df to something more useful:
hypso.rename(columns={"value": "median"}, inplace=True)
and then adding the standard deviation as a column:
hypso["std"] = stds["value"]
The glacier outlines from Geyman et al., 2022 (for some reason labelled GAO) seem to be saved as if it was a NetCDF instead of a zipfile. Also its name is is2
which of course sounds like IceSat-2!
SvalbardSurges/download_urls.py
Line 14 in aa4f21d
It's great that you're caching the result of the DEM differencing! There's however no cache retrieval logic so the process is rerun every time!
Lines 190 to 191 in 448be05
This is an easy fix by simply adding something almost identical to what's in subset_is2
:
Lines 67 to 71 in 448be05
Basically, make the cache_path
at the start of the function, then add the cache check just after.
A good dH filtering step is now in the hypsometric interpolation function. This is the first time it's needed (right now at least), but I think this would fit better inside the elevation change function instead, as this filtering should always apply, not just when hypsometric binning is performed.
Lines 285 to 286 in d543748
Speaking of filtering, the whole reason for that line was that nodata values are sampled. I wonder if this can be fixed by changing the sampling strategy a bit. The Raster.value_at_coords
method has a masked
keyword that should mask out any nodata.
https://geoutils.readthedocs.io/en/stable/gen_modules/geoutils.Raster.value_at_coords.html#geoutils.Raster.value_at_coords
A masked array cannot (to my knowledge) be set to an xarray dataset, so we need the .filled(np.nan)
method to replace masked values with nans:
dem.value_at_coords(is2.easting, is2.northing, masked=True).filled(np.nan)
This might solve the nodata problem completely if we're lucky!
I checked the hypsometric plotting function and it saves to a figures/
directory which is great! This will however fail if the figures/
directory doesn't already exist.
SvalbardSurges/svalbardsurges/plotting.py
Lines 61 to 63 in 4430e91
There are two ways of fixing this:
Path("figures/").mkdir(exist_ok=True)
if not os.path.isdir("figures/"):
os.mkdir("figures/")
The hypsometric binning seems to work for all data! Hooray!
The next step would be to bin it separately for each year (we can discuss other time periods like 6/3 months, but one year is commonplace).
Lines 293 to 295 in d543748
This is quite easy using the xarray.DataArray.dt
accessor that allows us to loop through every year easily:
for year, data_subset in data.groupby(data["date"].dt.year):
# run the hypsometric binning on `data_subset` instead of `data`
# Save the data by appending to a list or a dataframe
Then we have to merge the results somehow into one nice dataframe or dataset. Making an xarray dataset might be easiest, as it's so easy to work with multidimensional data. We could have a dataset with three dimensions: (glacier, time, elevation band). What do you think of that, @eliskasieglova?
Another consideration is that for nice comparison between years, we might have to pre-set the bins to use. Otherwise, the estimated bins may be slightly different depending on the year, which means we'll have a really hard time comparing different years. We could pre-define the bins to be, say, percentiles:
bins = np.nanpercentile(data["dem_elevation"], np.linspace(0, 100, 11))
This will generate bins for the 0-10th percentile, 10-20th, etc.
We can use these bins by just supplying them bins=bins
to the hypsometry function.
You bring up a good point with this comment!:
Line 105 in 806d6e2
There might be cases of duplicate names. This isn't great, and as you say, requires some kind of handling.
I see two solutions, where the choice between them almost needs a real-world example:
If you want to filter by the largest, it's pretty fast in (geo-)pandas:
gao_glacier["area"] = gao_glacier.geometry.area
return gao_glacier.sort_values("area").iloc[-1]
This first creates an area column (maybe it already exists? But best to make a new one in case we change outlines). Then it sorts by area, and extracts the last (largest) row. Note that this will fail with a very confusing error if the gao.query(...)
call was empty. Therefore, after the query, I suggest to add a check that there is at least one outline in the filtered data:
assert gao_glacier.shape[0] > 0, f"Query NAME=='{glacier_name}' returned zero features."
If you want to merge all touching outlines, you can run these two:
gao_glacier = gao_glacier.dissolve().explode()
This will first take the unary union of all outlines (dissolve()
) and then separate all outlines that don't touch into separate rows (explode()
). Note that most metadata are lost using this approach. But we don't need that right now, so I guess it's fine!
It's great that all functions are documented! There's a slightly better way to do this than to just add #
-comments after the definition like now:
Lines 64 to 65 in 806d6e2
Docstrings are special types of strings that are automatically parsed in the function creation when placed correctly. They also really help with larger code bases, as your editor will show them automatically when you start typing the function.
There are many different styling standards for docstrings. I personally prefer numpydoc as it is very easy to read and can also be parsed automatically by documentation tools. In the above case, the docstring could look something like:
def subset_is2(is2_data, bounds, label):
"""
Subset the IceSat-2 data using the specified bounds.
The easting and northing variables need to be loaded in memory, so this is a computationally expensive task.
The function is cached using the label argument to speed up later calls.
Parameters
----------
- is2_data
The IceSat-2 data to subset
- bounds
The bounding box to use (requires the keys "left", "right", "bottom", "top")
- label
A label to assign when caching the subset result
Returns
-------
A subset IS2 dataset within the given bounds.
"""
# The actual code goes here
This might seem awfully long for something that right now is obvious, but trust me, it's worth it. In a few months, you might have so many different functions that these docstrings will really help us to keep track of what's doing what!
Note the somewhat rigid structure of the docstring convention:
- arg
, followed by a new line, a tab, and then a description of the argumentThere are also more optional titles like "Raises" (for potential errors), "Examples", "See also" etc. We don't need them right now, but they might come in handy later.
There's nothing wrong with how the hypsometric binning function takes a path instead of a dataset object for the hypsometric binning. However, it is inconsistent with how the other functions work, which might make for problems further down the line.
Lines 264 to 283 in d543748
I suggest simply changing the argument of the function to take a dataset instead of a path. That way, in the main function, you don't have to hard-code the cache path either and can just give it the is2_dh
variable!
For each glacier, we only want the points that intersect with the particular glacier, of course. I recommended one approach which ended up being very computationally demanding. The current approach is:
This technically works, but requires that a clipped DEM is saved (for caching) each time the analysis is run. If we do this on all Svalbard's glaciers, we'd end up with about 1500 DEMs! There are some other issues:
A much simpler approach is to just clip the IS2 points:
The caching could be one of two approaches:
The latter might be a better option because xarray indexing on coordinates works really fast. So for each glacier, it'd just be a table (or a simple text file) of indices that should be used. Thus, there would be no data duplication and the cache would be much smaller.
We are soon (or already) at the point where we can start working on ways to properly detect glacier surges in your data! The easiest would just be to threshold the elevation change rate, say >2 m/a is a surge. This would work for all big surges, but perhaps not the small ones.
The best approach might be to go in a supervised fashion, i.e. by collecting a surge/no surge inventory that we can use to refine the best classification from. Kääb et al. (2023) made an inventory from 2017-2022 that includes Svalbard for large surges, which overlaps nicely with your timeframe. Note that it's only the large ones, and I unfortunately don't know of any other inventory (apart from a personal one which only covers Heer Land) with the smaller ones. But perhaps we can use Kääb et al., (2023) for now! The supplementary is just a table, but I have an equivalent shapefile which will be easier to work with.
The cache should preferably not be in the repo, as it's by definition temporary state that should be reproducible. Also, while there are currently only very small VRTs in there, it could blow up like crazy in the future, so we should be wary of committing the cache. It's very easy to configure git to ignore the directory completely.
It's usually made through a .gitignore
file:
cache/
It's a bit more convoluted when the cache is already in the history. First, you have to remove it from the history by either removing it or moving it. In a CLI, this would be:
>>> mv cache cache.backup # Move the cache somewhere else. It needs to look like it doesn't exist anymore.
>>> git add cache/ # Contradictory to what it looks like, this REMOVES the cache (basically it adds the change, which is a removal)
>>> git add .gitignore # This needs to be created and filled in
>>> git commit -m "Removed cache directory from history and added to .gitignore"
>>> mv cache.backup cache # Restore the cache
I think there are convoluted git commands to do all of this, but the above commands are the most straightforward.
For now, having just one main.py
works fine. In the long-term, however, this will become quite cumbersome and hard to navigate. When the time comes, I recommend splitting the code up into smaller modules that will increase readability.
For example:
svalbardsurges/
analysis.py
For hypsometric binning, change detection, etc.plotting.py
For making figuresutilities.py
For utility functions like generic downloading, cache handling, etc.inputs
A submodule for all the different inputs
dems.py
is2.py
outlines.py
Then the main.py
function simply imports the svalbardsurges
package and calls the necessary functions from within.
See some projects I've worked on for inspiration:
main.py
to a package! The main.py
in comparison looks awful now!We need good ways to visualize dH over time, and there are many different ways of doing that. Since your data are hypsometrically binned, a time series is technically two-dimensional (one is time, one is elevation). A few publications visualize this as a heatmap, i.e. an image where one axis is time and another is distance/elevation.
Here's a not-too-great example from an ongoing project of mine. I binned elevation change in percentile bins (0-10%, 10-20%, etc.) and plotted dH/dt and d²H/dt² (acceleration) against time:
The glacier surged in 2016 but the data are too coarse to even see it! What we do see, however, is that the lower part of the glacier gets more and more negative over time, leading up to the surge, and the upper part stops having net positive accumulation (it's not blue any more). The acceleration part is a bit interesting because it shows a relatively constant acceleration for that entire period, suggesting that the events leading up to the surge ultimately were quite predictable for a long time.
I may have some more exciting examples, but what do you think of this way of visualizing it?
Sevestre et al., (2018) used this visualization for velocities, which quite nicely show some intricacies of the surges they study.
In the long term (a few weeks), we should strive for making the code as repeatable as possible. What I mean by that is:
An example of non-repeatability is the hardcoding of paths that we've downloaded manually.
I want to underscore that we've of course done this for brevity and making sure we could get a fast initial result, so this is just a natural evolution from that!
Lines 172 to 173 in 2fa3e53
What would instead be better to do is:
Downloading data with python is a bit of a mess. To my knowledge, there's no one-liner for good downloading. In one of my projects, I wrote this script:
https://github.com/erikmannerfelt/ADSvalbard/blob/091b231c3eea25239ec40ad6c58155ff51826b19/adsvalbard/utilities.py#L105-L131
A modification of this may work directly in here:
# toplevel
import tempfile
import shutil
import os
from pathlib import Path
import requests # Might need to be installed.
# / toplevel
def download_large_file(url: str, filename: str | None = None, directory: Path | str | None = None) -> Path:
"""
Download a file from the requested URL.
Parameters
----------
- url:
The URL to download the file from.
- filename:
The output filename of the file. Defaults to the basename of the URL.
- directory:
The directory to save the file in. Defaults to `cache/`
Returns
-------
A path to the downloaded file.
"""
# If `directory` is defined, make sure it's a path. If it's not defined, default to `cache/`
if isinstance(directory, (str, Path)):
out_dir = Path(directory)
else:
out_dir = Path("cache/")
if filename is not None:
out_path = out_dir.joinpath(filename)
else:
out_path = out_dir.joinpath(os.path.basename(url))
# If the file already exists, skip downloading it.
if not out_path.is_file():
# Open a data stream from the URL. This means not everything has to be kept in memory.
with requests.get(url, stream=True) as request:
# Stop and raise an exception if there's a problem.
request.raise_for_status()
# Save the file to a temporary directory first. The file is constantly appended to due to the streaming
# Therefore, if the stream is cancelled, the file is not complete and should be dropped.
with tempfile.TemporaryDirectory() as temp_dir:
temp_path = Path(temp_dir).joinpath("temp.file")
# Write the data stream as it's received
with open(temp_path, "wb") as outfile:
shutil.copyfileobj(request.raw, outfile)
# When the download is done, move the file to its ultimate location.
shutil.move(temp_path, out_path)
return out_path
Note that I'm using type annotations for the arguments and return types. This is optional but really helps in documenting the function. For python versions before 3.10 (I think), the first line in the script has to be this to support type annotations:
from __future__ import annotations
But since we specified the environment.yml
to have a minimum supported python version of 3.10 (like my foreshadowing? 🤓), this is not needed. I just wanted to add this in case we run into issues with older python versions!
There are three relevant processing levels of IS2 as I understand it:
ATL03
: Photon data. A huge dataset but with the least processing that may otherwise skew our resultsATL06
: Gridded photon data: A much smaller dataset but with some processing that we may not wantATL08
: Interpolated and heavily filtered data: The smallest dataset with lots of processing that we may not want.The advantage of ATL08 (the one we currently use) is the small size, and the fact that we don't need a statistically robust approach to filter it, since filtering has already been done. I've seen some strange effects with it, however, such as mountains appearing too soon in the data (it looks like they're much wider than they really are), so this might bite us in the long term.
ATL03 would not have these issues, but would require much much more processing power for each glacier. If we can avoid this without issues, that'd be great. The problem right now, however, is that we don't know how different they are and which are more advantageous.
Maybe we can do a small case-study with all three on a few reference glaciers? For speed, we could ask Desiree or Marco to help us with the downloading and merging of different strips if that would take a long time.
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.