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

See what happens if we default to utm #1734

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
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
2 changes: 1 addition & 1 deletion oggm/params.cfg
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,7 @@ lru_maxsize = 100
### Glacier map and grid determination

# Map projection: 'tmerc' or 'utm'
map_proj = 'tmerc'
map_proj = 'utm'

# Decision on grid spatial resolution for each glacier
# 'fixed': dx (meters) = fixed_dx
Expand Down
3 changes: 3 additions & 0 deletions oggm/tests/funcs.py
Original file line number Diff line number Diff line change
Expand Up @@ -409,6 +409,7 @@ def init_hef(reset=False, border=40, logging_level='INFO', rgi_id=None,
cfg.PARAMS['prcp_fac'] = 2.5
cfg.PARAMS['temp_bias_min'] = -10
cfg.PARAMS['temp_bias_max'] = 10
cfg.PARAMS['map_proj'] = 'tmerc'
hef_file = get_demo_file('Hintereisferner_RGI5.shp')
entity = gpd.read_file(hef_file).iloc[0]

Expand Down Expand Up @@ -515,6 +516,7 @@ def init_columbia(reset=False):
cfg.PARAMS['prcp_fac'] = 2.5
cfg.PARAMS['baseline_climate'] = 'CRU'
cfg.PARAMS['evolution_model'] = 'FluxBased'
cfg.PARAMS['map_proj'] = 'tmerc'

entity = gpd.read_file(get_demo_file('01_rgi60_Columbia.shp')).iloc[0]
gdir = oggm.GlacierDirectory(entity, reset=reset)
Expand Down Expand Up @@ -556,6 +558,7 @@ def init_columbia_eb(dir_name, reset=False):
cfg.PARAMS['prcp_fac'] = 2.5
cfg.PARAMS['baseline_climate'] = 'CRU'
cfg.PARAMS['evolution_model'] = 'FluxBased'
cfg.PARAMS['map_proj'] = 'tmerc'

entity = gpd.read_file(get_demo_file('01_rgi60_Columbia.shp')).iloc[0]
gdir = oggm.GlacierDirectory(entity)
Expand Down
2 changes: 1 addition & 1 deletion oggm/tests/test_benchmarks.py
Original file line number Diff line number Diff line change
Expand Up @@ -383,7 +383,7 @@ def to_optimize(x):
ds['ref'].data[df['j'], df['i']] = df['thick']

rmsd = ((df.oggm - df.thick) ** 2).mean() ** .5
assert rmsd < 30
assert rmsd < 32

dfm = df.mean()
np.testing.assert_allclose(dfm.thick, dfm.oggm, 10)
Expand Down
5 changes: 5 additions & 0 deletions oggm/tests/test_graphics.py
Original file line number Diff line number Diff line change
Expand Up @@ -182,6 +182,7 @@ def test_multiple_inversion():
cfg.PARAMS['use_winter_prcp_fac'] = False
cfg.PARAMS['use_temp_bias_from_file'] = False
cfg.PARAMS['prcp_fac'] = 2.5
cfg.PARAMS['map_proj'] = 'tmerc'
cfg.PATHS['working_dir'] = testdir

# Get the RGI ID
Expand Down Expand Up @@ -269,6 +270,7 @@ def test_multiple_models():
cfg.PARAMS['use_temp_bias_from_file'] = False
cfg.PARAMS['prcp_fac'] = 2.5
cfg.PARAMS['border'] = 40
cfg.PARAMS['map_proj'] = 'tmerc'

# Get the RGI ID
hef_rgi = gpd.read_file(get_demo_file('divides_hef.shp'))
Expand Down Expand Up @@ -379,6 +381,7 @@ def test_chhota_shigri():
cfg.PARAMS['use_winter_prcp_fac'] = False
cfg.PARAMS['use_temp_bias_from_file'] = False
cfg.PARAMS['prcp_fac'] = 2.5
cfg.PARAMS['map_proj'] = 'tmerc'

hef_file = get_demo_file('divides_RGI50-14.15990.shp')
df = gpd.read_file(hef_file)
Expand Down Expand Up @@ -424,6 +427,7 @@ def test_ice_cap():
cfg.PARAMS['use_winter_prcp_fac'] = False
cfg.PARAMS['use_temp_bias_from_file'] = False
cfg.PARAMS['prcp_fac'] = 2.5
cfg.PARAMS['map_proj'] = 'tmerc'

df = gpd.read_file(get_demo_file('divides_RGI50-05.08389.shp'))
df['Area'] = df.Area * 1e-6 # cause it was in m2
Expand Down Expand Up @@ -467,6 +471,7 @@ def test_coxe():
cfg.PARAMS['use_winter_prcp_fac'] = False
cfg.PARAMS['use_temp_bias_from_file'] = False
cfg.PARAMS['prcp_fac'] = 2.5
cfg.PARAMS['map_proj'] = 'tmerc'

hef_file = get_demo_file('rgi_RGI50-01.10299.shp')
entity = gpd.read_file(hef_file).iloc[0]
Expand Down
40 changes: 22 additions & 18 deletions oggm/tests/test_prepro.py
Original file line number Diff line number Diff line change
Expand Up @@ -121,7 +121,7 @@ def test_define_region(self):
# From string
gdir = oggm.GlacierDirectory(gdir.rgi_id, base_dir=self.testdir)
# This is not guaranteed to be equal because of projection issues
np.testing.assert_allclose(extent, gdir.extent_ll, atol=1e-5)
np.testing.assert_allclose(extent, gdir.extent_ll, atol=5e-4)
with warnings.catch_warnings():
warnings.filterwarnings("ignore", category=FutureWarning)
# Warning in salem
Expand Down Expand Up @@ -263,7 +263,7 @@ def test_repr(self):
Status: Glacier or ice cap
Area: 8.036 km2
Lon, Lat: (10.7584, 46.8003)
Grid (nx, ny): (159, 114)
Grid (nx, ny): (158, 116)
Grid (dx, dy): (50.0, -50.0)
""")

Expand Down Expand Up @@ -368,11 +368,12 @@ def test_compute_hypsometry_attributes(self):

np.testing.assert_allclose(dfh['slope_deg'], entity.Slope, atol=0.5)
np.testing.assert_allclose(dfh['aspect_deg'], entity.Aspect, atol=5)
np.testing.assert_allclose(dfh['zmed_m'], entity.Zmed, atol=20)
np.testing.assert_allclose(dfh['zmax_m'], entity.Zmax, atol=20)
np.testing.assert_allclose(dfh['zmin_m'], entity.Zmin, atol=20)
np.testing.assert_allclose(dfh['zmax_m'], entity.Zmax, atol=20)
np.testing.assert_allclose(dfh['zmin_m'], entity.Zmin, atol=20)
atol = 25
np.testing.assert_allclose(dfh['zmed_m'], entity.Zmed, atol=atol)
np.testing.assert_allclose(dfh['zmax_m'], entity.Zmax, atol=atol)
np.testing.assert_allclose(dfh['zmin_m'], entity.Zmin, atol=atol)
np.testing.assert_allclose(dfh['zmax_m'], entity.Zmax, atol=atol)
np.testing.assert_allclose(dfh['zmin_m'], entity.Zmin, atol=atol)
# From google map checks
np.testing.assert_allclose(dfh['terminus_lon'], 10.80, atol=0.01)
np.testing.assert_allclose(dfh['terminus_lat'], 46.81, atol=0.01)
Expand Down Expand Up @@ -427,6 +428,7 @@ def test_rasterio_glacier_masks(self):
hef_file = get_demo_file('Hintereisferner_RGI5.shp')
entity = gpd.read_file(hef_file).iloc[0]

cfg.PARAMS['map_proj'] = 'tmerc'
gdir = oggm.GlacierDirectory(entity, base_dir=self.testdir)
gis.define_glacier_region(gdir)

Expand Down Expand Up @@ -1160,7 +1162,7 @@ def test_width(self):
50.)
h1, b = np.histogram(hgt, weights=harea, density=True, bins=bins)
h2, b = np.histogram(rhgt, density=True, bins=bins)
assert utils.rmsd(h1*100*50, h2*100*50) < 1
assert utils.rmsd(h1*100*50, h2*100*50) < 1.1

# Check that utility function is doing what is expected
hh, ww = gdir.get_inversion_flowline_hw()
Expand Down Expand Up @@ -1915,8 +1917,8 @@ def to_optimize(x):
tol=1e-4)['x']
self.assertTrue(out[0] > 0.1)
self.assertTrue(out[1] > 0.1)
self.assertTrue(out[0] < 1.1)
self.assertTrue(out[1] < 1.1)
self.assertTrue(out[0] < 1.2)
self.assertTrue(out[1] < 1.2)
glen_a = glen_a * out[0]
fs = fs * out[1]
v = inversion.mass_conservation_inversion(gdir, fs=fs,
Expand Down Expand Up @@ -1998,7 +2000,7 @@ def test_invert_hef_from_consensus(self):
df = workflow.calibrate_inversion_from_consensus(gdir,
a_bounds=a,
error_on_mismatch=False)
np.testing.assert_allclose(df.vol_itmix_m3, df.vol_oggm_m3, rtol=0.08)
np.testing.assert_allclose(df.vol_itmix_m3, df.vol_oggm_m3, rtol=0.1)

# With fs it can work
a = (0.1, 3)
Expand Down Expand Up @@ -2217,22 +2219,24 @@ def test_distribute(self):
inversion.prepare_for_inversion(gdir)

ref_v = 0.573 * 1e9
cfg.PARAMS['inversion_fs'] = 5.7e-20

def to_optimize(x):
glen_a = cfg.PARAMS['inversion_glen_a'] * x[0]
fs = cfg.PARAMS['inversion_fs'] * x[1]
v = inversion.mass_conservation_inversion(gdir, fs=fs,
glen_a=glen_a)
glen_a=glen_a)
return (v - ref_v)**2

import scipy.optimize as optimization
out = optimization.minimize(to_optimize, [1, 1],
bounds=((0.01, 10), (0.01, 10)),
tol=1e-1)['x']
glen_a = cfg.PARAMS['inversion_glen_a'] * out[0]
fs = cfg.PARAMS['inversion_fs'] * out[1]
v = inversion.mass_conservation_inversion(gdir, fs=fs,
glen_a=glen_a,
write=True)
glen_a=glen_a,
write=True)
np.testing.assert_allclose(ref_v, v)

inversion.distribute_thickness_interp(gdir, varname_suffix='_interp')
Expand Down Expand Up @@ -2281,7 +2285,7 @@ def to_optimize(x):
glen_a = cfg.PARAMS['inversion_glen_a'] * x[0]
fs = 0.
v = inversion.mass_conservation_inversion(gdir, fs=fs,
glen_a=glen_a)
glen_a=glen_a)
return (v - ref_v)**2

import scipy.optimize as optimization
Expand All @@ -2290,13 +2294,13 @@ def to_optimize(x):
tol=1e-4)['x']

self.assertTrue(out[0] > 0.1)
self.assertTrue(out[0] < 10)
self.assertTrue(out[0] < 15)

glen_a = cfg.PARAMS['inversion_glen_a'] * out[0]
fs = 0.
v = inversion.mass_conservation_inversion(gdir, fs=fs,
glen_a=glen_a,
write=True)
glen_a=glen_a,
write=True)
np.testing.assert_allclose(ref_v, v)

cls = gdir.read_pickle('inversion_output')
Expand Down
6 changes: 3 additions & 3 deletions oggm/tests/test_shop.py
Original file line number Diff line number Diff line change
Expand Up @@ -138,7 +138,7 @@ def test_thickvel_to_glacier(self, class_case_dir, monkeypatch):
thick = ds.millan_ice_thickness.where(mask).data

# Simply some coverage and sanity checks
assert np.isfinite(thick).sum() / mask.sum() > 0.98
assert np.isfinite(thick).sum() / mask.sum() > 0.97
assert np.nanmax(thick) > 800
assert np.nansum(thick) * gdir.grid.dx**2 * 1e-9 > 174

Expand All @@ -152,14 +152,14 @@ def test_thickvel_to_glacier(self, class_case_dir, monkeypatch):
vy = ds.millan_vy.where(mask).data

# Simply some coverage and sanity checks
assert np.isfinite(v).sum() / mask.sum() > 0.98
assert np.isfinite(v).sum() / mask.sum() > 0.97
assert np.nanmax(v) > 2000
assert np.nanmax(vx) > 500
assert np.nanmax(vy) > 400

# Stats
df = millan22.compile_millan_statistics([gdir]).iloc[0]
assert df['millan_avg_vel'] > 180
assert df['millan_avg_vel'] > 175
assert df['millan_max_vel'] > 2000
assert df['millan_perc_cov'] > 0.95
assert df['millan_vel_perc_cov'] > 0.97
Expand Down
1 change: 1 addition & 0 deletions oggm/tests/test_utils.py
Original file line number Diff line number Diff line change
Expand Up @@ -1686,6 +1686,7 @@ def test_start_from_prepro(self):
'baseline_climate': 'CRU',
'use_winter_prcp_fac': False,
'use_temp_bias_from_file': False,
'map_proj': 'tmerc',
'prcp_fac': 2.5,
}
)
Expand Down
1 change: 1 addition & 0 deletions oggm/tests/test_workflow.py
Original file line number Diff line number Diff line change
Expand Up @@ -103,6 +103,7 @@ def up_to_climate(reset=False, use_mp=None):
cfg.PARAMS['baseline_climate'] = 'CRU'
cfg.PARAMS['evolution_model'] = 'FluxBased'
cfg.PARAMS['downstream_line_shape'] = 'parabola'
cfg.PARAMS['map_proj'] = 'tmerc'

# Go
gdirs = workflow.init_glacier_directories(rgidf)
Expand Down
Loading