- Data from Nico Jourdain processed from Mathiot (2023)
- Rework to match other products generated for the Schmidt freshwater paper
<xarray.Dataset> Size: 226MB
Dimensions: (time: 12, region: 18, latitude: 360, longitude: 720)
Coordinates:
* longitude (longitude) float64 6kB -179.8 -179.2 ... 179.2 179.8
* latitude (latitude) float64 3kB -89.75 -89.25 ... 89.25 89.75
* region (region) int32 72B 1 2 3 4 5 6 7 ... 13 14 15 16 17 18
* time (time) int8 12B 1 2 3 4 5 6 7 8 9 10 11 12
spatial_ref int8 1B ...
Data variables:
melt (time, region, latitude, longitude) float32 224MB ...
msk_nemo (latitude, longitude) float32 1MB ...
region_name (region) <U5 360B ...
region_map (latitude, longitude) int16 518kB ...
region_map_expanded (latitude, longitude) int16 518kB ...
Attributes: (12/13)
description: Annual JRA55 climatology
original_data: NEMO 0.25° simulations by Anna Olive-Abello (in prep...
script_used: remap_per_basin.py
processed_by: Nicolas Jourdain <nicolas.jourdain@univ-grenoble-alp...
geospatial_lat_min: -89.75
geospatial_lat_max: 89.75
... ...
geospatial_lon_max: 179.75
date_created: 20250220T153857Z
title: Normalised iceberg melt climatology in the Southern ...
history: Processed for Schmidt (YYYY; in prep); by Ken Mankoff
Conventions: CF-1.8
DOI: https://doi.org/10.5281/zenodo.14020895
Warning
The highest melt rates (largest meltwater injection) occurs in near-coastal cells. If your model has land covering some of these cells, you may lose large melt inputs. Rescaling the melt so melt*area sums to 1 for your ocean is a good idea, but this also redistributes the largest melt points over the entire melt region. It may be better to re-scale the melt by increasing only the largest cell or largest few cells (which are hopefully nearby the coast and the high melt rate cells that were covered by land)
100% 19/19 [29:50<00:00, 94.22s/it]
import rioxarray as rxr
import rasterio as rio
import xarray as xr
import numpy as np
import datetime
ds = xr.open_dataset('~/data/Mathiot_2023/iceberg_melt_pattern_SH_per_basin.nc')
# add projection metadata
ds = ds.rio.write_crs('epsg:4326') # create ds['spatial_ref']
ds = ds.rio.set_spatial_dims(x_dim='longitude', y_dim='latitude') # or ('lon','lat') and only maybe needed
ds['spatial_ref'] = ds['spatial_ref'].astype(np.byte)
# provide cell center values at all coordinates
ds = ds.pad(latitude=(0, 1), longitude=(0,1)) # Add one column at the end
ds['latitude'] = np.linspace(-89.75, 89.75, num=360)
ds['longitude'] = np.linspace(-179.75, 179.75, num=720)
ds['time'] = (('time'), np.arange(12).astype(np.int8)+1)
# Rignot basins are 1 through 18, not 0 through 17
# ds = ds.rename({'basin':'region'})
ds['region'] = (ds['region']).astype(np.int32)
ds['region_name'] = (('region'), ['A-Ap', 'Ap-B', 'B-C', 'C-Cp', 'Cp-D',
'D-Dp', 'Dp-E', 'E-Ep', 'Ep-F', 'F-G',
'G-H', 'H-Hp', 'Hp-I', 'I-Ipp', 'Ipp-J',
'J-Jpp', 'Jpp-K', 'K-A'])
ds.attrs['description'] = 'Annual JRA55 climatology'
ds = ds.rename_vars({'melt_pattern':'melt'})
# 'pattern_SH_allbasins':'melt_AQ'})
ds['melt'].attrs['units'] = 'm-2'
ds['melt'].attrs['grid_mapping'] = 'spatial_ref'
ds['melt'].attrs['standard_name'] = 'water_flux_into_sea_water_from_icebergs'
ds['melt'].attrs['long_name'] = 'Normalised iceberg melt climatology per region of calving'
rt = rio.open('./tmp/regions.tif').read(1)[::-1,:]
rt[rt < 0] = 0
ds['region_map'] = (('latitude','longitude'), rt)
rt = rio.open('./tmp/regions_expanded.tif').read(1)[::-1,:]
rt[rt < 0] = 0
ds['region_map_expanded'] = (('latitude','longitude'), rt)
ds['spatial_ref'].attrs['horizontal_datum_name'] = 'WGS 84'
ds['region'].attrs['long_name'] = 'Region IDs'
ds['time'].attrs['standard_name'] = 'time'
ds['longitude'].attrs['standard_name'] = 'longitude'
ds['longitude'].attrs['long_name'] = 'longitude'
ds['longitude'].attrs['axis'] = 'X'
ds['longitude'].attrs['units'] = 'degrees_east'
ds['latitude'].attrs['standard_name'] = 'latitude'
ds['latitude'].attrs['long_name'] = 'latitude'
ds['latitude'].attrs['axis'] = 'Y'
ds['latitude'].attrs['units'] = 'degrees_north'
ds['region_map'].attrs['long_name'] = 'IMBIE regions'
ds['region_name'].attrs['long_name'] = 'IMBIE regions'
ds['region_name'].attrs['standard_name'] = 'region'
ds.attrs['geospatial_lat_min'] = ds['latitude'].values.min()
ds.attrs['geospatial_lat_max'] = ds['latitude'].values.max()
ds.attrs['geospatial_lon_min'] = ds['longitude'].values.min()
ds.attrs['geospatial_lon_max'] = ds['longitude'].values.max()
ds.attrs['date_created'] = datetime.datetime.now(datetime.timezone.utc).strftime("%Y%m%dT%H%M%SZ")
ds.attrs['title'] = 'Normalised iceberg melt climatology in the Southern Hemisphere per month and region of calving'
ds.attrs['history'] = 'Processed for Schmidt (YYYY; in prep); by Ken Mankoff'
ds.attrs['Conventions'] = 'CF-1.8'
ds.attrs['DOI'] = 'https://doi.org/10.5281/zenodo.14020895'
comp = dict(zlib=True, complevel=5)
encoding = {var: comp for var in ds.drop_vars(['region_name']).data_vars}
!rm ./dat/AQ_iceberg_melt.nc
ds.to_netcdf('./dat/AQ_iceberg_melt.nc', encoding=encoding)
!ncdump -h ./dat/AQ_iceberg_melt.nc
netcdf AQ_iceberg_melt {
dimensions:
time = 12 ;
region = 18 ;
latitude = 360 ;
longitude = 720 ;
variables:
float melt(time, region, latitude, longitude) ;
melt:_FillValue = NaNf ;
melt:long_name = "Normalised iceberg melt climatology per region of calving" ;
melt:comment = "The spatial integral on the spherical Earth summed over the 12 months and all regions is equal to 1.0" ;
melt:units = "m-2" ;
melt:grid_mapping = "spatial_ref" ;
melt:standard_name = "water_flux_into_sea_water_from_icebergs" ;
melt:coordinates = "spatial_ref" ;
float msk_nemo(latitude, longitude) ;
msk_nemo:_FillValue = NaNf ;
msk_nemo:long_name = "Original land/sea mask in the NEMO simulation" ;
msk_nemo:coordinates = "spatial_ref" ;
string region_name(region) ;
region_name:long_name = "IMBIE regions" ;
region_name:standard_name = "region" ;
region_name:coordinates = "spatial_ref" ;
double longitude(longitude) ;
longitude:_FillValue = NaN ;
longitude:standard_name = "longitude" ;
longitude:long_name = "longitude" ;
longitude:axis = "X" ;
longitude:units = "degrees_east" ;
double latitude(latitude) ;
latitude:_FillValue = NaN ;
latitude:standard_name = "latitude" ;
latitude:long_name = "latitude" ;
latitude:axis = "Y" ;
latitude:units = "degrees_north" ;
int region(region) ;
region:long_name = "Region IDs" ;
region:comment = "IMBIE2 basin (https://doi.org/10.1038/s41586-018-0179-y)" ;
byte time(time) ;
time:standard_name = "time" ;
byte spatial_ref ;
spatial_ref:crs_wkt = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]" ;
spatial_ref:semi_major_axis = 6378137. ;
spatial_ref:semi_minor_axis = 6356752.31424518 ;
spatial_ref:inverse_flattening = 298.257223563 ;
spatial_ref:reference_ellipsoid_name = "WGS 84" ;
spatial_ref:longitude_of_prime_meridian = 0. ;
spatial_ref:prime_meridian_name = "Greenwich" ;
spatial_ref:geographic_crs_name = "WGS 84" ;
spatial_ref:horizontal_datum_name = "WGS 84" ;
spatial_ref:grid_mapping_name = "latitude_longitude" ;
spatial_ref:spatial_ref = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]" ;
short region_map(latitude, longitude) ;
region_map:long_name = "IMBIE regions" ;
region_map:coordinates = "spatial_ref" ;
short region_map_expanded(latitude, longitude) ;
region_map_expanded:coordinates = "spatial_ref" ;
// global attributes:
:description = "Annual JRA55 climatology" ;
string :original_data = "NEMO 0.25° simulations by Anna Olive-Abello (in preparation)" ;
:script_used = "remap_per_basin.py" ;
:processed_by = "Nicolas Jourdain <nicolas.jourdain@univ-grenoble-alpes.fr>" ;
:geospatial_lat_min = -89.75 ;
:geospatial_lat_max = 89.75 ;
:geospatial_lon_min = -179.75 ;
:geospatial_lon_max = 179.75 ;
:date_created = "20250220T153857Z" ;
:title = "Normalised iceberg melt climatology in the Southern Hemisphere per month and region of calving" ;
:history = "Processed for Schmidt (YYYY; in prep); by Ken Mankoff" ;
:Conventions = "CF-1.8" ;
:DOI = "https://doi.org/10.5281/zenodo.14020895" ;
}
import xarray as xr
import numpy as np
ds = xr.open_dataset('dat/AQ_iceberg_melt.nc')
llon,llat = np.meshgrid(ds['longitude'].values, ds['latitude'].values)
earth_rad = 6.371e6 # Earth radius in m
resdeg = 0.5 # output grid resolution in degrees
cell_area = np.cos(np.deg2rad(llat)) * earth_rad**2 * np.deg2rad(resdeg)**2
ds['area'] = (('latitude','longitude'), cell_area)
# print(ds)
print( 'melt', (ds['melt']*ds['area']).sum().values )
times = (ds['melt']*ds['area']).sum(dim=['latitude','longitude','region'])
print( 'melt times', times.values, times.sum().values)
rois = (ds['melt']*ds['area']).sum(dim=['latitude','longitude','time'])
print( 'melt rois', rois.values, rois.sum().values)
melt 0.989742102609713 melt times [0.21008673 0.20930613 0.1380628 0.07019004 0.0393484 0.03732747 0.02741233 0.02707993 0.02636447 0.03068494 0.05382992 0.12004895] 0.9897421026097116 melt rois [0.04957716 0.02489026 0.04744073 0.06769404 0.08279573 0.06925033 0.00440376 0.04834739 0.09674372 0.06607943 0.12372214 0.01435447 0.01112696 0.03560035 0.01393721 0.14893319 0.04256935 0.04227588] 0.9897421026097117