Add GEBCO elevation bootstrap data. (#51)

* Add ZipDownloader class.

* Chunk downloads to avoid overloading memory.

* Allow trailing slashes in archive url.

* Add dependency to rioxarray for reprojections.

* Add hyoga.open.surface function and module.

* Use decode_coords='all' for CF compatibility.

* Set increasing y-coord for PISM compatibility.

* Add hyoga.open.surface docstring.

* Rename hyoga.open.surface to hyoga.open.bootstrap.

* Update api docs, roadmap, whatsnew.

doc/api/open.rst

@@ -7,23 +7,26 @@ Input - ``hyoga.open``
 .. automodule:: hyoga.open
 .. currentmodule:: hyoga.open
 
-Opening example data
---------------------
+Opening local data
+------------------
 
 .. autosummary::
    :toctree: generated/
 
-   example
+   dataset
+   mfdataset
+   subdataset
+
 
-Opening local datasets
-----------------------
+Opening online data
+-------------------
 
 .. autosummary::
    :toctree: generated/
 
-   dataset
-   mfdataset
-   subdataset
+   bootstrap
+   example
+
 
 Opening vector data
 -------------------

doc/roadmap.rst

@@ -47,7 +47,7 @@ v0.3.x Inputs
 -------------
 
 - |-| :func:`hyoga.open.atmosphere`
-- |-| :func:`hyoga.open.surface`
+- |x| :func:`hyoga.open.bootstrap`
 - |-| :func:`hyoga.open.timeseries`
 
 v0.2.x Cartography

doc/whatsnew.rst

@@ -21,6 +21,17 @@
 What's new
 ==========
 
+.. _v0.3.0:
+
+v0.3.0 (unreleased)
+-------------------
+
+New features
+~~~~~~~~~~~~
+
+- Add :func:`hyoga.open.bootstrap` to open global elevation data from GEBCO as
+  bootstrapping data for PISM (:issue:`1`, :pull:`51`).
+
 .. _v0.2.2:
 
 v0.2.2 (16 Dec. 2022)

hyoga/open/__init__.py

@@ -5,12 +5,14 @@
 Hyoga input tools to open glacier modelling datasets.
 """
 
+from .bootstrap import bootstrap
 from .example import example
 from .local import dataset, mfdataset, subdataset
 from .naturalearth import natural_earth
 from .paleoglaciers import paleoglaciers
 
 __all__ = [
+    'bootstrap',
     'example',
     'dataset', 'mfdataset', 'subdataset',
     'natural_earth',

hyoga/open/bootstrap.py

@@ -0,0 +1,68 @@
+# Copyright (c) 2022, Julien Seguinot (juseg.github.io)
+# GNU General Public License v3.0+ (https://www.gnu.org/licenses/gpl-3.0.txt)
+
+"""
+This module contains a function to open surface-level bootstrapping datasets
+for predefined domains, including altitude, ice thickness and geothermal heat
+flux data.
+"""
+
+import xarray as xr
+import rioxarray  # noqa pylint: disable=unused-import
+
+import hyoga.open.downloader
+
+
+def _download_gebco():
+    """Download GEBCO sub-ice bathymetric and topographic data."""
+    downloader = hyoga.open.downloader.ZipDownloader()
+    filepath = downloader(
+        'https://www.bodc.ac.uk/data/open_download/gebco/'
+        'gebco_2022_sub_ice_topo/zip/', 'gebco/GEBCO_2022_sub_ice_topo.nc')
+    return filepath
+
+
+def bootstrap(crs, extent, resolution=1e3):
+    """
+    Open bootstrapping data from online datasets for PISM.
+
+    Currently a single dataset (GEBCO) is supported.
+
+    Parameters
+    ----------
+    crs : str
+        Coordinate reference system for the resulting dataset as OGC WKT or
+        Proj.4 string, will be passed to Dataset.rio.reproject.
+    extent : (west, east, south, north)
+        Extent for the resulting dataset in projected coordinates given by
+        ``crs``, will be passed to Dataset.rio.clip_box.
+    resolution : float, optional
+        Resolution for the output dataset in projected coordinates given by
+        ``crs``, will be passed to Dataset.rio.reproject.
+
+    Returns
+    -------
+    ds : Dataset
+        The resulting dataset containing surface variables with the requested
+        ``crs``, ``extent``, and ``resolution``. Use ``ds.to_netcdf()`` to
+        export as PISM bootstrapping file.
+    """
+
+    # open global data (use decode_coords='all' to read grid_mapping attribute)
+    filepath = _download_gebco()
+    ds = xr.open_dataset(filepath, decode_coords='all')
+
+    # clip, reproject and clip again
+    west, east, south, north = extent
+    ds = ds.rio.clip_box(west, south, east, north, crs=crs)
+    ds = ds.rio.reproject(crs, resolution=resolution)
+    ds = ds.rio.clip_box(west, south, east, north)
+
+    # flip data to increasing y-coord (needed by PISM)
+    ds = ds.isel(y=slice(None, None, -1))
+
+    # set better standard name
+    ds.elevation.attrs.update(standard_name='bedrock_altitude')
+
+    # return projected dataset
+    return ds

hyoga/open/downloader.py

@@ -73,10 +73,15 @@ def get(self, url, path):
         # create directory if missing
         os.makedirs(os.path.dirname(path), exist_ok=True)
 
-        # download file
-        with open(path, 'wb') as binaryfile:
+        # open url and raise any http error
+        with requests.get(url, stream=True, timeout=5) as request:
+            request.raise_for_status()
+
+            # download file chunks
             print(f"downloading {url}...")
-            binaryfile.write(requests.get(url, timeout=5).content)
+            with open(path, 'wb') as binaryfile:
+                for chunk in request.iter_content(chunk_size=1024**2):
+                    binaryfile.write(chunk)
 
 
 class CacheDownloader(Downloader):
@@ -130,7 +135,8 @@ def get(self, url, path, member=None):
 
         # save archive as named online
         outdir, basename = os.path.split(path)
-        archivepath = os.path.join(outdir, url.split('/')[-1])
+        # FIXME GEBCO arvhive is just named 'zip'
+        archivepath = os.path.join(outdir, url.rstrip('/').split('/')[-1])
 
         # download it only if missing
         if not super().check(archivepath):
@@ -147,8 +153,27 @@ def deflate(self, archivepath, member, outdir):
         raise NotImplementedError("This should be implemented in subclasses.")
 
 
+class ZipDownloader(ArchiveDownloader):
+    """
+    Download a zip archive and extract a single file.
+
+    Call parameters
+    ---------------
+    url : str
+        The url of the file to download
+    path : str
+        The path of the extracted file relative to the cache directory.
+    member : str, optional
+        Member file to extract from , default to the basename of ``path``.
+    """
+
+    def deflate(self, archivepath, member, outdir):
+        with zipfile.ZipFile(archivepath, 'r') as archive:
+            archive.extract(member, path=outdir)
+
+
 class ShapeZipDownloader(ArchiveDownloader):
-    """A downloader that extract shapefiles and metafiles from zip archives.
+    """A downloader that extracts shapefiles and metafiles from zip archives.
 
     Call parameters
     ---------------

setup.cfg

@@ -26,6 +26,7 @@ install_requires =
     geopandas
     matplotlib
     requests
+    rioxarray
     scipy
     xarray
 python_requires = >=3.8