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.

Sign up for GitHub

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

Jump to bottom

Bedmap2 restructure redo #123

Merged
merged 3 commits into from
Nov 6, 2022
Merged

Bedmap2 restructure redo #123

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
40d893c
diagnosing bedmap2 issues
mdtanker Nov 3, 2022
23f0065
bedmap2 preprocesses as Zarr files
mdtanker Nov 4, 2022
a275acc
fixed geoid conversion
mdtanker Nov 6, 2022
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
278 changes: 180 additions & 98 deletions antarctic_plots/fetch.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -1333,10 +1333,24 @@ def bedmap2(
**kwargs,
) -> xr.DataArray:
"""
Load bedmap2 data. 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.
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 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`.
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
----------
Expand All @@ -1346,8 +1360,8 @@ def bedmap2(
"lakemask_vostok", "rockmask", "surface", "thickness",
"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
Expand All @@ -1356,30 +1370,58 @@ 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
-------
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())

# 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
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]
initial_region = [-3399000.0, 3401000.0, -3400000.0, 3400000.0]
initial_spacing = 5e3
initial_registration = "p"
initial_registration = "g"

else:
# y lims are differnt if inputting fname straight to utils.get_grid_info()
initial_region = [-3333500, 3333500, -3332500, 3332500]
elif layer in [
"bed",
"coverage",
"grounded_bed_uncertainty",
"icemask_grounded_and_shelves",
"rockmask",
"surface",
"thickness",
"gl04c_geiod_to_WGS84",
"icebase",
]:
initial_region = [-3333000, 3333000, -3333000, 3333000]
initial_spacing = 1e3
initial_registration = "p"
initial_registration = "g"

else:
raise ValueError("invalid layer string")

if region is None:
region = initial_region
Expand All @@ -1388,46 +1430,66 @@ def bedmap2(
if registration is None:
registration = initial_registration

# 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,
)
def preprocessing(fname, action, pooch2):
"Unzip the folder, convert the tiffs to compressed .zarr files"
# extract each layer to it's own folder
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)

# Rename to the file to ***.zarr
fname_processed = fname.with_suffix(".zarr")

# 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)

# Save to disk
grid.to_zarr(fname_processed)

return str(fname_processed)

# 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,
# 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]

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,
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]

# this changes the registration from pixel to gridline
resampled = surface - thickness
# calculate icebase with the resampled surface and thickness
grid = surface - thickness

# reset layer variable
layer = "icebase"

elif layer in [
"bed",
Expand All @@ -1441,65 +1503,80 @@ def bedmap2(
"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,
# 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]

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()

elif reference not in ["ellipsoid", "geoid"]:
raise ValueError("invalid reference string")

else:
final_grid = resampled

# replace nans with 0's
# 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
final_grid = final_grid.fillna(0)
grid = grid.fillna(0)

# 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,
)
# 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(
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)
resampled.plot(robust=True)
if info is True:
print(pygmt.grdinfo(final_grid))
print(pygmt.grdinfo(resampled))

return final_grid
return resampled


def REMA(
Expand Down Expand Up @@ -2036,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
Expand Down
Loading