Skip to content

Commit

Permalink
Bedmap2 restructure redo (#123)
Browse files Browse the repository at this point in the history 
* diagnosing bedmap2 issues

* bedmap2 preprocesses as Zarr files

* fixed geoid conversion
  • Loading branch information
@mdtanker
mdtanker authored Nov 6, 2022
1 parent 7cc8b09 commit b89040d
Show file tree
Hide file tree
Showing 2 changed files with 260 additions and 131 deletions.
278 changes: 180 additions & 98 deletions antarctic_plots/fetch.py
View file
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

0 comments on commit b89040d

Please sign in to comment.