-
Notifications
You must be signed in to change notification settings - Fork 0
/
bootstrap.py
76 lines (60 loc) · 2.63 KB
/
bootstrap.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
# 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 affine
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 _reproject_data_array(da, crs, extent, resolution):
"""Reproject data array to exact bounds via affine transform."""
west, east, south, north = extent
da = da.rio.clip_box(west, south, east, north, crs=crs)
bounds = da.rio.transform_bounds(crs)
xoffset = bounds[0] - bounds[0] % resolution
yoffset = bounds[1] - bounds[1] % resolution
transform = affine.Affine(resolution, 0, xoffset, 0, resolution, yoffset)
da = da.rio.reproject(crs, transform=transform, resampling=1)
da = da.rio.clip_box(west, south, east, north)
return da
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()
da = xr.open_dataarray(filepath, decode_coords='all', decode_cf=True)
# reproject to match bounds
da = _reproject_data_array(da, crs, extent, resolution)
# set better standard name
da.attrs.update(standard_name='bedrock_altitude')
# return as dataset
return da.to_dataset()