juseg / hyoga Goto Github PK

Allow the user to configure the location of hyoga's cache directory (currently $XDG_CACHE_HOME/hyoga).

Add annotation helpers from cartowik, improve them, and document them. Facilitate accurate placement of matplotlib annotations labels, especially on geographic axes.

• Add annotations with auto label anchor.
• Add annotations by angle and offset.
• Add annotations example.
• Add annotated scatter plot function.
• Add annotated scatter plot example.

• Set up documentation with Sphinx
• Complete docstrings in hyoga.hyoga
• Complete docstrings in hyoga.open
• Complete docstrings in hyoga.plot
• Write a page on opening model output
• Write a page on plotting model output
• Upload to readthedocs.io.

Something like pismx.boot.etopo1().to_netcdf() hyoga.open.surface().to_netcdf() hyoga.open.bootstrap().to_netcdf().

Downloads from multiple sources and potentially formats will be needed in v0.3.0 onward. I think the best way to solve this is object-oriented. With this issue I intend to get a head start on a new downloader, and subclass it for Natural Earth data, so that cartopy is no longer required. I plan to do it in v0.2.x and thereby get readier for v0.3.0.

• Add generic downloader class.
• Add Natural Earth downloader.
• Remove cartopy dependence.

Most examples contain the following code.

ax.xaxis.set_visible(False)                                                 
ax.yaxis.set_visible(False)

I think we want to hide axes by default instead. The x-coordinate... labels are typically projection coordinates which are hard to read anyway. A scale bar (#16) will be more useful.

Interpolation to axes coordinate with ds.hyoga.interp(ax=ax) fails on matplotlib Axes due to an internal call of ax.get_extent() which is only available on cartopy Geoaxes. To reproduce:

import matplotlib.pyplot as plot
import hyoga

ax = plt.axes()
with hyoga.open.example('pism.alps.out.2d.nc') as ds:
    ds = ds.hyoga.interp(hyoga.open.example('pism.alps.vis.refined.nc'), ax=ax)

Readthedocs builds to fail since latest scipy release.

This is related to scipy/scipy#17718, only affecting scipy==1.10.0 (will be fixed in 1.10.1).

Profile interpolation crashes on mixed precision, for instance mixing float32 PISM output data and float64 profile coordinates. While not an ideal solution for large arrays, the bug can be worked around by upcasting everything to float64.

Add Natural Earth functions from cartowik. This is essentially a thin wrapper around cartopy, adding a little convenience, and a lot of speed in some cases.

• Add shapefile partial plotting and docstring.
• Add Natural Earth functions and docstrings.
• Use annotated scatter for custom cities maps.
• Update examples, cartopy section in docs.
• New doc page on adding map elements.

I have been working on a logo during early development on Anafi, time to use it.

• Commit old Anafi logo script.
• Decide main layout, colours.
• Other layouts for favicon etc.

Make dataset plot methods Dataset.hyoga.plot.method(self, ...) available as functions hyoga.plot.function(ds, ...).

• Add dataset plot functions in hyoga.plot.datasets.
• Add dataset calc functions in hyoga.calc.interp etc.

I was looking for a way to combine shapefile geometries in a single artist, decided to give a quick try to GeoPandas and it blew my mind. Here is a quick speed test merging Ehlers et al. LGM global and alps paleoglaciers, removing duplicates, and plotting in Orthographic projection (as in the logo) with:

• Cartopy BasicReader in 2m41s,
• Cartopy FionaReader in 2m45s,
• GeoPandas read_file in 5s!

Besides the increase in speed, GeoPandas will simplify hyoga internals and provide users with more flexible objects for further manipulation or plotting styles. However this will also introduces changes in API. To avoid future breaking change or deprecation warnings, I think I want to get this right before v0.2.

• Use GeoPandas for Natural Earth.
• Use GeoPandas for Paleoglaciers.
• Add hyoga.open.naturalearth?
• Add hyoga.open.paleoglaciers?
• Update setup.cfg, doc/environment.yml.

Add optional atmospheric buffer to prevent ice flow across domain boundaries. Related to juseg/r.out.pism@377de4d and pism/pism#279.

Add five colormaps from cartowik, document them, and enable them in plot methods.

• Add altitude and relief shading colormaps.
• Enable new colormaps in plot methods as imshow.
• Add color table example as in matplotlib.
• Add colormap reference as in matplotlib.
• Document new colormaps on plotting page.

Currently aspect ratios are not set and they are disturbed with every new image plot. Adding a naturalearth() overlay hides this problem, because GeoPandas sets the aspect ratio to 'equal'. Xarray provides an aspect keyword, but this is not satisfactory because it also requires setting size and can't be used with the ax keyword. We probably need wrappers around imshow, contour etc. But at least in the case imshow, ax.set_aspect needs to come after the call to imshow, something like:

def _imshow(self, **kwargs):
    img = var.plot.imshow(**style)
    kwargs.get('ax', plt.gca()).set_aspect('equal')
    return img

I think the wrappers could also hide coordinate ticks and labels, and handle colorbars.

Add shaded relief functions from cartowik, and enable hillshades in either existing, or new, accessor plot methods. Try to keep reasonably low numbers of arguments in plot method.

• Add hillshade functions from cartowik and docstrings.
• Rework *_altitude plot methods to include hillshade, or
• Add {bedrock,surface}_hillshade plot methods.
• Add mat and shiny shaded relief examples.
• Include a shaded relief on plotting page.

I found out that typing hyoga all the time is not so practical (at least on a Swiss keyboard). I am thinking to rename the dataset accessor as ice which is certainly more memorable from v0.2.0. The name may be used in other projects though, I am not sure. The name hyoga would remain usable until at least v0.4.0. If you are watching and have any thoughts about that, feel free to share your opinion.

Something like pismx.boot.srtm().to_netcdf().

Currently masking is achieved using xarray.DataArray.where (in hyoga.hyoga.HyogaDataset.where). Because it uses nans, one drawback with this approach is that there is no "turning back" masked values into non-masked values. This is an issue at least for hyoga.hyoga.HyogaDataset.interp, which needs to refine the mask. An alternative or perhaps complementary masking method would be to use a separate variable (standard name land_ice_area_fraction) and to make sure plotting methods use this variable if it is present (see commit message in 4120b79). Example method names:

• Add hyoga.hyoga.HyogaDataset.assign_icemask.
• Add hyoga.hyoga.HyogaDataset.where_icemask.
• Plotting methods will use the latter xarray.DataArray.where.

A few more tasks before closing this issue.

• Add hyoga.hyoga.HyogaDataset.assign.
• Use new mask variable in interp method.
• Add configurable thickness threshold.
• Document changes in whatsnew.
• Update plot examples.

Add class for modelling domains with coordinate reference system, bounds and resolution.

This could be done with pylint (which I use locally) or black (used by xarray) or both.

Add basic tests for plot methods. This should be quite straightforward now. If plot methods work, there is a good chance the rest of the package works.

Store files currently stored in root of cache directory (example and paleoglacier data) into subdirectories (just like Natural Earth data). Add code to move old files to their new location so we don't pollute the users' cache with duplicates.

Generic scripts:

• Plotting with xarray
• Plotting with cartopy
• Plotting with cartowik

Hyoga plot presets:

• Bedrock uplift preset
• Glacier erosion preset
• Surface velocity preset
• Surface streamplot preset

I saw your tweet and this looks like a cool project.

cf_xarray should help you remove the CF parsing code, it's intended to allow indexing by standard_name etc. I'd be happy to consider changes that would make its use in this project easier.

Add scale bar method from cartowik.

• Add scale_bar function and docstring.
• Remove subfig_labels, will be moved to absplots.
• Add scale bars in docs, example(s).

• Add method to mask glacier variables.
• Add method to compute bedrock uplift.
• Clean-up method to interpolate data.

With a small change we can make CHELSA input elevation data an option for bootstrapping. This would reduce the need for a separate elevation dataset (i.e. GEBCO) in case the user does not need bathymetry.

• Add bedrock, surface kwargs to bootstrap.
• Add 'chelsa', 'gebco' options for both.

Add default set of modelling domains for glaciated mountain ranges worldwide.

Add profile coordinate builder from cartowik, and implement profile interpolation in new accessor method.

• Add profile tools from cartowik and docstrings.
• Add accessor method ds.hyoga.profile.
• Add demo profile coords, example(s).
• Add docs, maybe alongside interp.
• Consider 3D demo data if not too heavy.

Not all hyoga functionality can be demonstrated with the Zenodo data. Besides, these data are already post-processed and thus not ideal to show what happens when opening a raw PISM (or other models) file. Finally, the file is very large, using a lot of resources on read the docs and for anyone who would like to use the demo data. Create and use a separate repository containing:

• PISM demo boot file
• PISM demo series file
• PISM demo gridded file
• PISM demo interp file

When this is done new examples will be added to

• Plot bedrock uplift
• Plot interpolated data

Merge pism-palwrapper into pismx.

Hyoga circumvents auto-scaling by freezing axes bounds when plotting global paleoglacier and Natural Earth vectors. A better approach will be to leave auto-scaling on but exclude the vectors using autolim or is_autoscale in future versions of matplotlib and geopandas:

• geopandas/geopandas#2602
• geopandas/geopandas#2817
• matplotlib/matplotlib#24416
• matplotlib/matplotlib#15595

This is the problematic code:

Potential candidates:

• Cleator et al., 2020: partial land pollen reconstruction.
• PMIP4 LGM experiments: ten models, no ensemble means.
• CHELSA-PMIP5: seven models, 1km downscaled.
• CHELSA-Trace21k: transient output, 1km downscaled, every 100 yr.

Add functions to download, cache, and easily plot global (Last Glacial Maximum) paleoglacier extent datasets.

• Support plotting broken shapefiles.
• Add Ehlers & al paleoglaciers.
• Add Batchelor & al paleoglaciers.

Something like pismx.atm[osphere].worldclim().to_netcdf() hyoga.open.atmosphere().to_netcdf().

Add type hints, and type checking with mypy, possibly through continuous integration.

Add a class for domain collections with methods to e.g. set their variable resolution or plot themselves on a map, or add such functionality as part of the domain class.

Hyoga v0.3.0 will use rioxarray for geographic transformations. Vector drawing methods currently only recognize a proj4 attribute, but instead should be compatible with rioxarray CRS. Rioxarray's documentation recommends using decode_coords="all" to read mapping attributes in netCDF.

To reproduce (in gebco branch for now, but soon in main too):

ds = hyoga.open.surface('epsg:32632', [150e3, 1050e3, 4820e3, 5420e3])
ds.hyoga.plot.natural_earth()

The new code produces bootstrapping data where the x and y coordinate values differ from the requested extent. They appear to align instead to the somewhat random coordinates of the NE corner of the clip_box prior to reprojection, which depend on the original grid, and the user-input extent.

This is a problem as data from different sources will project to different grids (e.g. bootstrapping and atmospheric data). I think the workaround is to provide ds.rio.reproject directly with a transform argument instead of crs and resolution.

I realize the documentation are dotted with center=False. This is because xarray defaults to centering any image plot where data span negative and positive values, however this is not really appropriate for glacier-depressed bedrock elevation typically ranging from slightly below sea level to the mountain tops. I can think of two ways to fix this: use sealevel=-120 and possibly the 'Elevational' instead of 'Topographic' colormap, or changing bedrock_altitude() and friends' default behaviour.

• Implement better defaults for altitude?
• Remove center=False in docs.
• Doc changes in whatsnew.

As far as I understand Lucazeau (2019) provides an update on the similarity method by Goutorbe et al. (2011) which I have used on the Alps, including newer point data and increased resolution of 0.5°. The IHFC provides a nice overview on recent studies.

I just realized Dataset.hyoga.plot.bedrock_altitude_contours is still missing a proper docstring.

Current syntax and API reference might be a bit bulky. I am considering to import hyoga.open input functions and the main class HyogaDataset at top level. Perhaps plotting methods can be then documented as HyogaDataset.plot.method.

• Top-level imports
• Update api
• Update examples
• Update documentation

The result is going to look ridiculous but I think it would be nice to have it in place early on.

The new downloader only checks for the .shp file. This is a problem for at least Ehlers et al. data as hyoga<=2.0 did not download the .prj file.

So we don't have to run twine upload for every release.

Explain hyoga's development workflow and add instructions for issues and pull requests.

Merge pism-palseries into pismx and cleanup.

I am opening this documentation-only issue to avoid overloading #1 and #3.

• Document opening online datasets.
• Add a minimal example for PISM.
• Perhaps rework datasets chapter outline.

I am thinking to replace opening and masking pages with:

• Opening datasets
• Masking and selection
• Spatial interpolation