Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix silly overflow issue in RGI7 rgi7g_to_complex task #1720

Merged
merged 1 commit into from
Aug 6, 2024
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Fix silly overflow issue in RGI7 rgi7g_to_complex task
fmaussion committed Aug 6, 2024
commit f660f3c91198379d37f67feca1ac05a85409347b
39 changes: 20 additions & 19 deletions oggm/core/gis.py
Original file line number Diff line number Diff line change
@@ -2048,34 +2048,35 @@ def rgi7g_to_complex(gdir, rgi7g_file=None, rgi7c_to_g_links=None):

# OK all good, now the mask
# Load the DEM
# TODO: this is really really unnecessary, we should have a better way
with rasterio.open(gdir.get_filepath('dem'), 'r', driver='GTiff') as ds:
data = ds.read(1).astype(rasterio.float32)
profile = ds.profile

# Initialize the mask with -1s, the same shape as the DEM
mask = np.full(data.shape, -1, dtype=np.int16)

# Iterate over each polygon, rasterizing it onto the mask
for index, geometry in enumerate(subset.geometry):
# Correct invalid polygons
geometry = geometry.buffer(0)
if not geometry.is_valid:
raise Exception('Invalid geometry.')

# Compute the glacier mask using rasterio
# Small detour as mask only accepts DataReader objects
profile['dtype'] = 'int16'
profile.pop('nodata', None)
with rasterio.io.MemoryFile() as memfile:
with memfile.open(**profile) as dataset:
dataset.write(data.astype(np.int16)[np.newaxis, ...])
dem_data = rasterio.open(memfile.name)
# Compute the glacier mask using rasterio
# Small detour as mask only accepts DataReader objects
profile['dtype'] = 'int16'
profile.pop('nodata', None)
with rasterio.io.MemoryFile() as memfile:
with memfile.open(**profile) as dataset:
dataset.write(data.astype(np.int16)[np.newaxis, ...])
dem_data = rasterio.open(memfile.name)

# Iterate over each polygon, rasterizing it onto the mask
for index, geometry in enumerate(subset.geometry):
# Correct invalid polygons
geometry = geometry.buffer(0)
if not geometry.is_valid:
raise Exception('Invalid geometry.')
masked_dem, _ = riomask(dem_data, [shpg.mapping(geometry)],
filled=False)
glacier_mask = ~masked_dem[0, ...].mask
glacier_mask = ~masked_dem[0, ...].mask

# Update the mask: only change -1 to the new index
mask = np.where((mask == -1) & glacier_mask, index, mask)
# Update the mask: only change -1 to the new index
mask = np.where((mask == -1) & glacier_mask, index, mask)

grids_file = gdir.get_filepath('gridded_data')
with ncDataset(grids_file, 'a') as nc:
@@ -2084,7 +2085,7 @@ def rgi7g_to_complex(gdir, rgi7g_file=None, rgi7c_to_g_links=None):
if vn in nc.variables:
v = nc.variables[vn]
else:
v = nc.createVariable(vn, 'i1', ('y', 'x', ))
v = nc.createVariable(vn, 'i2', ('y', 'x', ))
v.units = '-'
v.long_name = 'Sub-entities glacier mask (number is index)'
v[:] = mask
17 changes: 17 additions & 0 deletions oggm/tests/test_workflow.py
Original file line number Diff line number Diff line change
@@ -3,6 +3,7 @@
import unittest
import pickle
import pytest
import json
import numpy as np
import xarray as xr
from numpy.testing import assert_allclose
@@ -506,6 +507,9 @@ def test_rgi7_complex_glacier_dirs():
# initialize
cfg.initialize()
cfg.PATHS['working_dir'] = _TEST_DIR
cfg.PATHS['dem_file'] = get_demo_file('hef_srtm.tif')
cfg.PARAMS['border'] = 5

# load and read test data
hef_rgi7_df = gpd.read_file(get_demo_file('rgi7c_hef.shp'))
# create GDIR
@@ -521,3 +525,16 @@ def test_rgi7_complex_glacier_dirs():
assert gdir.rgi_version == '70C'
assert gdir.rgi_dem_source is None
assert gdir.utm_zone == 32

# Also test rgi7g_to_complex
tasks.define_glacier_region(gdir)
tasks.glacier_masks(gdir)
rgi7g_file = gpd.read_file(get_demo_file('rgi7g_hef_complex.shp'))
with open(get_demo_file('hef-CtoG_links.json'), 'r') as f:
rgi7c_to_g_links = json.load(f)
tasks.rgi7g_to_complex(gdir,
rgi7g_file=rgi7g_file,
rgi7c_to_g_links=rgi7c_to_g_links)

with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds:
assert ds.sub_entities.max().item() == (len(rgi7c_to_g_links[gdir.rgi_id]) - 1)
4 changes: 3 additions & 1 deletion oggm/utils/_downloads.py
Original file line number Diff line number Diff line change
@@ -62,7 +62,7 @@
# The given commit will be downloaded from github and used as source for
# all sample data
SAMPLE_DATA_GH_REPO = 'OGGM/oggm-sample-data'
SAMPLE_DATA_COMMIT = 'c6b64892b6da37b8b43c83e2cb148e9872d34f48'
SAMPLE_DATA_COMMIT = '281ccfae3c0e82a2cae34ccb47f4ef32eb5bf6f5'

CHECKSUM_URL = 'https://cluster.klima.uni-bremen.de/data/downloads.sha256.hdf'
CHECKSUM_VALIDATION_URL = CHECKSUM_URL + '.sha256'
@@ -1917,6 +1917,8 @@ def _get_rgi_dir_unlocked(version=None, reset=False):
test_file = os.path.join(rgi_dir, 'RGI2000-v7.0-C-01_alaska', 'README.md')
dfile = 'https://daacdata.apps.nsidc.org/pub/DATASETS/nsidc0770_rgi_v7/'
dfile += 'global_files/RGI2000-v7.0-C-global.zip'
else:
raise InvalidParamsError(f'RGI version not valid: {version}')

if len(glob.glob(test_file)) == 0:
# if not there download it