# -*- coding: utf-8 -*-

# ---

# jupyter:

# jupytext:

# formats: ipynb,py:hydrogen

# text_representation:

# extension: .py

# format_name: hydrogen

# format_version: '1.3'

# jupytext_version: 1.13.4

# kernelspec:

# display_name: Python 3 (ipykernel)

# language: python

# name: python3

# ---

# %% [markdown]

# # **ICESat-2 detected active subglacial lakes**

#

# Making pretty maps of active subglacial lakes in Antarctica,

# a companion jupyter notebook to https://github.com/weiji14/deepicedrain.

# Uses [PyGMT](https://www.pygmt.org) for illustration,

# heavily inspired by https://github.com/mrsiegfried/Venturelli2020-GRL

# %%

import os

import re

import geopandas as gpd

import numpy as np

import pandas as pd

import pooch

import pygmt

# %% [markdown]

# # Get data files

#

# Data for background basemaps:

#

# - Haran, T. M., Bohlander, J., Scambos, T. A., Painter, T. H., & Fahnestock, M. A. (2014). MODIS Mosaic of Antarctica 2008-2009 (MOA2009) Image Map. U.S. Antarctic Program Data Center (USAP-DC), via National Snow and Ice Data Center (NSIDC). https://doi.org/10.7265/N5KP8037

# - Depoorter, M. A., Bamber, J. L., Griggs, J. A., Lenaerts, J. T. M., Ligtenberg, S. R. M., van den Broeke, M. R., & Moholdt, G. (2013). Antarctic masks (ice-shelves, ice-sheet, and islands), link to shape file (p. 15.8 MBytes) [Application/zip]. PANGAEA - Data Publisher for Earth & Environmental Science. https://doi.org/10.1863/PANGAEA.819147

# - Mouginot, J., Rignot, E., & Scheuchl, B. (2019). MEaSUREs Phase Map of Antarctic Ice Velocity, Version 1 [Data set]. NASA National Snow and Ice Data Center DAAC. https://doi.org/10.5067/PZ3NJ5RXRH10

# %%

# Using Quantarctica3 from https://www.npolar.no/quantarctica/

datafold: str = os.getenv("DATAHOME") or os.path.abspath("Quantarctica3")

os.makedirs(name=datafold, exist_ok=True)

# %%

# MODIS Mosaic of Antarctica

moa_no_nan: str = pooch.retrieve(

)

moa = f"{datafold}/SatelliteImagery/MODIS/moa750_2009_hp1_v01.1.tif"

try:

assert os.path.exists(path=moa)

except AssertionError:

with pygmt.clib.Session() as lib:

# !gmt grdmath $moa_no_nan 0 NAN = $moa

lib.call_module(module="grdmath", args=f"{moa_no_nan} 0 NAN = {moa}")

# %%

# Scripps Grounding Line

shapefiles: list = pooch.retrieve(

)

groundingline: str = [file for file in shapefiles if file.endswith(".shp")][0]

gdf_groundingline: gpd.GeoDataFrame = gpd.read_file(filename=groundingline)

# %%

# MEaSUREs Phase Map of Antarctic Ice Velocity

vel_file = f"{datafold}/Glaciology/MEaSUREs_PhaseBased_Velocity/antarctic_ice_vel_phase_map_v01"

vel: str = f"{vel_file}-vmag.nc"

try:

assert os.path.exists(vel)

except AssertionError:

# Note, download require a .netrc file containing 'machine urs.earthdata.nasa.gov login <uid> password <password>'

# see https://nsidc.org/support/how/how-do-i-programmatically-access-data-spatial-temporal

try:

# Run processing of velocity magnitude grid locally

vel_file: str = pooch.retrieve(

url="https://n5eil01u.ecs.nsidc.org/MEASURES/NSIDC-0754.001/1996.01.01/antarctic_ice_vel_phase_map_v01.nc",

known_hash="fa0957618b8bd98099f4a419d7dc0e3a2c562d89e9791b4d0ed55e6017f52416",

fname="antarctic_ice_vel_phase_map_v01.nc",

path=f"{datafold}/Glaciology/MEaSUREs_PhaseBased_Velocity",

)

with pygmt.clib.Session() as lib:

#! gmt grdmath ${vel_file}.nc?VX 2 POW ${vel_file}.nc?VY 2 POW POW 0.5 = ${vel_file}-vmag.nc

lib.call_module(

module="grdmath",

args=f"{vel_file}.nc?VX 2 POW {vel_file}.nc?VY 2 POW POW 0.5 = {vel}",

)

except:

# Just download pre-processed velocity magnitude grid from GitHub

vel: str = pooch.retrieve(

url="https://github.com/weiji14/nzasc2021/releases/download/v0.0.0/antarctic_ice_vel_phase_map_v01-vmag.nc",

known_hash="ed6393275d8d8475c2162a838d6b9220cd529d28a2b5d674a6bf6dbda4971049",

fname="antarctic_ice_vel_phase_map_v01-vmag.nc",

path=f"{datafold}/Glaciology/MEaSUREs_PhaseBased_Velocity",

)

# %%

# DeepIceDrain active subglacial lake outlines

lakes = "https://raw.githubusercontent.com/weiji14/deepicedrain/0cec859288b2add98a42b095955c8ec644dd616f/antarctic_subglacial_lakes_3031.geojson"

gdf_lakes: gpd.GeoDataFrame = gpd.read_file(filename=lakes)

# %% [markdown]

# ### Make color maps

#

# For MOA and ice velocity and vertical elevation trend (dhdt).

# Scientific color maps are from http://doi.org/10.5281/zenodo.1243862

# %%

pygmt.makecpt(

series=[15000, 17000, 1],

cmap="grayC",

continuous=True,

output="cmap_moa.cpt",

reverse=True,

)

with pygmt.config(COLOR_FOREGROUND="240/249/33", COLOR_BACKGROUND="13/8/135"):

pygmt.makecpt(series=[0, 800, 1], cmap="batlow", output="cmap_vel.cpt")

pygmt.makecpt(

cmap="berlin",

series=[-3, 3, 1],

reverse=True,

continuous=True,

output="cmap_dhdt.cpt",

)

# %%

# %% [markdown]

# # Figure of Siple Coast active subglacial lakes

# %%

# We're making this a specific height

figheight = 115 # in mm

# Region in PS71 for main part of the figure

sip_xl, sip_xh, sip_yl, sip_yh = sipreg = [-800_000, 25_000, -1_000_000, -400_000]

# Calculate the figure width and map scale

figwidth = figheight * (sip_xh - sip_xl) / (sip_yh - sip_yl)

sipratio = (sip_yh - sip_yl) / (figheight / 1000)

# Make a GMT region string and projection strings in both PS71 and Lon/Lat

sipreg = [sip_xl, sip_xh, sip_yl, sip_yh]

sipproj = f"x1:{sipratio}"

sipproj_ll = f"s0/-90/-71/1:{sipratio}"

# %%

# Initialize figure and plot MOA as the basemap with ticks every 200 km in xy directions

fig = pygmt.Figure()

with pygmt.config(MAP_FRAME_TYPE="inside"):

fig.basemap(

region=sipreg, projection=sipproj, frame=["nwse", "xf200000", "yf200000"]

)

fig.grdimage(grid=moa, cmap="cmap_moa.cpt", nan_transparent=True)

# Plot graticules overtop, every 2° latitude and 15° longitude

with pygmt.config(

MAP_ANNOT_OFFSET_PRIMARY="-2p",

MAP_FRAME_TYPE="inside",

MAP_ANNOT_OBLIQUE=0,

FONT_ANNOT_PRIMARY="8p,grey",

MAP_GRID_PEN_PRIMARY="grey",

MAP_TICK_LENGTH_PRIMARY="-10p",

MAP_TICK_PEN_PRIMARY="thinnest,grey",

FORMAT_GEO_MAP="dddF",

MAP_POLAR_CAP="88/90", # less longitude graticules at >88°S

):

fig.basemap(

projection=sipproj_ll, region=sipreg, frame=["NSWE", "xa15g15", "ya2g2"]

)

# Plot the grounding line in white

fig.plot(data=gdf_groundingline, region=sipreg, projection=sipproj, pen="0.15p,white")

# %%

# Overlay ice velocity with 70% transparency

fig.grdimage(grid=vel, cmap="cmap_vel.cpt", transparency=70, nan_transparent=True)

# Overlay dhdt with 30% transparency

# pygmt.makecpt(cmap="berlin", series=[-1.0, 1.0, 0.25], continuous=True, reverse=True)

# fig.grdimage(

# grid="ATLXI/ds_grid_dhdt_siple_coast.nc",

# cmap=True,

# # cmap="cmap_dhdt.cpt",

# transparency=30,

# nan_transparent=True,

# )

fig.show()

# %%

# Plot lakes in PS71 as blobs (red for draining, blue for filling)

fig.plot(

data=gdf_lakes,

pen="thinnest,yellow,-",

cmap="cmap_dhdt.cpt",

color=" z",

close=True,

# transparency=30,

aspatial="Z=inner_dhdt",

)

# %%

# Siple Coast placename labels

gdf = gpd.read_file("antarctic_subglacial_lakes_3031.gmt")

with open("place_labels_siple_coast.tsv", mode="w") as file:

# Ice Streams A to E

font = "7p,Helvetica-Narrow-Oblique,white"

print(f"-320000\t-440000\t-65\t{font}\tCM\tMercer Ice Stream", file=file)

print(f"-385000\t-555000\t-5\t{font}\tCM\tWhillans Ice Stream", file=file)

print(f"-470000\t-450000\t-55\t{font}\tCM\tVan der Veen Ice Stream", file=file)

print(f"-550000\t-625000\t-40\t{font}\tCM\tKamb Ice Stream", file=file)

print(f"-700000\t-700000\t-45\t{font}\tCM\tBindschadler Ice Stream", file=file)

print(f"-700000\t-850000\t-37\t{font}\tCM\tMacAyeal Ice Stream", file=file)

# Ice ridges, rises and domes

print(f"-370000\t-480000\t-30\t{font}\tCM\tConway Ice Ridge", file=file)

print(f"-400000\t-600000\t0\t{font}\tCM\tEngelhardt Ice Ridge", file=file)

print(f"-100000\t-750000\t-45\t{font}\tCM\tCrary Ice Rise", file=file)

print(f"-450000\t-780000\t-35\t{font}\tCM\tSiple Dome", file=file)

print(f"-650000\t-800000\t-25\t{font}\tCM\tShabtaie Ice Ridge", file=file)

print(f"-650000\t-950000\t-15\t{font}\tCM\tHarrison Ice Ridge", file=file)

# Abbreviated lakes

abbrev_dict: dict = {

"Subglacial": "S",

"Lake": "L",

"Conway": "C",

"Engelhardt": "E",

"Kamb": "K",

"MacAyeal": "Mac",

"Mercer": "M",

# "Recovery", "R",

# "Slessor","S"

"Whillans": "W",

}

pattern = re.compile(pattern=rf"\b({'|'.join(abbrev_dict.keys())})\b")

for idx, row in gdf.dissolve(by="lake_name", as_index=False).iterrows():

x, y = row.geometry.centroid.xy

label = pattern.sub(

repl=lambda name: abbrev_dict[name.group()], string=row.lake_name

).replace(" ", "")

justify = (

"BL"

if label in ("L12", "L78", "SLW", "W7", "WX", "WXI")

else "TC"

if ("*" in label or label == "SLM")

else "TR"

)

print(f"{x[0]:.0f}\t{y[0]:.0f}\t0\t6p,white\t{justify}\t{label}", file=file)

# %%

# Plot labels for Siple Coast ice streams, active subglacial lakes, etc

fig.text(

textfiles="place_labels_siple_coast.tsv",

angle=True,

font=True,

justify=True,

offset="j0.12c",

# frame=["WsNe", "af10000g50000"],

)

fig.show()

# %%

# Plot the color bar once with a transparent box, then again with no box and no transparency

with pygmt.config(

FONT_ANNOT_PRIMARY="6p,white",

FONT_LABEL="6p,white",

MAP_ANNOT_OFFSET_PRIMARY="2p",

MAP_TICK_PEN_PRIMARY="0.25p,white",

MAP_TICK_LENGTH_PRIMARY="3p",

MAP_FRAME_PEN="0.5p,white",

MAP_LABEL_OFFSET="4p",

):

colorbar_kwargs = dict(

cmap="cmap_dhdt.cpt",

position="jBR jBR w1.6c/0.18c o1.2c/0.3c v e",

frame=["xaf", 'y l"dhdt (m/yr)"'],

)

fig.colorbar(

box=" gblack c-9p/3p",

# box = ' gblack p0.5p,black c3p'

transparency=70,

**colorbar_kwargs,

)

fig.colorbar(**colorbar_kwargs)

# Add a scalebar

fig.basemap(

projection=sipproj_ll, region=sipreg, map_scale="jBR o2.2c/0.3c w50k uk f"

)

fig.show()

# %%

# Make inset overview map of Antarctica

antwidth: int = 3

with fig.inset(position=f"jTR w{antwidth}c"):

# Plot the inset map

fig.basemap(region=vel, projection=f"S0/-90/71/{antwidth}c", frame=" n")

fig.coast(area_thresh=" a", land="white") # ice shelf in white

fig.coast(area_thresh=" ag", land="gray") # grounded ice in gray

fig.plot(

projection=f"X{antwidth}c",

x=[sip_xl, sip_xl, sip_xh, sip_xh, sip_xl],

y=[sip_yl, sip_yh, sip_yh, sip_yl, sip_yl],

pen="1p,black", # map location in black

)

fig.show()

# %%

# Save the figure

fig.savefig(fname="siple_coast_lakes.pdf")

fig.savefig(fname="siple_coast_lakes.png", dpi=1200)

# %%

# %% [markdown]

# # Figure of Antarctic active subglacial lake map

# %%

# We're making this a specific height

figheight = 115 # in mm

# Region in PS71 for main part of the figure

ais_xl, ais_xh, ais_yl, ais_yh = aisreg = [-2700000, 2800000, -2200000, 2300000]

# Calculate the figure width and map scale

figwidth = figheight * (ais_xh - ais_xl) / (ais_yh - ais_yl)

aisratio = (ais_yh - ais_yl) / (figheight / 1000)

# Make a GMT region string and projection strings in both PS71 and Lon/Lat

aisreg = [ais_xl, ais_xh, ais_yl, ais_yh]

aisproj = "x1:" str(aisratio)

aisproj_ll = "s0/-90/-71/1:" str(aisratio)

# %%

# Initialize figure and plot MOA as the base map with ticks every 200 km both directions

fig = pygmt.Figure()

with pygmt.config(MAP_FRAME_TYPE="inside"):

# fig.coast(

# projection=aisproj_ll, region=aisreg, land="lightblue", water="royalblue2"

# )

fig.coast(region=aisreg, projection=aisproj_ll, resolution="c", water=True)

fig.grdimage(

grid="@earth_relief_03m",

projection=aisproj_ll,

region=aisreg,

cmap="fes",

shading=True,

)

fig.coast(Q=True) # end water clip path

fig.basemap(

projection=aisproj, region=aisreg, frame=["nwse", "xf200000", "yf200000"]

)

fig.grdimage(grid=moa, cmap="cmap_moa.cpt", nan_transparent=True)

# Plot graticules overtop, every 10° latitude and 45° longitude

with pygmt.config(

MAP_ANNOT_OFFSET_PRIMARY="-2p",

MAP_FRAME_TYPE="inside",

MAP_ANNOT_OBLIQUE=0,

FONT_ANNOT_PRIMARY="8p,grey",

MAP_GRID_PEN_PRIMARY="grey",

MAP_TICK_LENGTH_PRIMARY="-10p",

MAP_TICK_PEN_PRIMARY="thinnest,grey",

FORMAT_GEO_MAP="dddF",

MAP_POLAR_CAP="88/90", # less longitude graticules at >88°S

):

fig.basemap(

projection=aisproj_ll, region=aisreg, frame=["NSWE", "xa45g45", "ya10g10"]

)

# Plot the grounding line in white

fig.plot(data=gdf_groundingline, region=aisreg, projection=aisproj, pen="0.15p,white")

# Plot bounding box of Siple Coast study area

fig.plot(

x=[sip_xl, sip_xl, sip_xh, sip_xh, sip_xl],

y=[sip_yl, sip_yh, sip_yh, sip_yl, sip_yl],

pen="0.5p,black", # map location in black

transparency=50,

)

# %%

# Overlay ice velocity with 70% transparency

fig.grdimage(grid=vel, cmap="cmap_vel.cpt", transparency=70, nan_transparent=True)

fig.show()

# %%

# Antarctica placename labels

with open("place_labels_antarctica.tsv", mode="w") as file:

font = "10p,Helvetica-Narrow,white"

print(f"-1070000\t-360000\t0\t{font}\tWest", file=file)

print(f"-1070000\t-460000\t0\t{font}\tAntarctica", file=file)

print(f"950000\t600000\t0\t{font}\tEast", file=file)

print(f"950000\t500000\t0\t{font}\tAntarctica", file=file)

font = "6p,Helvetica-Narrow-Oblique,white"

print(f"0\t-1000000\t0\t{font}\tRoss", file=file)

print(f"0\t-1060000\t0\t{font}\tIce Shelf", file=file)

print(f"-1000000\t700000\t0\t{font}\tRonne-Filchner", file=file)

print(f"-1000000\t640000\t0\t{font}\tIce Shelf", file=file)

print(f"2100000\t760000\t0\t{font}\tAmery", file=file)

print(f"2100000\t700000\t0\t{font}\tIce Shelf", file=file)

font = "4p,Helvetica-Narrow-Oblique,white"

print(f"-400000\t1150000\t45\t{font}\tSlessor Glacier", file=file)

print(f"-95000\t850000\t0\t{font}\tRecovery Glacier", file=file)

print(f"-340000\t200000\t-30\t{font}\tFoundation Ice Stream", file=file)

print(f"-810000\t95000\t-60\t{font}\tInstitute Ice Stream", file=file)

# print(f"-1500000\t50000\t20\t{font}\tRutford Ice Stream", file=file)

print(f"-1300000\t-300000\t30\t{font}\tThwaites Glacier", file=file)

print(f"-500000\t-700000\t40\t{font}\tSiple Coast", file=file)

print(f"500000\t-440000\t0\t{font}\tNimrod 2", file=file)

print(f"550000\t-750000\t35\t{font}\tByrd Glacier", file=file)

print(f"700000\t-1500000\t5\t{font}\tDavid Glacier", file=file)

print(f"870000\t-1700000\t0\t{font}\tCook E2", file=file)

# print(f"2200000\t-850000\t-70\t{font}\tTotten Glacier", file=file)

print(f"1400000\t650000\t15\t{font}\tLambert Glacier", file=file)

# %%

# Plot labels for Antarctic ice shelves, ice streams, etc

fig.text(

region=aisreg,

projection=aisproj,

textfiles="place_labels_antarctica.tsv",

angle=True,

font=True,

# frame=["WsNe", "af100000g500000"],

)

# %%

# Plot lakes in PS71 as cyan blobs with 60% transparency

fig.plot(data=gdf_lakes, pen="0.5p,cyan", color="cyan", transparency=60)

fig.show()

# %%

# Add the finishing touches by putting a panel label in the bottom left corner,

# and add a legend with a 70% transparent box behind it (so you can see it

# clearly) that includes:

# - a cyan ellipse the cyan lakes

# - a color bar for the ice surface velocity

barwidth = 0.3 * figwidth / 10 # legend width is 30% of the map width, in cm

# Position string for feeding to the colorbar call. A bit coded, but the pyGMT

# documents are great: https://www.pygmt.org/dev/api/index.html

legend_pos: str = f"jBL jBL w{barwidth}c o0.2c/0.2c"

# Plot the color bar once with a transparent box, then again with no box and no transparency

with pygmt.config(

FONT_ANNOT_PRIMARY="8p,white",

FONT_LABEL="8p,white",

MAP_ANNOT_OFFSET_PRIMARY="2p",

MAP_TICK_PEN_PRIMARY="0.5p,white",

MAP_TICK_LENGTH_PRIMARY="3p",

MAP_FRAME_PEN="0.5p,white",

MAP_LABEL_OFFSET="4p",

):

fig.legend(

spec="legend.txt", position=legend_pos, box=" gblack c1.2p", transparency="70"

)

fig.legend(spec="legend.txt", position=legend_pos)

fig.show()

# %%

# Save the figure

fig.savefig(fname="antarctica_lakes.pdf")

fig.savefig(fname="antarctica_lakes.png", dpi=600)

fig.savefig(fname="antarctica_lakes.jpg", dpi=900)

# %%

# %% [markdown] tags=[]

# ## Name lake clusters

#

# Manually selecting some subglacial lakes by their id(s) to give them a name.

# Export this set to an OGR GMT format that can be plotted easily later.

# %%

# Name the lakes!

import deepicedrain

import pandas as pd

import pygmt

import geopandas as gpd

import os

_ = pygmt.which(fname=lakes, download=True)

lake_catalog = deepicedrain.catalog.subglacial_lakes()

antarctic_lakes = lake_catalog.read()

antarctic_lakes.insert(loc=0, column="lake_name", value=pd.NA)

for _, lake_data in pd.json_normalize(lake_catalog.metadata["lakedict"]).iterrows():

for id in lake_data.ids:

antarctic_lakes.at[id, "lake_name"] = lake_data.lakename

subset_lakes = antarctic_lakes.dropna(subset=["lake_name"]).set_index(

keys=["lake_name"]

)

subset_lakes: gpd.GeoDataFrame = subset_lakes.sort_values(

by=["basin_name", "lake_name"]

)

# %%

# Save subset of lakes

# subset_lakes = gdf.query("basin_name == 'Whillans'")

try:

os.remove(f"antarctic_subglacial_lakes_3031.gmt")

except FileNotFoundError:

pass

subset_lakes.to_file(

filename=f"antarctic_subglacial_lakes_3031.gmt", driver="OGR_GMT", index=True

)