Skip to content

Latest commit

 

History

History
995 lines (800 loc) · 35.6 KB

GL_iceberg_melt.org

File metadata and controls

995 lines (800 loc) · 35.6 KB
Raw

Table of contents

Introduction

Data

Printout

<xarray.Dataset> Size: 177MB
Dimensions:      (region: 7, time: 12, longitude: 720, latitude: 360)
Coordinates:
  * region       (region) int8 7B 1 2 3 4 5 6 7
  * time         (time) int8 12B 1 2 3 4 5 6 7 8 9 10 11 12
  * longitude    (longitude) float64 6kB -179.8 -179.2 -178.8 ... 179.2 179.8
  * latitude     (latitude) float64 3kB -89.75 -89.25 -88.75 ... 89.25 89.75
Data variables:
    melt         (region, time, latitude, longitude) float64 174MB ...
    melt_GL      (latitude, longitude) float64 2MB ...
    region_map   (latitude, longitude) int8 259kB ...
    region_name  (region) <U2 56B ...
    spatial_ref  int8 1B ...
Attributes:
    geospatial_lat_min:  -90
    geospatial_lat_max:  90
    geospatial_lon_min:  -180
    geospatial_lon_max:  180
    date_created:        20241117T144923Z
    title:               Normalised iceberg melt climatology per region of ca...
    history:             Processed for Schmidt (YYYY; in prep); by Ken Mankoff
    Conventions:         CF-1.8
    DOI:                 https://doi.org/10.5281/zenodo.14020895

Information

Warning

Iceberg tracks come from Marson (2024) http://doi.org/10.1029/2023jc020697. 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)

Figure

100% 8/8 [03:41<00:00, 27.66s/it]

./fig/GL_berg_melt.png

Processing

Provenance

Data from Marson (2024) http://doi.org/10.1029/2023jc020697

Tests and checks

Algorithm demonstration

  • Find the first and last iceberg track indices, derive melt from that, check inputs equal outputs
  • For some reason, last time often has large melt. Maybe we shouldn’t include it. TBD.
# synthetic ice mass array, dimesions [x=time, y=mass]
mass = np.zeros((5,5)) * np.nan
mass[0,0] = 5
mass[1,[1,2]] = [10,9]
mass[2,[1,2]] = [10,5]
mass[3,[2,3,4]] = [20,10,5]
mass[4,[3,4]] = [1.1,1]
print('mass', mass)

# flag the first time the iceberg appears
first = (~np.isnan(mass) & np.isnan(np.roll(mass,1,axis=1)))
print('first', first)

# flag the last time the iceberg appears
last = (~np.isnan(mass) & np.isnan(np.roll(mass,-1,axis=1)))
print('last', last)

# the water mass is just the derivative of the ice mass in time
h2o = -1 * np.diff(mass, axis=1, append=np.nan)
h2o[last] = mass[last]
print('h2o', h2o)

h2o_nolast = -1 * np.diff(mass, axis=1, append=np.nan) # NO LAST
print('h2o no last', h2o_nolast)

print('ice mass: ', mass[first].sum(), mass[first])
print('water mass: ', np.nansum(h2o), np.nansum(h2o,axis=1))
print('water mass no last: ', np.nansum(h2o_nolast), np.nansum(h2o_nolast,axis=1))

mass [[ 5.   nan  nan  nan  nan]
 [ nan 10.   9.   nan  nan]
 [ nan 10.   5.   nan  nan]
 [ nan  nan 20.  10.   5. ]
 [ nan  nan  nan  1.1  1. ]]
first [[ True False False False False]
 [False  True False False False]
 [False  True False False False]
 [False False  True False False]
 [False False False  True False]]
last [[ True False False False False]
 [False False  True False False]
 [False False  True False False]
 [False False False False  True]
 [False False False False  True]]
h2o [[ 5.   nan  nan  nan  nan]
 [ nan  1.   9.   nan  nan]
 [ nan  5.   5.   nan  nan]
 [ nan  nan 10.   5.   5. ]
 [ nan  nan  nan  0.1  1. ]]
h2o no last [[ nan  nan  nan  nan  nan]
 [ nan  1.   nan  nan  nan]
 [ nan  5.   nan  nan  nan]
 [ nan  nan 10.   5.   nan]
 [ nan  nan  nan  0.1  nan]]
ice mass:  46.1 [ 5.  10.  10.  20.   1.1]
water mass:  46.1 [ 5.  10.  10.  20.   1.1]
water mass no last:  21.1 [ 0.   1.   5.  15.   0.1]

Tests on real data

Load a subset

import xarray as xr
import numpy as np

root = "~/data/Marson_2024/"

mass = xr.open_mfdataset(root+'from_email/mass_01.nc')
bits = xr.open_mfdataset(root+'from_email/mass_of_bits_01.nc')
scale = xr.open_mfdataset(root+'from_email/mass_scaling_01.nc')

# xarray needs things named the same in order to multiply them together.
bits = bits.rename({'mass_of_bits':'mass'})
scale = scale.rename({'mass_scaling':'mass'})

ds = xr.merge([(mass+bits)*scale])
ds = ds.rename({'timestep':'time'})

# %time ds = ds.isel({'particle':np.arange(1000), 'time':np.arange(1000)}).load()

ds['time'].attrs['calendar'] = 'noleap'
ds['time'].attrs['units'] = 'days since 2000-01-01'
ds['time'] = pd.Timestamp("2000-01-01") + pd.to_timedelta(ds['time'].values, unit='D')
ds['particle'] = ds['particle'].astype(np.int32)

print(ds)

<xarray.Dataset> Size: 467MB
Dimensions:   (time: 5840, particle: 10000)
Coordinates:
  * time      (time) datetime64[ns] 47kB 2000-01-02 2000-01-03 ... 2000-12-31
  * particle  (particle) int32 40kB 117 118 128 129 ... 205896 205897 205916
Data variables:
    mass      (particle, time) float64 467MB dask.array<chunksize=(910, 531), meta=np.ndarray>

Initial ice mass should equal final water mass

# flag the first time the iceberg appears
first = (~np.isnan(ds['mass'].values) & np.isnan(np.roll(ds['mass'].values,1,axis=1)))
last = (~np.isnan(ds['mass'].values) & np.isnan(np.roll(ds['mass'].values,-1,axis=1)))

# the water mass is just the derivative of the ice mass in time
h2o = -1 * np.diff(ds['mass'].values, axis=1, append=np.nan)
h2o[last] = ds['mass'].values[last]

h2oX = -1 * np.diff(ds['mass'].values, axis=1, append=np.nan) # NO LAST

ds['h2o'] = (('particle','time'), h2o)
ds['h2oX'] = (('particle','time'), h2oX)
print(ds)


ice_mass = ds['mass'].values[first].sum(); print('ice mass: ', ice_mass)
water_mass = ds['h2o'].sum().values; print('water mass: ', water_mass)
print(f'diff: {water_mass - ice_mass} ({water_mass / ice_mass * 100} %)')
water_mass_X = ds['h2oX'].sum().values; print('water mass X: ', water_mass_X)
print(f'diff: {water_mass_X - ice_mass} ({water_mass_X / ice_mass * 100} %)')

<xarray.Dataset> Size: 1GB
Dimensions:   (time: 5840, particle: 10000)
Coordinates:
  * time      (time) datetime64[ns] 47kB 2000-01-02 2000-01-03 ... 2000-12-31
  * particle  (particle) int32 40kB 117 118 128 129 ... 205896 205897 205916
Data variables:
    mass      (particle, time) float64 467MB dask.array<chunksize=(910, 531), meta=np.ndarray>
    h2o       (particle, time) float64 467MB nan nan nan nan ... nan nan nan nan
    h2oX      (particle, time) float64 467MB nan nan nan nan ... nan nan nan nan
ice mass:  1886583699309968.5
water mass:  1886583699309959.8
diff: -8.75 (99.99999999999953 %)
water mass X:  1410154213770365.2
diff: -476429485539603.25 (74.74644322889779 %)

From the above, the algorithm appears to work, and water mass, computed from derivative of ice mass, matches. If we drop the last time which is always artificially large we lose 25 % of mass. However, because this is a WEIGHTED MASK, not a flux product, the magnitude does not matter, and the weights are probably more realistic by removing the last element. Unless icebergs catastrophically fail at the end and should have a large melt pulse at their final location.

Per Marson (2021) http://doi.org/10.1029/2021jc017542

The annual mass loss (hereafter referred as discharge) from the Greenland Ice Sheet (GrIS) is currently estimated to be around 1,100 Gt/yr, half of which is attributed to liquid runoff and the other half to solid discharge (Bamber et al., 2012, 2018)

Greenland discharge was provided by Bamber et al. (2012) on a 5 × 5 km grid and was remapped to the ANHA4 grid. According to the averages estimated in Bamber et al. (2012), we divided the total discharge into 46% liquid runoff and 54% solid discharge.

So discharge should be ~1100*0.54 = 594 Gt/yr

Greenland ROIs

grass ./G_3413/PERMANENT
g.mapset PERMANENT
v.import input=${DATADIR}/Mouginot_2019/GL_regions.gpkg output=ROIs
v.db.select map=ROIs
v.to.rast input=ROIs output=ROIs use=attr attribute_column=cat_

Load data

In addition to loading the public data from Marson (2024) http://doi.org/10.1029/2023jc020697 we need to add in the bergy bits (personal communication). Also, the provided mass is particles (groups of bergs) and needs to be scaled by Martin (2010) http://doi.org/10.1016/j.ocemod.2010.05.001 Table 1 to convert particle mass to ice mass.

import xarray as xr
import pandas as pd
import numpy as np

root='~/data/Marson_2024/'

lon = xr.open_mfdataset(root+'lon_*.nc', join='override', concat_dim='particle', combine='nested')
lat = xr.open_mfdataset(root+'lat_*.nc', join='override', concat_dim='particle', combine='nested')
mass = xr.open_mfdataset([root+'from_email/mass_01.nc',
                          root+'from_email/mass_02.nc',
                          root+'from_email/mass_03.nc',
                          root+'from_email/mass_04.nc'],
                         join='override', concat_dim='particle', combine='nested')
bits = xr.open_mfdataset(root+'from_email/mass_of_bits_*.nc', join='override', concat_dim='particle', combine='nested')
scale = xr.open_mfdataset(root+'from_email/mass_scaling_*.nc', join='override', concat_dim='particle', combine='nested')

# xarray needs things named the same in order to multiply them together.
bits = bits.rename({'mass_of_bits':'mass'})
scale = scale.rename({'mass_scaling':'mass'})

%time ds = xr.merge([lon,lat,(mass+bits)*scale])

ds = ds.rename({'timestep':'time'})
ds['time'].attrs['calendar'] = 'noleap'
ds['time'].attrs['units'] = 'days since 2000-01-01'
ds['time'] = pd.Timestamp("2000-01-01") + pd.to_timedelta(ds['time'].values, unit='D')
ds['particle'] = ds['particle'].astype(np.int32)

print(ds)

CPU times: user 6.31 ms, sys: 0 ns, total: 6.31 ms
Wall time: 6.31 ms
<xarray.Dataset> Size: 5GB
Dimensions:   (time: 5840, particle: 34025)
Coordinates:
  * time      (time) datetime64[ns] 47kB 2000-01-02 2000-01-03 ... 2000-12-31
  * particle  (particle) int32 136kB 117 118 128 129 ... 1806577 1806831 1807085
Data variables:
    lon       (particle, time) float64 2GB dask.array<chunksize=(1667, 974), meta=np.ndarray>
    lat       (particle, time) float64 2GB dask.array<chunksize=(1667, 974), meta=np.ndarray>
    mass      (particle, time) float64 2GB dask.array<chunksize=(910, 531), meta=np.ndarray>

Compute mass loss

# flag the first time the iceberg appears
first = (~np.isnan(ds['mass'].values) & np.isnan(np.roll(ds['mass'].values,1,axis=1)))

# the water mass is just the derivative of the ice mass in time
h2o = -1 * np.diff(ds['mass'].values, axis=1, append=np.nan) # drop last time as diff does naturally

ds['h2o'] = (('particle','time'), h2o)
ds['first'] = (('particle','time'), first)
print(ds)

<xarray.Dataset> Size: 7GB
Dimensions:   (time: 5840, particle: 34025)
Coordinates:
  * time      (time) datetime64[ns] 47kB 2000-01-02 2000-01-03 ... 2000-12-31
  * particle  (particle) int32 136kB 117 118 128 129 ... 1806577 1806831 1807085
Data variables:
    lon       (particle, time) float64 2GB dask.array<chunksize=(1667, 974), meta=np.ndarray>
    lat       (particle, time) float64 2GB dask.array<chunksize=(1667, 974), meta=np.ndarray>
    mass      (particle, time) float64 2GB dask.array<chunksize=(910, 531), meta=np.ndarray>
    h2o       (particle, time) float64 2GB nan nan nan nan ... nan nan nan nan
    first     (particle, time) bool 199MB False False False ... False False True

Save snapshot

comp = dict(zlib=True, complevel=2)
encoding = {var: comp for var in ds.data_vars}

delayed_obj = ds.to_netcdf('tmp/bergs.nc', encoding=encoding, compute=False)
from dask.diagnostics import ProgressBar
with ProgressBar():
    results = delayed_obj.compute()

# saves as 250 MB file. Takes a few minutes...

[########################################] | 100% Completed | 231.73 s

Load snapshot

import xarray as xr
import numpy as np
import pandas as pd

%time ds = xr.open_dataset('tmp/bergs.nc').load() # load everything into memory
# Takes a while...

CPU times: user 13.9 s, sys: 18 s, total: 31.8 s
Wall time: 43.3 s

Test

%time ice_mass = ds['mass'].values[ds['first'].values].sum()
print('ice mass: ', ice_mass * 1E-12 / 16) # total kg over 16 years -> Gt/yr
%time water_mass = np.nansum(ds['h2o'].values)
print('water mass: ', water_mass * 1E-12 / 16)

CPU times: user 92.9 ms, sys: 0 ns, total: 92.9 ms
Wall time: 91.4 ms
ice mass:  407.2388163829433
CPU times: user 996 ms, sys: 4.79 s, total: 5.79 s
Wall time: 7.51 s
water mass:  296.9381105981331

The difference between the Marson (2024) http://doi.org/10.1029/2023jc020697 407 Gt/year and the Mankoff (2020) http://doi.org/10.5194/essd-12-1367-2020 ~500 Gt/year (subject to change with each version) is not important. It can represent a lot of things, most likely that Mankoff (2020) is discharge across flux gates upstream from the terminus, so 100 - 407/500 % = 18.6 % is submarine melt, and the remainder is the Marson icebergs.

Additional melting occurs in the fjord and must be handled if the model does not resolve fjords.

This product should be shared as one and several weighted masks that sum to 1, and then users can scale by their own estimated discharge.

The difference between ice inputs and water outputs is described above - we drop the last melt event of each berg which seems anomalously high.

Iceberg meltwater locations

Export each particle to file

  • Warning: 34k files generated here.
from tqdm import tqdm
for p in tqdm(range(ds['particle'].values.size)):
    df = ds.isel({'particle':p})\
           .to_dataframe()\
           .dropna()
    if df.size == 0: continue
    df['first'] = df['first'].astype(int)
    df[['particle','lon','lat','mass','h2o','first']]\
        .to_csv(f"./Marson_2024_tmp/{str(p).zfill(5)}.csv", header=None)

100% 34025/34025 [03:16<00:00, 173.55it/s]

Ingest each track and organize by source

Set up domain

[[ -e G_3413 ]] || grass -ec EPSG:3413 ./G_3413
grass ./G_3413/PERMANENT
g.mapset -c Marson_2024
export GRASS_OVERWRITE=1

Load ice ROIs

ogr2ogr ./tmp/Mouginot.gpkg -t_srs "EPSG:3413" ${DATADIR}/Mouginot_2019/Greenland_Basins_PS_v1.4.2.shp
v.import input=${DATADIR}/Mouginot_2019/GL_regions.gpkg output=GL
v.db.select map=GL
g.region vector=GL res=10000 -pa
v.to.rast input=GL output=GL use=attr attribute_column=cat_

Import each track and find closest ice ROI for initial location

  • Take time (month) into account.
  • 84 maps: 7 roi * 12 month: CE_01, CE_02, etc…
# reorder from "cat,id,lon,lat,ice mass,water mass" to lon,lat,water,id,time
cat Marson_2024_tmp/*.csv | awk -F, '{OFS=",";print $3,$4,$6,$7,$2,$1}' > tmp/tracks.csv

# head -n100 tmp/tracks.csv \
cat tmp/tracks.csv \
  | m.proj -i input=- separator=comma \
  | tr ' ' ',' \
  | v.in.ascii -n input=- output=bergs sep=, \
               columns='x double,y double,water double,first int,id int,time varchar(10)'

g.region vector=bergs res=25000 -pa
g.region save=iceberg_region

r.mapcalc "x = x()"
r.mapcalc "y = y()"
r.mapcalc "area = area()"

# Record nearest region at all times, by finding the region nearest the 1st time
v.db.addcolumn map=bergs columns="roi VARCHAR(3)"

# v.db.select map=bergs|head

v.extract input=bergs where='(first == 1)' output=t0
v.distance from=t0 to=GL upload=to_attr to_column=label column=roi
db.select table=t0|head| column -s"|" -t
db.select table=bergs|head| column -s"|" -t

roi=NO # debug
for roi in NO NE SE SW CW NW CE; do
  echo "Processing ROI: ${roi}"
  ids=$(db.select -c sql="select id from t0 where roi == '${roi}'")
  ids=$(echo ${ids}| tr ' ' ',')
  db.execute sql="update bergs set roi = \"${roi}\" where id in (${ids})"
done

db.select table=bergs | head -n 10 | column -s"|" -t

# # convert to raster, binned by melt per cell (a.k.a density or heat or quilt map)
# roi=NO; month=03 # debug
# # this loop takes a few minutes per ROI. Could use GNU parallel.
# for roi in NO NE SE SW CW NW CE; do
#   echo "Processing ROI: ${roi}"
#   for month in $(seq -w 12); do
#     echo "Processing month: ${month}"
#     v.out.ascii input=bergs output=- format=point columns=water \
# 		where="roi == \"${roi}\" AND time LIKE \"2000-${month}%\""\
# 		| r.in.xyz input=- z=4 output=${roi}_${month} method=sum
#     r.colors -g map=${roi}_${month} color=viridis
#     # Convert from kg/16years to kg/s
#     r.mapcalc "${roi}_${month} = ${roi}_${month} / 16 / 365 / 86400"
#   done
# done

rois="NO NE SE SW CW NW CE"
months=$(seq -w 12)
# convert to raster, binned by melt per cell (a.k.a density or heat or quilt map)
parallel -j4 --bar "v.out.ascii input=bergs output=- format=point columns=water \
		where=\"roi == '{1}' AND time LIKE '2000-{2}%'\" \
		| r.in.xyz input=- z=4 output={1}_{2} method=sum" ::: $rois ::: $months
# Convert from kg/16years to kg/s
parallel --bar "r.mapcalc \"{1}_{2} = {1}_{2} / 16 / 365 / 86400\"" ::: $rois ::: $months
parallel --bar "r.colors -g map={1}_{2} color=viridis" ::: $rois ::: $months

Sanity check: Gt/year/sector

tot=0
for roi in CE CW NE NO NW SE SW; do
  roimonths=$(g.list type=raster pattern=${roi}_?? sep=,)
  eval $(r.univar -g ${roimonths})
  # convert from kg/s to Gt/year
  roi_gt=$(echo "${sum} * 86400 * 365 * 10^(-12)" | bc -l)
  echo "${roi}: ${roi_gt}"
  tot=$(echo "${tot} + ${roi_gt}" | bc -l)
done
echo ""
echo "total: " ${tot}

My estimates of discharge by ROI?

import xarray as xr
dd = xr.open_dataset('/home/kdm/data/Mankoff_2020/ice/latest/region.nc')\
       .sel({'time':slice('2000-01-01','2019-12-31')})\
       .resample({'time':'YS'})\
       .mean()\
       .mean(dim='time')\
       ['discharge']

print(dd.sum())
dd.to_dataframe()

<xarray.DataArray 'discharge' ()> Size: 8B
array(476.48053387)
regiondischarge
CE77.8964
CW86.1499
NE25.9822
NO25.329
NW103.127
SE139.048
SW18.9477

Clean up isolated tracks and smooth ROI boundary

Display problem:

g.mapset Marson_2024
d.mon wx0
r.colors -g map=CE color=viridis
d.rast CE_01 values=0.001-1000
d.rast CE_06 values=0.001-1000

Load global coastlines for eventual cropping

r.mask -r
v.import extent=region input=${DATADIR}/NaturalEarth/ne_50m_land.shp output=land
v.to.rast input=land output=land use=value value=1
r.mask -i land

Problem: Small isolate tracks

Algorithm:

  • Smooth (if this is all we do, it would remove the high melt cells near the coast)
  • Crop by coastlines (Greenland, Canada, Iceland, Svalbard, etc.)
  • Set new area equal to old melt values where both exist (undo smooth where there was data)
  • Set new area (from smooth) equal to some low amount (median) where new area but no old melt
  • Rescale so total melt is unchanged
export GRASS_OVERWRITE=1

rois="NO NE SE SW CW NW CE"
months=$(seq -w 12)

# g.gui.mapswipe first=CE_08 second=CE_08_n

roi=NW; month=08 # debug
for roi in $rois; do
  for month in $(seq -w 12); do

    r.mapcalc "r_main = if(${roi}_${month} > 0, ${roi}_${month}, null())"
    eval $(r.univar -g -e r_main percentile=95)
    roisum=${sum}
    roimedian=${median}

    # buffer by X m to get a smooth border
    r.grow.distance -m input=r_main distance=distance
    r.mapcalc "r_buffer = distance < 250000" # 250 km
    # at 25 km resolution this is 10 cells

    # Expand original melt map with average of neighbors
    # SIZE here should take into account buffer distance above
    r.neighbors input=${roi}_${month} output=r_neighbor method=average size=15 weighting_function=gaussian weighting_factor=3

    # Convert new area to melt w/ median values where filled in
    r.mask -i land
    r.mapcalc "r_newmelt = if(r_buffer, if((${roi}_${month} > ${percentile_95}), ${roi}_${month}, r_neighbor))"
    r.mask -r

    # rescale to total melt from original map
    eval $(r.univar -g r_newmelt) # get $sum
    r.mapcalc "${roi}_${month}_smooth = (r_newmelt / ${sum}) * ${roisum}"
    r.colors -g map=${roi}_${month}_smooth,${roi}_${month} color=viridis
  done
done

# x=${roi}_${month}; g.gui.mapswipe first=$x second=r_newmelt
# x=${roi}_${month}; g.gui.mapswipe first=$x second=${x}_smooth # mode=mirror
x=NO_08; g.gui.mapswipe first=$x second=${x}_smooth mode=mirror

# # Generate GL-wide map
# r.mapcalc "GL_notail = 0"
# for roi in NO NE SE SW CW NW CE; do
#   r.mapcalc "GL_notail = GL_notail + if(isnull(${roi}_notail), 0, ${roi}_notail)"
# done
# r.null map=GL_notail setnull=0

# convert to m-2 prior to reproject
for roi in $rois; do
  for mm in $months; do
    r.mapcalc "${roi}_${mm}_m2 = ${roi}_${mm}_smooth / area"
    r.colors map=${roi}_m2 color=viridis -g
  done
done

Reproject from 3413 to 4326

  • GRASS raster reprojection uses nearest neighbor.
  • We need to convert to points in 3413, and then sum multiple 3413 points that fall within a 4326 grid cell.
  • Work at a high resolution in 3413 so that there are no cells in 4326 that are left empty.
  • If we reproject points from the current 25 km resolution there will be gaps because in N Greenland at EPSG:4326 and 0.5 degree resolution grid cells are ~15 km wide.
  • Therefore, resample to 5 km and then reproject points at that resolution.
  • Check: The sum of flux (kg) and the sum of flux (m-2) should be the same after reprojection.

3413:

# grass ./G_3413/Marson_2024

g.region res=5000 -pa # Divide value by 25 because they were calculated on a grid 5x5 larger

roi=NO; mm=09 # debug
for roi in $rois; do
  for mm in $months; do
    echo $roi $mm
    r.mapcalc --q "tmp_5km = ${roi}_${mm}_smooth / 25"
    r.out.xyz --q input=tmp_5km output=- \
      | grep -v "|0$" \
      | m.proj --q -o -d input=- > ./tmp/${roi}_${mm}.xyz
  done
done
# r.univar -gr tmp_5km  | grep -E "mean|sum"
# r.univar -gr NO_09_m2  | grep -E "mean|sum"
g.region -pa region=iceberg_region # reset region

4326:

grass ./G_4326/Marson_2024

rois="NO NE SE SW CW NW CE"
months=$(seq -w 12)
export GRASS_OVERWRITE=1

r.proj project=G_3413 input=area output=area_3413

roi=NO; mm=09 # debug
for roi in $rois; do
  for mm in $months; do
    echo $roi $mm
    r.in.xyz --q input=./tmp/${roi}_${mm}.xyz z=3 output=${roi}_${mm} method=sum
    #r.univar -g ${roi}_${mm} | grep -E "mean|sum" # matches SUM [kg] from 3413
    # r.mapcalc "${roi}_${mm}_m2 = ${roi}_${mm} / area_3413"
    # r.univar -g ${roi}_${month}_m2 | grep -E "mean|sum" # matches SUM [m-2] from 3413
  done
done

Sanity check: Gt/year/sector

tot=0
for roi in $rois; do
  roimonths=$(g.list type=raster pattern=${roi}_?? sep=,)
  eval $(r.univar -g ${roimonths})
  # convert from kg/s to Gt/year
  roi_gt=$(echo "${sum} * 86400 * 365 * 10^(-12)" | bc -l)
  echo "${roi}: ${roi_gt}"
  tot=$(echo "${tot} + ${roi_gt}" | bc -l)
done
echo ""
echo "total: " ${tot}
CE: 44.76374035974432864000
CW: 55.79715155801872896000
NE: 19.41732943638206054400
NO: 20.40425666413177809600
NW: 62.87261995622087280000
SE: 80.38475519025573360000
SW: 13.37355447290342323200
total:  297.01340763765692587200

These are the results without last berg removed:

CE: 60.83514096745521600000
CW: 77.10383458384345872000
NE: 25.40014065127645070400
NO: 28.68057265378820385600
NW: 85.30203183944359392000
SE: 111.19281504778338912000
SW: 18.72428302446387782400

total:  407.23881876805419014400

Export to NetCDF

import numpy as np
import xarray as xr
import rioxarray as rxr
from tqdm import tqdm
import datetime

from grass_session import Session
from grass.script import core as gcore
import grass.script as gscript
# import grass.script.setup as gsetup
# import grass python libraries
from grass.pygrass.modules.shortcuts import general as g
from grass.pygrass.modules.shortcuts import raster as r
from grass.pygrass.modules.shortcuts import vector as v
from grass.pygrass.modules.shortcuts import temporal as t
from grass.script import array as garray

S = Session()
S.open(gisdb=".", location="G_4326", mapset="Marson_2024", create_opts=None)
lon = garray.array("x")[::-1,:]
lat = garray.array("y")[::-1,:]

melt = np.zeros((7, 12, lon.shape[0], lat.shape[1]))
for roi_i,roi in enumerate(['CE','CW','NE','NO','NW','SE','SW']):
    for month in range(12):
        mm = str(month+1).zfill(2)
        melt[roi_i,month,:,:] = garray.array(roi+'_'+mm)[::-1,:]

ds = xr.Dataset({
    'melt': xr.DataArray(data = melt,
                         dims = ['region','time','latitude','longitude'],
                         coords = {'region':np.arange(7).astype(np.int8)+1,
                                   'time': np.arange(12).astype(np.int8)+1,
                                   'longitude':lon[0,:],
                                   'latitude':lat[:,0]})})

ds['melt_GL'] = ds['melt'].sum(dim=['region','time'])

# normalize from kg m-2 to m-2
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['melt'] = ds['melt'] / ds['melt_GL'].sum() / cell_area
ds['melt_GL'] = ds['melt_GL'] / ds['melt_GL'].sum() / cell_area

ROIs = garray.array("ROIs")[::-1,:]
ds['region_map'] = (('latitude','longitude'), ROIs.astype(np.byte))

S.close() # Done with GRASS

ds['region_name'] = (('region'), ['CE','CW','NE','NO','NW','SE','SW'])

ds = ds.rio.write_crs('epsg:4326')
ds['spatial_ref'] = ds['spatial_ref'].astype(np.byte)
ds = ds.rio.set_spatial_dims('longitude','latitude')

ds['latitude'].attrs['long_name'] = 'latitude'
ds['latitude'].attrs['axis'] = 'Y'
ds['latitude'].attrs['units'] = 'degrees_north'
ds['latitude'].attrs['standard_name'] = 'latitude'
ds['longitude'].attrs['long_name'] = 'longitude'
ds['longitude'].attrs['axis'] = 'X'
ds['longitude'].attrs['units'] = 'degrees_east'
ds['longitude'].attrs['standard_name'] = 'longitude'

ds['time'].attrs['standard_name'] = 'time'

ds['melt'].attrs['long_name'] = 'Normalised iceberg melt climatology per region of calving'
ds['melt'].attrs['units'] = 'm-2'
ds['melt'].attrs['standard_name'] = 'water_flux_into_sea_water_from_icebergs'
ds['melt'].attrs['comment'] = 'This value summed by area and time and multiplied by cell area should sum to 1'
ds['melt_GL'].attrs['long_name'] = 'Normalised iceberg melt climatology for all Greenland'
ds['melt_GL'].attrs['units'] = 'm-2'
ds['melt_GL'].attrs['standard_name'] = 'water_flux_into_sea_water_from_icebergs'
ds['melt_GL'].attrs['comment'] = 'This value multiplied by cell area should sum to 1'

ds['region'].attrs['long_name'] = 'Region IDs'
ds['region_name'].attrs['long_name'] = 'Mouginot (2019) region'
ds['region_name'].attrs['standard_name'] = 'region'
ds['region_map'].attrs['long_name'] = 'Region IDs'

ds['spatial_ref'].attrs['horizontal_datum_name'] = 'WGS 84'

ds.attrs['geospatial_lat_min'] = -90 # ds['latitude'].values.min()
ds.attrs['geospatial_lat_max'] = 90 # ds['latitude'].values.max()
ds.attrs['geospatial_lon_min'] = -180 # ds['longitude'].values.min()
ds.attrs['geospatial_lon_max'] = 180 # 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 per region of calving from Marson (2024)'
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=2) # Internal NetCDF compression
encoding = {var: comp for var in ds.drop_vars('region_name').data_vars}

!rm ./dat/GL_iceberg_melt.nc
ds.to_netcdf('./dat/GL_iceberg_melt.nc', encoding=encoding)
print(ds)

<xarray.Dataset> Size: 177MB
Dimensions:      (region: 7, time: 12, longitude: 720, latitude: 360)
Coordinates:
  * region       (region) int8 7B 1 2 3 4 5 6 7
  * time         (time) int8 12B 1 2 3 4 5 6 7 8 9 10 11 12
  * longitude    (longitude) float64 6kB -179.8 -179.2 -178.8 ... 179.2 179.8
  * latitude     (latitude) float64 3kB -89.75 -89.25 -88.75 ... 89.25 89.75
    spatial_ref  int8 1B 0
Data variables:
    melt         (region, time, latitude, longitude) float64 174MB 0.0 ... 0.0
    melt_GL      (latitude, longitude) float64 2MB 0.0 0.0 0.0 ... 0.0 0.0 0.0
    region_map   (latitude, longitude) int8 259kB 0 0 0 0 0 0 0 ... 0 0 0 0 0 0
    region_name  (region) <U2 56B 'CE' 'CW' 'NE' 'NO' 'NW' 'SE' 'SW'
Attributes:
    geospatial_lat_min:  -90
    geospatial_lat_max:  90
    geospatial_lon_min:  -180
    geospatial_lon_max:  180
    date_created:        20241117T144923Z
    title:               Normalised iceberg melt climatology per region of ca...
    history:             Processed for Schmidt (YYYY; in prep); by Ken Mankoff
    Conventions:         CF-1.8
    DOI:                 https://doi.org/10.5281/zenodo.14020895

Units check

import xarray as xr
import numpy as np

ds = xr.open_dataset('dat/GL_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_GL', (ds['melt_GL']*ds['area']).sum(dim=['latitude','longitude']) )

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_GL <xarray.DataArray ()> Size: 8B
array(1.)
melt times [0.05709769 0.05647682 0.05988504 0.0497022  0.05714598 0.07278505
 0.12313315 0.15238323 0.13335883 0.09943356 0.07442078 0.06417767] 1.0
melt rois [0.15071286 0.18786072 0.06537526 0.0686981  0.21168277 0.27064352
 0.04502677] 0.9999999999999999