diff --git a/build_Bamber_greenland.py b/build_Bamber_greenland.py index b13eaf3..8db825b 100755 --- a/build_Bamber_greenland.py +++ b/build_Bamber_greenland.py @@ -18,6 +18,8 @@ from data import insar from data import icebridge from data import ice2sea +from data import csatho + """ Build a CISM dataset @@ -38,6 +40,10 @@ # lc_massCon = "data/IceBridge/Greenland/MCdataset-2014-11-19.nc" lc_mask = "data/Ice2Sea/ice2sea_Greenland_geometry_icesheet_mask_Zurich.nc" +lc_csatho = ( + "data/Csatho2014/GreenlandIceSheetdhdt_csatho/" + "greenland_ice_sheet_dhdt_icesat_atm_l1a_to_l2f.cdf" +) # ==== SETUP ===== # get args, time @@ -116,6 +122,9 @@ speak.verbose(args, " Found Zurich mask") nc_mask.close() +nc_csatho = get_nc_file(lc_csatho, "r") +speak.verbose(args, " Found Csatho data") + speak.verbose(args, "\n All data files found!") @@ -202,6 +211,16 @@ nc_bamber.close() nc_massCon.close() +# == Csatho dHdt data needs to follow thickness interpolation) == +# Approx. 8km dataset +# ===================== +speak.notquiet(args, "\nGetting dhdt from the Csatho data.") +csatho.dhdt_bamber( + args, nc_csatho, nc_base, base, proj_epsg3413, proj_eigen_gl04c +) +speak.notquiet(args, " Done!") +nc_csatho.close() + # ==== Zurich mask ===== # apply mask, and get # new surface variable diff --git a/data/csatho.py b/data/csatho.py new file mode 100644 index 0000000..a61ddab --- /dev/null +++ b/data/csatho.py @@ -0,0 +1,194 @@ +""" +data.csatho : Csatho 2014 data import module. + +This module provides functions to import surface elevation change (rates) data +from Csatho 2014 into a CISM dataset. + +Functions list: + * + +Notes +----- +This data was downloaded from the citation below. + +The data uses the EPSG:32624 (WGS84 / UTM zone 24N) projection: + * Tranverse_Mercator + * WGS84 ellipsoid + * Latitude of projection origin = 0 degrees + * Central meridian = -39 degrees + * false eastings = 500000 + * false northings = 0 + * scale factor: 0.9996 + * Proj4: +proj=utm +zone=24 +datum=WGS84 +units=m +no_defs + +Because the original data files are written in terms of ICESat Mission phases, +the data have been processed into a more obvious form. The processed data has +been distributed into analysis year NetCDF files, where the extension on the +.nc files indicate the years over which the rates of change were calculated +(e.g. "0304" were calculated between 2003 and 2004). Because the coordinates +appear to be transposed from standard convention in the original data, the +processed data's coordinates have been transposed (common features were +spot checked) and the thickness change variable has been renamed to dhdt. + +References +---------- +Beata M. Csatho, Anton F. Schenk, Cornelis J. van der Veen, Gregory Babonis, +Kyle Duncan, Soroush Rezvanbehbahani, Michiel R. van den Broeke, Sebastian B. +Simonsen, Sudhagar Nagarajan, Jan H. van Angelen: Greenland ice dynamics from +laser altimetry, Proceedings of the National Academy of Sciences Dec 2014, +111 (52) 18478-18483; DOI: 10.1073/pnas.1411680112 + +Data: https://arcticdata.io/catalog/view/doi:10.5065/D6HM56KS +""" + +import numpy as np +import scipy +import scipy.interpolate + +from util import speak +from util import projections +from util.ncfunc import copy_atts_add_fill + + +# ICESat campaign mid-dates in form of: {name: MMDDYY} +ICESAT_CAMPAIGN_MIDDATES = { + "L1A": "030503", + "L1B": "032503", + "L2A": "100703", + "L2B": "030504", + "L2C": "060504", + "L3A": "101904", + "L3B": "030705", + "L3C": "060605", + "L3D": "110705", + "L3E": "031106", + "L3F": "060906", + "L3G": "111006", + "L3H": "032907", + "L3I": "101907", + "L3J": "030508", + "L2D": "120608", + "L2E": "031509", + "L2F": "101509", +} + +MISSING_VAL = np.float32(2.0e36) + + +def dhdt_epsg3413(args, nc_csatho, nc_base, base, proj_epsg3413): + in_var = "L1A_to_L2F_Avg_dhdt" + csatho = projections.DataGrid() + csatho.lat = nc_csatho.variables["lat"] + csatho.lon = nc_csatho.variables["lon"] + csatho.dhdt = nc_csatho.variables[in_var] + + csatho.lat_grid, csatho.lon_grid = scipy.meshgrid( + csatho.lat[:], csatho.lon[:], indexing="ij" + ) + csatho.x, csatho.y = proj_epsg3413( + csatho.lon_grid.ravel(), csatho.lat_grid.ravel() + ) + + speak.verbose(args, " Interpolating dhdt and writing to base.") + vy = csatho.y[~csatho.dhdt[:].mask.ravel()] + vx = csatho.x[~csatho.dhdt[:].mask.ravel()] + vz = csatho.dhdt[:].ravel()[~csatho.dhdt[:].mask.ravel()] + + dhdt = np.ma.masked_invalid( + scipy.interpolate.griddata( + np.column_stack((vy, vx)), + vz, + (base.y_grid, base.x_grid), + method="linear", + ) + ) + + dhdt.mask = dhdt.mask | np.isclose(base.thk[:], 0.0) + + base.dhdt = nc_base.createVariable("dhdt", "f4", ("y", "x",)) + base.dhdt[:] = dhdt[:] + # Copy the attributes, somehow it's missing the missing value setter + copy_atts_add_fill(nc_csatho.variables[in_var], base.dhdt, MISSING_VAL) + + base.dhdt.long_name = "average {}--{} surface elevation change rate".format( + ICESAT_CAMPAIGN_MIDDATES["L1A"][-2:], + ICESAT_CAMPAIGN_MIDDATES["L2F"][-2:], + ) + base.dhdt.standard_name = "tendency_of_land_ice_thickness" + base.dhdt.units = "m year-1" + base.dhdt.grid_mapping = "epsg_3413" + base.dhdt.coordinates = "lon lat" + base.dhdt.source = ( + "https://arcticdata.io/catalog/view/doi:10.5065/D6HM56KS" + ) + base.dhdt.reference = ( + "Beata M. Csatho, Anton F. Schenk, Cornelis J. van der " + "Veen, Gregory Babonis, Kyle Duncan, Soroush Rezvanbehbahani, " + "Michiel R. van den Broeke, Sebastian B. Simonsen, Sudhagar Nagarajan, " + "Jan H. van Angelen: Greenland ice dynamics from laser altimetry, " + "Proceedings of the National Academy of Sciences Dec 2014, 111 (52) " + "18478-18483; DOI: 10.1073/pnas.1411680112" + ) + + +def dhdt_bamber( + args, nc_csatho, nc_base, base, proj_epsg3413, proj_eigen_gl04c +): + in_var = "L1A_to_L2F_Avg_dhdt" + csatho = projections.DataGrid() + csatho.lat = nc_csatho.variables["lat"] + csatho.lon = nc_csatho.variables["lon"] + csatho.dhdt = nc_csatho.variables[in_var] + + csatho.lat_grid, csatho.lon_grid = np.meshgrid( + csatho.lat[:], csatho.lon[:], indexing="ij" + ) + + csatho.x, csatho.y = proj_eigen_gl04c( + csatho.lon_grid.ravel(), csatho.lat_grid.ravel() + ) + + speak.verbose(args, " Interpolating dhdt and writing to base.") + vy = csatho.y[~csatho.dhdt[:].mask.ravel()] + vx = csatho.x[~csatho.dhdt[:].mask.ravel()] + vz = csatho.dhdt[:].ravel()[~csatho.dhdt[:].mask.ravel()] + + dhdt = np.ma.masked_invalid( + scipy.interpolate.griddata( + np.column_stack((vy, vx)), + vz, + (base.y_grid, base.x_grid), + method="linear", + ) + ) + + dhdt.mask |= np.isclose(base.thk[:], 0.0) + + dhdt = np.ma.masked_greater(dhdt, 1e10) + dhdt = dhdt.filled(MISSING_VAL) + base.dhdt = nc_base.createVariable("dhdt", "f4", ("y", "x",)) + base.dhdt[:] = dhdt[:] + + # Copy the attributes, somehow it's missing the missing value setter + copy_atts_add_fill(nc_csatho.variables[in_var], base.dhdt, MISSING_VAL) + + # Now replace all the metadata except for the fill/missing value + base.dhdt.long_name = "average {}--{} surface elevation change rate".format( + ICESAT_CAMPAIGN_MIDDATES["L1A"][-2:], + ICESAT_CAMPAIGN_MIDDATES["L2F"][-2:], + ) + base.dhdt.standard_name = "tendency_of_land_ice_thickness" + base.dhdt.units = "m year-1" + base.dhdt.grid_mapping = "mcb" + base.dhdt.coordinates = "lon lat" + base.dhdt.source = ( + "https://arcticdata.io/catalog/view/doi:10.5065/D6HM56KS" + ) + base.dhdt.reference = ( + "Beata M. Csatho, Anton F. Schenk, Cornelis J. van der " + "Veen, Gregory Babonis, Kyle Duncan, Soroush Rezvanbehbahani, " + "Michiel R. van den Broeke, Sebastian B. Simonsen, Sudhagar Nagarajan, " + "Jan H. van Angelen: Greenland ice dynamics from laser altimetry, " + "Proceedings of the National Academy of Sciences Dec 2014, 111 (52) " + "18478-18483; DOI: 10.1073/pnas.1411680112" + )