-
Notifications
You must be signed in to change notification settings - Fork 108
/
Copy pathrgitopo.py
165 lines (135 loc) · 5.58 KB
/
rgitopo.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
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
import logging
import os
import shutil
import warnings
try:
import salem
except ImportError:
pass
try:
import rasterio
except ImportError:
pass
import numpy as np
from oggm import utils, workflow
from oggm.exceptions import InvalidParamsError
# Module logger
log = logging.getLogger(__name__)
DEMS_URL = 'https://cluster.klima.uni-bremen.de/data/gdirs/dems_v2/default'
DEMS_HR_URL = 'https://cluster.klima.uni-bremen.de/data/gdirs/dems_v1/highres/'
def init_glacier_directories_from_rgitopo(rgidf=None, dem_source=None,
resolution='default',
keep_dem_folders=False,
reset=False,
force=True):
"""Initialize a glacier directory from an RGI-TOPO DEM.
Wrapper around :func:`workflow.init_glacier_directories`, which selects
a default source for you (if not set manually with `dem_source`).
The default source is NASADEM for all latitudes between 60N and 56S,
and COPDEM90 elsewhere.
Parameters
----------
rgidf : GeoDataFrame or list of ids, optional for pre-computed runs
the RGI glacier outlines. If unavailable, OGGM will parse the
information from the glacier directories found in the working
directory. It is required for new runs.
reset : bool
delete the existing glacier directories if found.
force : bool
setting `reset=True` will trigger a yes/no question to the user. Set
`force=True` to avoid this.
dem_source : str
the source to pick from (default: NASADEM and COPDEM90)
keep_dem_folders : bool
the default is to delete the other DEM directories to save space.
Set this to True to prevent that (e.g. for sensitivity tests)
Returns
-------
gdirs : list of :py:class:`oggm.GlacierDirectory` objects
the initialised glacier directories
"""
if resolution == 'default':
base_url = DEMS_URL
elif resolution == 'lr':
base_url = DEMS_URL
elif resolution == 'hr':
log.warning('High-res DEMs not available yet in version 2 with COPDEM30')
base_url = DEMS_HR_URL
else:
raise InvalidParamsError('`resolution` should be of `lr` or `hr`')
gdirs = workflow.init_glacier_directories(rgidf, reset=reset, force=force,
prepro_base_url=base_url,
from_prepro_level=1,
prepro_rgi_version='62')
workflow.execute_entity_task(select_dem_from_dir, gdirs,
dem_source=dem_source,
keep_dem_folders=keep_dem_folders)
return gdirs
@utils.entity_task(log, writes=['dem'])
def select_dem_from_dir(gdir, dem_source=None, keep_dem_folders=False):
"""Select a DEM from the ones available in an RGI-TOPO glacier directory.
Parameters
----------
gdir : :py:class:`oggm.GlacierDirectory`
the glacier directory
dem_source : str
the source to pick from. If 'RGI', we assume that there is a
`dem_source` attribute in the RGI file. If 'BY_RES', we use
COPDEM30 for all gdirs with resolution smaller than 60m
keep_dem_folders : bool
the default is to delete the other DEM directories to save space.
Set this to True to prevent that (e.g. for sensitivity tests)
"""
if dem_source == 'RGI':
dem_source = gdir.rgi_dem_source
if dem_source == 'BY_RES':
dem_source = 'COPDEM30' if gdir.grid.dx < 60 else 'COPDEM90'
sources = [f.name for f in os.scandir(gdir.dir) if f.is_dir()
and not f.name.startswith('.')]
if dem_source not in sources:
raise RuntimeError('source {} not in folder'.format(dem_source))
gdir.add_to_diagnostics('dem_source', dem_source)
shutil.copyfile(os.path.join(gdir.dir, dem_source, 'dem.tif'),
gdir.get_filepath('dem'))
shutil.copyfile(os.path.join(gdir.dir, dem_source, 'dem_source.txt'),
gdir.get_filepath('dem_source'))
if not keep_dem_folders:
for source in sources:
shutil.rmtree(os.path.join(gdir.dir, source))
def _fallback_dem_quality_check(gdir):
return dict(rgi_id=gdir.rgi_id)
@utils.entity_task(log, writes=[], fallback=_fallback_dem_quality_check)
def dem_quality_check(gdir):
"""Run a simple quality check on the rgitopo DEMs
Parameters
----------
gdir : GlacierDirectory
the glacier directory
Returns
-------
a dict of DEMSOURCE:frac pairs, where frac is the percentage of
valid DEM grid points on the glacier.
"""
with rasterio.open(gdir.get_filepath('glacier_mask'), 'r',
driver='GTiff') as ds:
mask = ds.read(1) == 1
area = mask.sum()
sources = [f.name for f in os.scandir(gdir.dir) if f.is_dir()
and not f.name.startswith('.')]
out = dict(rgi_id=gdir.rgi_id)
for s in sources:
try:
with rasterio.open(os.path.join(gdir.dir, s, 'dem.tif'), 'r',
driver='GTiff') as ds:
topo = ds.read(1).astype(rasterio.float32)
topo[topo <= -999.] = np.NaN
topo[ds.read_masks(1) == 0] = np.NaN
valid_mask = np.isfinite(topo) & mask
if np.all(~valid_mask):
continue
if np.nanmax(topo[valid_mask]) == 0:
continue
out[s] = valid_mask.sum() / area
except BaseException:
pass
return out