Large Area Processing - Polar Hackathon

This notebook shows how you can run larger scale analysis on the cube. You will have to update the code for your own use case.

Since the cube is too large to load into memory, the notebook uses Dask reductions, preserving the native Zarr chunks, and returns small tables.

This notebook is designed to run on a cloud workspace where Dask Gateway is deployed and configured. The Gateway cells below will raise ValueError: No dask-gateway address provided or found in configuration on a local machine or any environment without Dask Gateway configuration.

If you are running locally, skip the Dask Gateway cells and run the local Dask cluster alternative instead. Local runs may need a smaller region or fewer datasets depending on your machine.

import dask.array as da
import dask.distributed
import warnings
import numpy as np
import xarray as xr
from dask.array.core import PerformanceWarning
from dask_gateway import Gateway

warnings.filterwarnings("ignore", message="In a future version of xarray the default value for join.*", category=FutureWarning)
warnings.filterwarnings("ignore", message="Increasing number of chunks.*", category=PerformanceWarning)

## run this line if you get a file system not found error!
# import os
# os.chdir("/tmp")
# print("cwd:", os.getcwd())

Dask Cluster Setup¶

Use the Dask Gateway cells on the cloud workspace. For local runs, skip the Gateway cells and run the local alternative cell instead.

# Local alternative: run this cell instead of the Dask Gateway cells below.
cluster = dask.distributed.LocalCluster(
    n_workers=2,
    threads_per_worker=2,
    memory_limit="4GB",
)
client = dask.distributed.Client(cluster)
client

Cloud Workspace¶

gateway = Gateway()
cluster_options = gateway.cluster_options()
cluster_options

cluster = gateway.new_cluster(cluster_options=cluster_options)
cluster.scale(2)
cluster

client = cluster.get_client()
client

Connect with Data Cube¶

cube_paths = [
    "https://s3.waw4-1.cloudferro.com/EarthCODE/OSCAssets/antarctica_cube/icetemp.zarr",
    "https://s3.waw4-1.cloudferro.com/EarthCODE/OSCAssets/antarctica_cube/sec.zarr",
    "https://s3.waw4-1.cloudferro.com/EarthCODE/OSCAssets/antarctica_cube/antarctica-combined.zarr",
    "https://s3.waw4-1.cloudferro.com/EarthCODE/OSCAssets/antarctica_cube/icemask_composite.zarr/",
    "https://s3.waw4-1.cloudferro.com/EarthCODE/OSCAssets/antarctica_cube/ice_velocity.zarr",
]

ds = xr.open_mfdataset(cube_paths, engine="zarr", chunks={}, compat="no_conflicts")
ds

dx = abs((ds.x.isel(x=1) - ds.x.isel(x=0)).load().item())
dy = abs((ds.y.isel(y=1) - ds.y.isel(y=0)).load().item())
km2 = dx * dy / 1e6

ds[[
    "bedrock_topography_mask",
    "bedrock_topography_bed",
    "bedrock_topography_thickness",
    "bedrock_topography_errbed",
    "ice_shelf_basal_melt_rate",
    "groundlines_mask",
    "subglacial_lakes_mask",
    "supra_glacial_lakes_mask",
]]

Area of the main surface classes¶

Mask codes: 0 ocean, 1 ice-free land, 2 grounded ice, 3 floating ice, 4 Lake Vostok.

mask = ds["bedrock_topography_mask"]

surface_area = xr.Dataset({
    "ocean": (mask == 0).sum() * km2,
    "ice_free_land": (mask == 1).sum() * km2,
    "grounded_ice": (mask == 2).sum() * km2,
    "floating_ice": (mask == 3).sum() * km2,
    "lake_vostok": (mask == 4).sum() * km2,
}).compute()

surface_area.to_array().to_series().rename("area_km2").to_frame()

Example 1: marine-based grounded ice exposure¶

This estimates how much grounded ice sits on bed below sea level, including deep basins below -1000 m.

bed = ds["bedrock_topography_bed"]
thickness = ds["bedrock_topography_thickness"]
err = ds["bedrock_topography_errbed"]
marine = (mask == 2) & (bed < 0)

marine_summary = xr.Dataset({
    "grounded_ice_bed_below_0m_km2": marine.sum() * km2,
    "grounded_ice_bed_below_minus_1000m_km2": (marine & (bed < -1000)).sum() * km2,
    "mean_bed_elevation_m": bed.where(marine).mean(),
    "mean_ice_thickness_m": thickness.where(marine).mean(),
    "mean_bed_error_m": err.where(marine).mean(),
}).compute()

marine_summary.to_array().to_series().rename("value").to_frame()

Example 2: floating-ice basal melt and a 5 km grounding-zone buffer¶

from scipy.ndimage import binary_dilation

melt = ds["ice_shelf_basal_melt_rate"]
floating = (mask == 3) & melt.notnull()

gl5 = da.map_overlap(
    binary_dilation,
    ds["groundlines_mask"].data.astype(bool),
    depth=5,
    boundary=False,
    dtype=bool,
    structure=np.ones((11, 11), bool),
)
grounding_zone = xr.DataArray(gl5, coords=mask.coords, dims=mask.dims) & floating

melt_summary = xr.Dataset({
    "floating_ice_with_melt_data_km2": floating.sum() * km2,
    "floating_ice_mean_melt_m": melt.where(floating).mean(),
    "floating_ice_net_melt_km3_per_dataset_unit": melt.where(floating).sum() * km2 / 1000,
    "grounding_zone_5km_area_km2": grounding_zone.sum() * km2,
    "grounding_zone_5km_mean_melt_m": melt.where(grounding_zone).mean(),
}).compute()

melt_summary.to_array().to_series().rename("value").to_frame()

client.close()