From 40d893c34416b82eae13a57972ef8fc8b3cef5f6 Mon Sep 17 00:00:00 2001 From: mdtanker Date: Fri, 4 Nov 2022 08:14:49 +1300 Subject: [PATCH 1/3] diagnosing bedmap2 issues --- antarctic_plots/fetch.py | 228 +++++++++++++++++++-------------------- 1 file changed, 113 insertions(+), 115 deletions(-) diff --git a/antarctic_plots/fetch.py b/antarctic_plots/fetch.py index 27d95b46..953d72aa 100644 --- a/antarctic_plots/fetch.py +++ b/antarctic_plots/fetch.py @@ -1323,14 +1323,14 @@ def bedmap_points( def bedmap2( layer: str, - reference: str = "geoid", - plot: bool = False, - info: bool = False, - region=None, - spacing=None, - registration=None, - fill_nans=False, - **kwargs, + # reference: str = "geoid", + # plot: bool = False, + # info: bool = False, + # region=None, + # spacing=None, + # registration=None, + # fill_nans=False, + # **kwargs, ) -> xr.DataArray: """ Load bedmap2 data. All grids are by default referenced to the g104c geoid. Use the @@ -1362,13 +1362,13 @@ def bedmap2( xr.DataArray Returns a loaded, and optional clip/resampled grid of Bedmap2. """ - # Declare initial grid values, - # use utils.get_grid_info(xr.load_dataarray(file).squeeze()) + # Declare initial grid values, of .nc files not .tiff files + # use utils.get_grid_info(xr.load_dataset(file).band_data # several of the layers have different values if layer == "lakemask_vostok": - initial_region = [1189500, 1470500, -401500, -291500] + initial_region = [1190000.0, 1470000.0, -402000.0, -291000.0] initial_spacing = 1e3 - initial_registration = "p" + initial_registration = "g" elif layer == "thickness_uncertainty_5km": initial_region = [-3401500, 3403500, -3397500, 3397500] @@ -1376,17 +1376,16 @@ def bedmap2( initial_registration = "p" else: - # y lims are differnt if inputting fname straight to utils.get_grid_info() - initial_region = [-3333500, 3333500, -3332500, 3332500] + initial_region = [-3333000, 3333000, -3333000, 3333000] initial_spacing = 1e3 - initial_registration = "p" + initial_registration = "g" - if region is None: - region = initial_region - if spacing is None: - spacing = initial_spacing - if registration is None: - registration = initial_registration + # if region is None: + # region = initial_region + # if spacing is None: + # spacing = initial_spacing + # if registration is None: + # registration = initial_registration # retrieve the specified layer file path = pooch.retrieve( @@ -1397,109 +1396,108 @@ def bedmap2( processor=pooch.Unzip(), progressbar=True, ) +# "https://ramadda.data.bas.ac.uk/repository/entry/show/Polar+Data+Centre/DOI/BEDMAP2+-+Ice+thickness%2C+bed+and+surface+elevation+for+Antarctica+-+gridding+products/bedmap2_tiff?entryid=synth%3Afa5d606c-dc95-47ee-9016-7a82e446f2f2%3AL2JlZG1hcDJfdGlmZg%3D%3D&output=zip.zipgroup" fname = [p for p in path if p.endswith(f"{layer}.tif")][0] + grid = xr.load_dataarray(fname).squeeze() + # # calculate icebase as surface-thickness + # if layer == "icebase": + # fname = [p for p in path if p.endswith("surface.tif")][0] + # grid = xr.load_dataarray(fname).squeeze() + # surface = resample_grid( + # grid, + # initial_spacing=initial_spacing, + # initial_region=initial_region, + # initial_registration=initial_registration, + # spacing=spacing, + # region=region, + # registration=registration, + # **kwargs, + # ) - # calculate icebase as surface-thickness - if layer == "icebase": - fname = [p for p in path if p.endswith("surface.tif")][0] - grid = xr.load_dataarray(fname).squeeze() - surface = resample_grid( - grid, - initial_spacing=initial_spacing, - initial_region=initial_region, - initial_registration=initial_registration, - spacing=spacing, - region=region, - registration=registration, - **kwargs, - ) - - fname = [p for p in path if p.endswith("thickness.tif")][0] - grid = xr.load_dataarray(fname).squeeze() - thickness = resample_grid( - grid, - initial_spacing=initial_spacing, - initial_region=initial_region, - initial_registration=initial_registration, - spacing=spacing, - region=region, - registration=registration, - **kwargs, - ) - - # this changes the registration from pixel to gridline - resampled = surface - thickness - - elif layer in [ - "bed", - "coverage", - "grounded_bed_uncertainty", - "icemask_grounded_and_shelves", - "lakemask_vostok", - "rockmask", - "surface", - "thickness", - "thickness_uncertainty_5km", - "gl04c_geiod_to_WGS84", - ]: - - fname = [p for p in path if p.endswith(f"{layer}.tif")][0] - grid = xr.load_dataarray(fname).squeeze() - resampled = resample_grid( - grid, - initial_spacing=initial_spacing, - initial_region=initial_region, - initial_registration=initial_registration, - spacing=spacing, - region=region, - registration=registration, - **kwargs, - ) + # fname = [p for p in path if p.endswith("thickness.tif")][0] + # grid = xr.load_dataarray(fname).squeeze() + # thickness = resample_grid( + # grid, + # initial_spacing=initial_spacing, + # initial_region=initial_region, + # initial_registration=initial_registration, + # spacing=spacing, + # region=region, + # registration=registration, + # **kwargs, + # ) - else: - raise ValueError("invalid layer string") + # # this changes the registration from pixel to gridline + # resampled = surface - thickness + + # elif layer in [ + # "bed", + # "coverage", + # "grounded_bed_uncertainty", + # "icemask_grounded_and_shelves", + # "lakemask_vostok", + # "rockmask", + # "surface", + # "thickness", + # "thickness_uncertainty_5km", + # "gl04c_geiod_to_WGS84", + # ]: + + # fname = [p for p in path if p.endswith(f"{layer}.tif")][0] + # grid = xr.load_dataarray(fname).squeeze() + # resampled = resample_grid( + # grid, + # initial_spacing=initial_spacing, + # initial_region=initial_region, + # initial_registration=initial_registration, + # spacing=spacing, + # region=region, + # registration=registration, + # **kwargs, + # ) - # change layer elevation to be relative to the ellipsoid instead of the geoid - if reference == "ellipsoid" and layer in [ - "surface", - "icebase", - "bed", - ]: - geoid_file = [p for p in path if p.endswith("gl04c_geiod_to_WGS84.tif")][0] - geoid = xr.load_dataarray(geoid_file).squeeze() - resampled_geoid = resample_grid( - geoid, - initial_spacing=initial_spacing, - initial_region=initial_region, - initial_registration=initial_registration, - spacing=spacing, - region=region, - registration=registration, - **kwargs, - ) + # else: + # raise ValueError("invalid layer string") + + # # change layer elevation to be relative to the ellipsoid instead of the geoid + # if reference == "ellipsoid" and layer in [ + # "surface", + # "icebase", + # "bed", + # ]: + # geoid_file = [p for p in path if p.endswith("gl04c_geiod_to_WGS84.tif")][0] + # geoid = xr.load_dataarray(geoid_file).squeeze() + # resampled_geoid = resample_grid( + # geoid, + # initial_spacing=initial_spacing, + # initial_region=initial_region, + # initial_registration=initial_registration, + # spacing=spacing, + # region=region, + # registration=registration, + # **kwargs, + # ) - final_grid = resampled + resampled_geoid - # geoid.close() - # resampled_geoid.close() + # final_grid = resampled + resampled_geoid - elif reference not in ["ellipsoid", "geoid"]: - raise ValueError("invalid reference string") + # elif reference not in ["ellipsoid", "geoid"]: + # raise ValueError("invalid reference string") - else: - final_grid = resampled + # else: + # final_grid = resampled - # replace nans with 0's - if fill_nans is True: - if layer in ["surface", "thickness", "icebase"]: - # pygmt.grdfill(final_grid, mode='c0') # doesn't work, maybe grid is too big - # this changes the registration from pixel to gridline - final_grid = final_grid.fillna(0) + # # replace nans with 0's + # if fill_nans is True: + # if layer in ["surface", "thickness", "icebase"]: + # # pygmt.grdfill(final_grid, mode='c0') # doesn't work, maybe grid is too big + # # this changes the registration from pixel to gridline + # final_grid = final_grid.fillna(0) - if plot is True: - final_grid.plot(robust=True) - if info is True: - print(pygmt.grdinfo(final_grid)) + # if plot is True: + # final_grid.plot(robust=True) + # if info is True: + # print(pygmt.grdinfo(final_grid)) - return final_grid + return grid #final_grid def REMA( From 23f00650322da68f8e5584f634f8a111ca188f31 Mon Sep 17 00:00:00 2001 From: mdtanker Date: Fri, 4 Nov 2022 13:35:06 +1300 Subject: [PATCH 2/3] bedmap2 preprocesses as Zarr files --- antarctic_plots/fetch.py | 296 ++++++++++++++++------------ antarctic_plots/tests/test_fetch.py | 79 +++++--- 2 files changed, 222 insertions(+), 153 deletions(-) diff --git a/antarctic_plots/fetch.py b/antarctic_plots/fetch.py index 953d72aa..5b72ceb2 100644 --- a/antarctic_plots/fetch.py +++ b/antarctic_plots/fetch.py @@ -1323,20 +1323,32 @@ def bedmap_points( def bedmap2( layer: str, - # reference: str = "geoid", - # plot: bool = False, - # info: bool = False, - # region=None, - # spacing=None, - # registration=None, - # fill_nans=False, - # **kwargs, + reference: str = "geoid", + plot: bool = False, + info: bool = False, + region=None, + spacing=None, + registration=None, + fill_nans=False, + **kwargs, ) -> xr.DataArray: """ - Load bedmap2 data. All grids are by default referenced to the g104c geoid. Use the + Load bedmap2 data as xarray.DataArrays + from Fretwell et a. 2022: BEDMAP2 - Ice thickness, bed and surface elevation for + Antarctica - gridding products (Version 1.0) [Data set]. NERC EDS UK Polar Data + Centre. + DOI: https://doi.org/10.5285/FA5D606C-DC95-47EE-9016-7A82E446F2F2 + accessed from https://ramadda.data.bas.ac.uk/repository/entry/show?entryid=fa5d606c-dc95-47ee-9016-7a82e446f2f2. # noqa + + All grids are by default referenced to the g104c geoid. Use the 'reference' parameter to convert to the ellipsoid. - Note, nan's in surface grid are set to 0. - from https://doi.org/10.5194/tc-7-375-2013. + + Unlike Bedmachine data, Bedmap2 surface and icethickness contain NaN's over the + ocean, instead of 0's. To fill these NaN's with 0's, set `fill_nans=True`. + Note, this only makes since if the reference is the geoid, therefore, if + `reference='ellipsoid` and `fill_nans=True`, the nan's will be filled before + converting the results to the geoid (just for surface, since thickness isn't + relative to anything). Parameters ---------- @@ -1344,7 +1356,7 @@ def bedmap2( choose which layer to fetch: "bed", "coverage", "grounded_bed_uncertainty", "icemask_grounded_and_shelves", "lakemask_vostok", "rockmask", "surface", "thickness", - "thickness_uncertainty_5km", "gl04c_geiod_to_WGS84", "icebase" + "thickness_uncertainty_5km", "icebase" reference : str choose whether heights are referenced to 'geoid' (g104c) or 'ellipsoid' (WGS84), by default is 'geoid' @@ -1362,6 +1374,15 @@ def bedmap2( xr.DataArray Returns a loaded, and optional clip/resampled grid of Bedmap2. """ + + # download url + url = ( + "https://ramadda.data.bas.ac.uk/repository/entry/show/Polar+Data+Centre/" + "DOI/BEDMAP2+-+Ice+thickness%2C+bed+and+surface+elevation+for+Antarctica" + "+-+gridding+products/bedmap2_tiff?entryid=synth%3Afa5d606c-dc95-47ee-9016" + "-7a82e446f2f2%3AL2JlZG1hcDJfdGlmZg%3D%3D&output=zip.zipgroup" + ) + # Declare initial grid values, of .nc files not .tiff files # use utils.get_grid_info(xr.load_dataset(file).band_data # several of the layers have different values @@ -1371,133 +1392,162 @@ def bedmap2( initial_registration = "g" elif layer == "thickness_uncertainty_5km": - initial_region = [-3401500, 3403500, -3397500, 3397500] + initial_region = [-3399000.0, 3401000.0, -3400000.0, 3400000.0] initial_spacing = 5e3 - initial_registration = "p" + initial_registration = "g" - else: + elif layer in [ + "bed", + "coverage", + "grounded_bed_uncertainty", + "icemask_grounded_and_shelves", + "rockmask", + "surface", + "thickness", + "icebase", + ]: initial_region = [-3333000, 3333000, -3333000, 3333000] initial_spacing = 1e3 initial_registration = "g" - # if region is None: - # region = initial_region - # if spacing is None: - # spacing = initial_spacing - # if registration is None: - # registration = initial_registration + else: + raise ValueError("invalid layer string") - # retrieve the specified layer file - path = pooch.retrieve( - url="https://secure.antarctica.ac.uk/data/bedmap2/bedmap2_tiff.zip", - fname="bedmap2_tiff.zip", - path=f"{pooch.os_cache('pooch')}/antarctic_plots/topography", - known_hash=None, - processor=pooch.Unzip(), - progressbar=True, - ) -# "https://ramadda.data.bas.ac.uk/repository/entry/show/Polar+Data+Centre/DOI/BEDMAP2+-+Ice+thickness%2C+bed+and+surface+elevation+for+Antarctica+-+gridding+products/bedmap2_tiff?entryid=synth%3Afa5d606c-dc95-47ee-9016-7a82e446f2f2%3AL2JlZG1hcDJfdGlmZg%3D%3D&output=zip.zipgroup" fname = [p for p in path if p.endswith(f"{layer}.tif")][0] - grid = xr.load_dataarray(fname).squeeze() - # # calculate icebase as surface-thickness - # if layer == "icebase": - # fname = [p for p in path if p.endswith("surface.tif")][0] - # grid = xr.load_dataarray(fname).squeeze() - # surface = resample_grid( - # grid, - # initial_spacing=initial_spacing, - # initial_region=initial_region, - # initial_registration=initial_registration, - # spacing=spacing, - # region=region, - # registration=registration, - # **kwargs, - # ) + if region is None: + region = initial_region + if spacing is None: + spacing = initial_spacing + if registration is None: + registration = initial_registration - # fname = [p for p in path if p.endswith("thickness.tif")][0] - # grid = xr.load_dataarray(fname).squeeze() - # thickness = resample_grid( - # grid, - # initial_spacing=initial_spacing, - # initial_region=initial_region, - # initial_registration=initial_registration, - # spacing=spacing, - # region=region, - # registration=registration, - # **kwargs, - # ) + def preprocessing(fname, action, pooch2): + "Unzip the folder, convert the tiffs to compressed .zarr files" + # extract each layer to it's own folder + fname = pooch.Unzip( + extract_dir=f"bedmap2_{layer}", + members=[f"bedmap2_tiff/bedmap2_{layer}.tif"], + )(fname, action, pooch2)[0] + # get the path to the layer's tif file + fname = Path(fname) - # # this changes the registration from pixel to gridline - # resampled = surface - thickness - - # elif layer in [ - # "bed", - # "coverage", - # "grounded_bed_uncertainty", - # "icemask_grounded_and_shelves", - # "lakemask_vostok", - # "rockmask", - # "surface", - # "thickness", - # "thickness_uncertainty_5km", - # "gl04c_geiod_to_WGS84", - # ]: - - # fname = [p for p in path if p.endswith(f"{layer}.tif")][0] - # grid = xr.load_dataarray(fname).squeeze() - # resampled = resample_grid( - # grid, - # initial_spacing=initial_spacing, - # initial_region=initial_region, - # initial_registration=initial_registration, - # spacing=spacing, - # region=region, - # registration=registration, - # **kwargs, - # ) + # Rename to the file to ***.zarr + fname_processed = fname.with_suffix(".zarr") - # else: - # raise ValueError("invalid layer string") - - # # change layer elevation to be relative to the ellipsoid instead of the geoid - # if reference == "ellipsoid" and layer in [ - # "surface", - # "icebase", - # "bed", - # ]: - # geoid_file = [p for p in path if p.endswith("gl04c_geiod_to_WGS84.tif")][0] - # geoid = xr.load_dataarray(geoid_file).squeeze() - # resampled_geoid = resample_grid( - # geoid, - # initial_spacing=initial_spacing, - # initial_region=initial_region, - # initial_registration=initial_registration, - # spacing=spacing, - # region=region, - # registration=registration, - # **kwargs, - # ) + # Only recalculate if new download or the processed file doesn't exist yet + if action in ("download", "update") or not fname_processed.exists(): + # load data + grid = xr.load_dataarray(fname).squeeze().drop_vars(["band", "spatial_ref"]) + grid = grid.to_dataset(name=layer) - # final_grid = resampled + resampled_geoid + # Save to disk + grid.to_zarr(fname_processed) - # elif reference not in ["ellipsoid", "geoid"]: - # raise ValueError("invalid reference string") + return str(fname_processed) - # else: - # final_grid = resampled + # calculate icebase as surface-thickness + if layer == "icebase": + # set layer variable so pooch retrieves correct file + layer = "surface" + fname = pooch.retrieve( + url=url, + fname="bedmap2_tiff.zip", + path=f"{pooch.os_cache('pooch')}/antarctic_plots/topography", + known_hash=None, + processor=preprocessing, + progressbar=True, + ) + # load zarr as a dataarray + surface = xr.open_zarr(fname)[layer] + + layer = "thickness" + # set layer variable so pooch retrieves correct file + fname = pooch.retrieve( + url=url, + fname="bedmap2_tiff.zip", + path=f"{pooch.os_cache('pooch')}/antarctic_plots/topography", + known_hash=None, + processor=preprocessing, + progressbar=True, + ) + # load zarr as a dataarray + thickness = xr.open_zarr(fname)[layer] + + # calculate icebase with the resampled surface and thickness + grid = surface - thickness + + # reset layer variable + layer = "icebase" + + elif layer in [ + "bed", + "coverage", + "grounded_bed_uncertainty", + "icemask_grounded_and_shelves", + "lakemask_vostok", + "rockmask", + "surface", + "thickness", + "thickness_uncertainty_5km", + ]: + # download/unzip all files, retrieve the specified layer file and convert to + # .zarr + fname = pooch.retrieve( + url=url, + fname="bedmap2_tiff.zip", + path=f"{pooch.os_cache('pooch')}/antarctic_plots/topography", + known_hash=None, + processor=preprocessing, + progressbar=True, + ) + # load zarr as a dataarray + grid = xr.open_zarr(fname)[layer] - # # replace nans with 0's - # if fill_nans is True: - # if layer in ["surface", "thickness", "icebase"]: - # # pygmt.grdfill(final_grid, mode='c0') # doesn't work, maybe grid is too big - # # this changes the registration from pixel to gridline - # final_grid = final_grid.fillna(0) + else: + raise ValueError("invalid layer string") + + # replace nans with 0's in surface, thickness or icebase grids + if fill_nans is True: + if layer in ["surface", "thickness", "icebase"]: + # pygmt.grdfill(final_grid, mode='c0') # doesn't work, maybe grid is too big + # this changes the registration from pixel to gridline + grid = grid.fillna(0) + + # change layer elevation to be relative to the ellipsoid instead of the geoid + if reference == "ellipsoid" and layer in [ + "surface", + "icebase", + "bed", + ]: + # get a grid of geoid values matching the user's input + geoid_correction = geoid( + spacing=initial_spacing, + region=initial_region, + registration=initial_registration, + ) + # convert grid to be referenced to the ellipsoid + grid = grid + geoid_correction + + elif reference not in ["ellipsoid", "geoid"]: + raise ValueError("invalid reference string") + + # resample grid to users input + resampled = resample_grid( + grid, + initial_spacing=initial_spacing, + initial_region=initial_region, + initial_registration=initial_registration, + spacing=spacing, + region=region, + registration=registration, + **kwargs, + ) - # if plot is True: - # final_grid.plot(robust=True) - # if info is True: - # print(pygmt.grdinfo(final_grid)) + if plot is True: + resampled.plot(robust=True) + if info is True: + print(pygmt.grdinfo(resampled)) - return grid #final_grid + return resampled def REMA( diff --git a/antarctic_plots/tests/test_fetch.py b/antarctic_plots/tests/test_fetch.py index d5004698..001f2339 100644 --- a/antarctic_plots/tests/test_fetch.py +++ b/antarctic_plots/tests/test_fetch.py @@ -18,6 +18,7 @@ def test_(): import os import geopandas as gpd +import numpy as np import pandas as pd import pytest from geopandas.testing import assert_geodataframe_equal @@ -529,59 +530,77 @@ def test_bedmachine_reference(): # utils.get_grid_info(grid) # %% bedmap2 -# test for all layers, but only test reference models with 1 layer +# test for all layers, but only test reference models with 1 layer and fill_nans with 1 +# layer test = [ ( - "icebase", - ("1000", [-3333000.0, 3333000.0, -3333000.0, 3333000.0], -2736.0, 3972.0, "g"), + dict(layer="bed"), + ("1000", [-3333000.0, 3333000.0, -3333000.0, 3333000.0], -7054.0, 3972.0, "g"), ), ( - "surface", - ("1000", [-3333500.0, 3333500.0, -3332500.0, 3332500.0], 1.0, 4082.0, "p"), + dict(layer="coverage"), + ("1000", [-3333000.0, 3333000.0, -3333000.0, 3333000.0], 1.0, 1.0, "g"), ), ( - "thickness", - ("1000", [-3333500.0, 3333500.0, -3332500.0, 3332500.0], 0.0, 4621.0, "p"), + dict(layer="grounded_bed_uncertainty"), + ("1000", [-3333000.0, 3333000.0, -3333000.0, 3333000.0], 0.0, 65535.0, "g"), ), ( - "bed", - ("1000", [-3333500.0, 3333500.0, -3332500.0, 3332500.0], -7054.0, 3972.0, "p"), + dict(layer="icemask_grounded_and_shelves"), + ("1000", [-3333000.0, 3333000.0, -3333000.0, 3333000.0], 0.0, 1.0, "g"), + ), + ( + dict(layer="rockmask"), + ("1000", [-3333000.0, 3333000.0, -3333000.0, 3333000.0], 0.0, 0.0, "g"), + ), + ( + dict(layer="surface"), + ("1000", [-3333000.0, 3333000.0, -3333000.0, 3333000.0], 1.0, 4082.0, "g"), + ), + ( + dict(layer="thickness"), + ("1000", [-3333000.0, 3333000.0, -3333000.0, 3333000.0], 0.0, 4621.0, "g"), + ), + ( + dict(layer="icebase"), + ("1000", [-3333000.0, 3333000.0, -3333000.0, 3333000.0], -2736.0, 3972.0, "g"), + ), + ( + dict(layer="lakemask_vostok"), + ("1000", [1190000.0, 1470000.0, -402000.0, -291000.0], 1.0, 1.0, "g"), ), ( - "gl04c_geiod_to_WGS84", + dict(layer="thickness_uncertainty_5km"), + ("5000", [-3399000.0, 3401000.0, -3400000.0, 3400000.0], 0.0, 65535.0, "g"), + ), + ( + dict(layer="icebase", reference="ellipsoid"), ( "1000", - [-3333500.0, 3333500.0, -3332500.0, 3332500.0], - -65.8680496216, - 36.6361198425, - "p", + [-3330000.0, 3330000.0, -3330000.0, 3330000.0], + -2784.4753418, + 3951.77758789, + "g", ), ), ] -@pytest.mark.working -def test_bedmap2_reference(): - grid = fetch.bedmap2(layer="surface", reference="ellipsoid") - expected = ( - "1000", - [-3333000.0, 3333000.0, -3333000.0, 3333000.0], - -50.4912605286, - 4090.53417969, - "g", - ) - assert utils.get_grid_info(grid) == pytest.approx(expected, rel=0.1) - - -@pytest.mark.working @pytest.mark.parametrize("test_input,expected", test) def test_bedmap2(test_input, expected): - grid = fetch.bedmap2(test_input) + grid = fetch.bedmap2(**test_input) assert utils.get_grid_info(grid) == pytest.approx(expected, rel=0.1) -# grid = fetch.bedmap2(layer="surface", reference="geoid") +def test_bedmap2_fill_nans(): + grid = fetch.bedmap2(layer="surface") + filled_grid = fetch.bedmap2(layer="surface", fill_nans=True) + assert np.nanmean(grid) == pytest.approx(1964.5349, rel=0.1) + assert np.nanmean(filled_grid) == pytest.approx(602.32306, rel=0.1) + + +# grid = fetch.bedmap2(layer="surface") # utils.get_grid_info(grid) # %% Bedmap points From a275accdecfacd15cb0df72ffe7c6ce4318c26d1 Mon Sep 17 00:00:00 2001 From: mdtanker Date: Mon, 7 Nov 2022 07:30:03 +1300 Subject: [PATCH 3/3] fixed geoid conversion --- antarctic_plots/fetch.py | 92 ++++++++++++++++++++--------- antarctic_plots/tests/test_fetch.py | 48 +++++++++++---- 2 files changed, 101 insertions(+), 39 deletions(-) diff --git a/antarctic_plots/fetch.py b/antarctic_plots/fetch.py index 5b72ceb2..5605cded 100644 --- a/antarctic_plots/fetch.py +++ b/antarctic_plots/fetch.py @@ -1323,7 +1323,7 @@ def bedmap_points( def bedmap2( layer: str, - reference: str = "geoid", + reference: str = "gl04c", plot: bool = False, info: bool = False, region=None, @@ -1340,8 +1340,10 @@ def bedmap2( DOI: https://doi.org/10.5285/FA5D606C-DC95-47EE-9016-7A82E446F2F2 accessed from https://ramadda.data.bas.ac.uk/repository/entry/show?entryid=fa5d606c-dc95-47ee-9016-7a82e446f2f2. # noqa - All grids are by default referenced to the g104c geoid. Use the - 'reference' parameter to convert to the ellipsoid. + + All grids are by default referenced to the gl04c geoid. Use the + reference='ellipsoid' to convert to the WGS-84 ellipsoid or reference='eigen' to + convert to the EIGEN-6c4 geoid. Unlike Bedmachine data, Bedmap2 surface and icethickness contain NaN's over the ocean, instead of 0's. To fill these NaN's with 0's, set `fill_nans=True`. @@ -1356,10 +1358,10 @@ def bedmap2( choose which layer to fetch: "bed", "coverage", "grounded_bed_uncertainty", "icemask_grounded_and_shelves", "lakemask_vostok", "rockmask", "surface", "thickness", - "thickness_uncertainty_5km", "icebase" + "thickness_uncertainty_5km", "gl04c_geiod_to_WGS84", "icebase" reference : str - choose whether heights are referenced to 'geoid' (g104c) or 'ellipsoid' - (WGS84), by default is 'geoid' + choose whether heights are referenced to the EIGEN-6c4 geoid 'eigen', the WGS84 + ellipsoid, 'ellipsoid', or by default the 'gl04c' geoid. plot : bool, optional choose to plot grid, by default False info : bool, optional @@ -1368,6 +1370,13 @@ def bedmap2( GMT-format region to clip the loaded grid to, by default doesn't clip spacing : str or int, optional grid spacing to resample the loaded grid to, by default 10e3 + registration : str, optional, + choose between 'g' (gridline) or 'p' (pixel) registration types, by default is + the original type of the grid + fill_nans : bool, optional, + choose whether to fill nans in 'surface' and 'thickness' with 0. If converting + to reference to the geoid, will fill nan's before conversion, by default is + False Returns ------- @@ -1404,6 +1413,7 @@ def bedmap2( "rockmask", "surface", "thickness", + "gl04c_geiod_to_WGS84", "icebase", ]: initial_region = [-3333000, 3333000, -3333000, 3333000] @@ -1423,10 +1433,13 @@ def bedmap2( def preprocessing(fname, action, pooch2): "Unzip the folder, convert the tiffs to compressed .zarr files" # extract each layer to it's own folder - fname = pooch.Unzip( - extract_dir=f"bedmap2_{layer}", - members=[f"bedmap2_tiff/bedmap2_{layer}.tif"], - )(fname, action, pooch2)[0] + if layer == "gl04c_geiod_to_WGS84": + member = ["bedmap2_tiff/gl04c_geiod_to_WGS84.tif"] + else: + member = [f"bedmap2_tiff/bedmap2_{layer}.tif"] + fname = pooch.Unzip(extract_dir=f"bedmap2_{layer}", members=member,)( + fname, action, pooch2 + )[0] # get the path to the layer's tif file fname = Path(fname) @@ -1488,6 +1501,7 @@ def preprocessing(fname, action, pooch2): "surface", "thickness", "thickness_uncertainty_5km", + "gl04c_geiod_to_WGS84", ]: # download/unzip all files, retrieve the specified layer file and convert to # .zarr @@ -1505,30 +1519,45 @@ def preprocessing(fname, action, pooch2): else: raise ValueError("invalid layer string") - # replace nans with 0's in surface, thickness or icebase grids + # replace nans with 0's in surface or thickness grids if fill_nans is True: - if layer in ["surface", "thickness", "icebase"]: + if layer in ["surface", "thickness"]: # pygmt.grdfill(final_grid, mode='c0') # doesn't work, maybe grid is too big # this changes the registration from pixel to gridline grid = grid.fillna(0) - # change layer elevation to be relative to the ellipsoid instead of the geoid - if reference == "ellipsoid" and layer in [ - "surface", - "icebase", - "bed", - ]: - # get a grid of geoid values matching the user's input - geoid_correction = geoid( - spacing=initial_spacing, - region=initial_region, - registration=initial_registration, + # change layer elevation to be relative to different reference frames. + if layer in ["surface", "icebase", "bed"]: + # set layer variable so pooch retrieves the geoid convertion file + layer = "gl04c_geiod_to_WGS84" + fname = pooch.retrieve( + url=url, + fname="bedmap2_tiff.zip", + path=f"{pooch.os_cache('pooch')}/antarctic_plots/topography", + known_hash=None, + processor=preprocessing, + progressbar=True, ) - # convert grid to be referenced to the ellipsoid - grid = grid + geoid_correction - - elif reference not in ["ellipsoid", "geoid"]: - raise ValueError("invalid reference string") + # load zarr as a dataarray + geoid_2_ellipsoid = xr.open_zarr(fname)[layer] + if reference == "ellipsoid": + # convert to the ellipsoid + grid = grid + geoid_2_ellipsoid + elif reference == "eigen": + # convert to the ellipsoid + grid = grid + geoid_2_ellipsoid + # get a grid of EIGEN geoid values matching the user's input + eigen_correction = geoid( + spacing=initial_spacing, + region=initial_region, + registration=initial_registration, + ) + # convert from ellipsoid to eigen geoid + grid = grid - eigen_correction + elif reference == "gl04c": + pass + else: + raise ValueError("invalid reference string") # resample grid to users input resampled = resample_grid( @@ -2084,8 +2113,13 @@ def geoid( **kwargs, ) -> xr.DataArray: """ - Loads a grid of Antarctic geoid height derived from the EIGEN-6C4 spherical + Loads a grid of Antarctic geoid heights derived from the EIGEN-6C4 spherical harmonic model of Earth's gravity field. Originally at 10 arc-min resolution. + Negative values indicate the geoid is below the ellipsoid surface and vice-versa. + To convert a topographic grid which is referenced to the ellipsoid to be referenced + to the geoid, add this grid. + To convert a topographic grid which is referenced to the geoid to be reference to the + ellipsoid, subtract this grid. orignally from https://dataservices.gfz-potsdam.de/icgem/showshort.php?id=escidoc:1119897 # noqa Accessed via the Fatiando data repository https://github.com/fatiando-data/earth-geoid-10arcmin # noqa diff --git a/antarctic_plots/tests/test_fetch.py b/antarctic_plots/tests/test_fetch.py index 001f2339..770d3c2f 100644 --- a/antarctic_plots/tests/test_fetch.py +++ b/antarctic_plots/tests/test_fetch.py @@ -574,16 +574,6 @@ def test_bedmachine_reference(): dict(layer="thickness_uncertainty_5km"), ("5000", [-3399000.0, 3401000.0, -3400000.0, 3400000.0], 0.0, 65535.0, "g"), ), - ( - dict(layer="icebase", reference="ellipsoid"), - ( - "1000", - [-3330000.0, 3330000.0, -3330000.0, 3330000.0], - -2784.4753418, - 3951.77758789, - "g", - ), - ), ] @@ -593,6 +583,44 @@ def test_bedmap2(test_input, expected): assert utils.get_grid_info(grid) == pytest.approx(expected, rel=0.1) +def test_bedmap2_reference(): + # fetch variations of grids and reference models + region = [-101e3, -100e3, -51e3, -50e3] + gl04c_grid = fetch.bedmap2( + layer="gl04c_geiod_to_WGS84", + region=region, + ) + eigen_grid = fetch.geoid( + region=region, + spacing=1e3, + ) + surface_eigen_grid = fetch.bedmap2( + layer="surface", + reference="eigen", + region=region, + ) + surface_ellipsoid_grid = fetch.bedmap2( + layer="surface", + reference="ellipsoid", + region=region, + ) + surface_gl04c_grid = fetch.bedmap2( + layer="surface", + reference="gl04c", + region=region, + ) + # get mean values + gl04c = np.nanmean(gl04c_grid) + eigen = np.nanmean(eigen_grid) + surface_eigen = np.nanmean(surface_eigen_grid) + surface_ellipsoid = np.nanmean(surface_ellipsoid_grid) + surface_gl04c = np.nanmean(surface_gl04c_grid) + assert surface_ellipsoid - eigen == pytest.approx(surface_eigen, rel=0.1) + assert surface_ellipsoid - gl04c == pytest.approx(surface_gl04c, rel=0.1) + assert surface_gl04c + gl04c == pytest.approx(surface_ellipsoid, rel=0.1) + assert surface_eigen + eigen == pytest.approx(surface_ellipsoid, rel=0.1) + + def test_bedmap2_fill_nans(): grid = fetch.bedmap2(layer="surface") filled_grid = fetch.bedmap2(layer="surface", fill_nans=True)