|
| 1 | +#!/usr/bin/env python |
| 2 | +import os |
| 3 | +import shutil |
| 4 | + |
| 5 | +import numpy as np |
| 6 | +import pyproj |
| 7 | +import xarray as xr |
| 8 | +from mpas_tools.cime.constants import constants as cime_constants |
| 9 | +from pyremap import Remapper, ProjectionGridDescriptor |
| 10 | +from pyremap.polar import get_antarctic_stereographic_projection |
| 11 | + |
| 12 | +from mpas_analysis.shared.io.download import download_files |
| 13 | + |
| 14 | + |
| 15 | +def download_paolo(out_filename): |
| 16 | + """ |
| 17 | + Remap the Paolo et al. (2023) melt rates at 1 km resolution to an MPAS |
| 18 | + mesh |
| 19 | +
|
| 20 | + Parameters |
| 21 | + ---------- |
| 22 | + out_filename : str |
| 23 | + The original Paolo et al. (2023) melt rates |
| 24 | + """ |
| 25 | + |
| 26 | + if os.path.exists(out_filename): |
| 27 | + return |
| 28 | + |
| 29 | + download_files(fileList=['ANT_G1920V01_IceShelfMelt.nc'], |
| 30 | + urlBase='https://its-live-data.s3.amazonaws.com/height_change/Antarctica/Floating', |
| 31 | + outDir='.') |
| 32 | + shutil.move('ANT_G1920V01_IceShelfMelt.nc', out_filename) |
| 33 | + |
| 34 | + |
| 35 | +def process_paolo(in_filename, out_filename): |
| 36 | + """ |
| 37 | + Convert Paolo et al. (2023) melt rates to freshwater equivalent |
| 38 | +
|
| 39 | + Parameters |
| 40 | + ---------- |
| 41 | + in_filename : str |
| 42 | + The original Paolo et al. (2023) melt rates |
| 43 | +
|
| 44 | + out_filename : str |
| 45 | + The Paolo et al. (2023) melt rates in freshwater equivalent |
| 46 | + """ |
| 47 | + if os.path.exists(out_filename): |
| 48 | + return |
| 49 | + |
| 50 | + print(f'Reading {in_filename}...') |
| 51 | + |
| 52 | + with xr.open_dataset(in_filename) as ds_in: |
| 53 | + |
| 54 | + x = ds_in.x |
| 55 | + y = ds_in.y |
| 56 | + melt_rate = ds_in.melt_mean |
| 57 | + melt_rate_uncertainty = np.sqrt((ds_in.melt_err**2).mean(dim='time')) |
| 58 | + |
| 59 | + print('done.') |
| 60 | + |
| 61 | + print('Creating an xarray dataset...') |
| 62 | + ds = xr.Dataset() |
| 63 | + |
| 64 | + projection = pyproj.Proj('+proj=stere +lat_ts=-71.0 +lat_0=-90 +lon_0=0.0 ' |
| 65 | + '+k_0=1.0 +x_0=0.0 +y_0=0.0 +ellps=WGS84') |
| 66 | + latlon_projection = pyproj.Proj(proj='latlong', datum='WGS84') |
| 67 | + |
| 68 | + print('Computing lat/lon...') |
| 69 | + x_2d, y_2d = np.meshgrid(x.values, y.values) |
| 70 | + transformer = pyproj.Transformer.from_proj(projection, |
| 71 | + latlon_projection) |
| 72 | + lon, lat = transformer.transform(x_2d, y_2d) |
| 73 | + print('done.') |
| 74 | + |
| 75 | + # Paolo et al. (2023) ice density |
| 76 | + rho_ice = 917. |
| 77 | + rho_fw = cime_constants['SHR_CONST_RHOFW'] |
| 78 | + ice_to_fw_equiv = rho_ice / rho_fw |
| 79 | + |
| 80 | + ds['x'] = x |
| 81 | + ds['y'] = y |
| 82 | + ds['lon'] = (('y', 'x'), lon) |
| 83 | + ds.lon.attrs['units'] = 'degrees' |
| 84 | + ds['lat'] = (('y', 'x'), lat) |
| 85 | + ds.lat.attrs['units'] = 'degrees' |
| 86 | + ds['meltRate'] = -ice_to_fw_equiv * melt_rate |
| 87 | + ds.meltRate.attrs['units'] = 'm/yr of freshwater' |
| 88 | + ds['meltRateUncertainty'] = ice_to_fw_equiv * melt_rate_uncertainty |
| 89 | + ds.meltRateUncertainty.attrs['units'] = 'm/yr of freshwater' |
| 90 | + print('Writing the dataset...') |
| 91 | + ds.to_netcdf(out_filename) |
| 92 | + print('done.') |
| 93 | + |
| 94 | + |
| 95 | +def remap_paolo(in_filename, out_prefix, date, task_count=128): |
| 96 | + """ |
| 97 | + Remap Paolo et al. (2023) melt rates to comparison grids |
| 98 | +
|
| 99 | + Parameters |
| 100 | + ---------- |
| 101 | + in_filename : str |
| 102 | + The Paolo et al. (2023) melt rates in NetCDF format |
| 103 | +
|
| 104 | + out_prefix : str |
| 105 | + A prefix for the file to contain the Paolo et al. (2023) melt |
| 106 | + rates and melt fluxes remapped to the comparison grid |
| 107 | +
|
| 108 | + date : str |
| 109 | + A date string to append to the file name. |
| 110 | +
|
| 111 | + task_count : int |
| 112 | + The number of MPI tasks to use to create the mapping file |
| 113 | + """ |
| 114 | + ds = xr.open_dataset(in_filename) |
| 115 | + |
| 116 | + melt_attrs = ds.meltRate.attrs |
| 117 | + uncert_attrs = ds.meltRateUncertainty.attrs |
| 118 | + |
| 119 | + mask = ds.meltRate.notnull() |
| 120 | + ds['meltRate'] = ds.meltRate.where(mask, 0.) |
| 121 | + ds['meltMask'] = mask.astype(float) |
| 122 | + mask = ds.meltRateUncertainty.notnull() |
| 123 | + ds['meltRateUncertSqr'] = (ds.meltRateUncertainty**2).where(mask, 0.) |
| 124 | + ds['uncertMask'] = mask.astype(float) |
| 125 | + ds = ds.drop_vars(['lat', 'lon', 'meltRateUncertainty']) |
| 126 | + |
| 127 | + in_x = ds.x.values |
| 128 | + in_y = ds.y.values |
| 129 | + lx = np.abs(1e-3 * (in_x[-1] - in_x[0])) |
| 130 | + ly = np.abs(1e-3 * (in_y[-1] - in_y[0])) |
| 131 | + dx = np.abs(1e-3 * (in_x[1] - in_x[0])) |
| 132 | + |
| 133 | + in_grid_name = f'{lx:g}x{ly:g}km_{dx:g}km_Antarctic_stereo' |
| 134 | + |
| 135 | + in_projection = pyproj.Proj('+proj=stere +lat_ts=-71.0 +lat_0=-90 ' |
| 136 | + '+lon_0=0.0 +k_0=1.0 +x_0=0.0 +y_0=0.0 ' |
| 137 | + '+ellps=WGS84') |
| 138 | + |
| 139 | + in_descriptor = ProjectionGridDescriptor.create( |
| 140 | + in_projection, in_x, in_y, in_grid_name) |
| 141 | + |
| 142 | + width = 6000. |
| 143 | + reses = [1., 4., 10.] |
| 144 | + |
| 145 | + for res in reses: |
| 146 | + x_max = 0.5 * width * 1e3 |
| 147 | + nx = int(width / res) + 1 |
| 148 | + out_x = np.linspace(-x_max, x_max, nx) |
| 149 | + |
| 150 | + out_grid_name = f'{width:g}x{width:g}km_{res:g}km_Antarctic_stereo' |
| 151 | + |
| 152 | + out_projection = get_antarctic_stereographic_projection() |
| 153 | + |
| 154 | + out_descriptor = ProjectionGridDescriptor.create( |
| 155 | + out_projection, out_x, out_x, out_grid_name) |
| 156 | + |
| 157 | + method = 'conserve' |
| 158 | + |
| 159 | + map_filename = f'map_{in_grid_name}_to_{out_grid_name}_{method}.nc' |
| 160 | + |
| 161 | + remapper = Remapper(in_descriptor, out_descriptor, map_filename) |
| 162 | + |
| 163 | + if not os.path.exists(map_filename): |
| 164 | + remapper.build_mapping_file(method=method, mpiTasks=task_count, |
| 165 | + esmf_parallel_exec='srun') |
| 166 | + |
| 167 | + ds_out = remapper.remap(ds) |
| 168 | + mask = ds_out.meltMask > 0. |
| 169 | + ds_out['meltRate'] = ds_out.meltRate.where(mask) |
| 170 | + ds_out.meltRate.attrs = melt_attrs |
| 171 | + mask = ds_out.uncertMask > 0. |
| 172 | + ds_out['meltRateUncertainty'] = \ |
| 173 | + (np.sqrt(ds_out.meltRateUncertSqr)).where(mask) |
| 174 | + ds_out.meltRateUncertainty.attrs = uncert_attrs |
| 175 | + ds_out = ds_out.drop_vars(['meltRateUncertSqr']) |
| 176 | + ds_out.to_netcdf(f'{out_prefix}_{out_grid_name}.{date}.nc') |
| 177 | + |
| 178 | + |
| 179 | +def main(): |
| 180 | + prefix = 'Paolo_2023_iceshelf_melt_rates_1992-2017_v1.0' |
| 181 | + date = '20240220' |
| 182 | + |
| 183 | + orig_filename = 'Paolo_2023_ANT_G1920V01_IceShelfMelt.nc' |
| 184 | + processed_filename = f'{prefix}.{date}.nc' |
| 185 | + |
| 186 | + download_paolo(orig_filename) |
| 187 | + process_paolo(orig_filename, processed_filename) |
| 188 | + remap_paolo(processed_filename, prefix, date) |
| 189 | + |
| 190 | + |
| 191 | +if __name__ == '__main__': |
| 192 | + main() |
0 commit comments