Define figure projection

Here we should the most minimal use of PolarToolkit. This example just creates a projection in EPSG:3031, based on a region and figure height (or width). The rest of the example uses standard PyGMT calls.

Import the packages

%%capture
%load_ext autoreload
%autoreload 2

# set Python's logging level to get information
import logging

import pygmt

from polartoolkit import fetch, maps, utils

logging.getLogger().setLevel(logging.INFO)

Define a region for the plot

# Options:

# 1) use the full extent of the grid file
# region = utils.get_grid_info(bed)[1]

# 2) use a preset region (polartoolkitgions())
# region = regions.antarctic_peninsula

# 3) define your own region, in meters e, w, n, s in EPSG:3031

region = (-2700e3, -2000e3, 1000e3, 2000e3)

print(f"region: {region}")

region: (-2700000.0, -2000000.0, 1000000.0, 2000000.0)

Fetch the data to plot

bed = fetch.bedmachine(layer="bed", region=region)
bed

INFO:root:returning grid with new region and/or registration, same spacing

<xarray.DataArray 'z' (y: 2001, x: 1401)>
array([[-3539.7139, -3541.865 , -3542.9233, ...,  1119.9469,  1091.0481,
         1066.7369],
       [-3545.0647, -3548.93  , -3551.0552, ...,  1101.8818,  1076.8281,
         1048.8257],
       [-3550.226 , -3555.766 , -3558.5515, ...,  1084.1982,  1062.6016,
         1036.5276],
       ...,
       [-3507.9219, -3491.199 , -3491.9875, ..., -4500.003 , -4489.7275,
        -4501.537 ],
       [-3520.2815, -3508.0112, -3504.5913, ..., -4513.788 , -4514.832 ,
        -4520.433 ],
       [-3527.554 , -3512.358 , -3506.0007, ..., -4521.106 , -4525.564 ,
        -4525.5522]], dtype=float32)
Coordinates:
  * x        (x) float64 -2.7e+06 -2.7e+06 -2.699e+06 ... -2e+06 -2e+06
  * y        (y) float64 1e+06 1e+06 1.001e+06 ... 1.999e+06 2e+06 2e+06
Attributes:
    long_name:     z
    actual_range:  [-4550.37695312  2708.81054688]

xarray.DataArray

'z'

• y: 2001
• x: 1401

• -3.54e+03 -3.542e+03 -3.543e+03 ... -4.521e+03 -4.526e+03 -4.526e+03

  array([[-3539.7139, -3541.865 , -3542.9233, ...,  1119.9469,  1091.0481,
           1066.7369],
         [-3545.0647, -3548.93  , -3551.0552, ...,  1101.8818,  1076.8281,
           1048.8257],
         [-3550.226 , -3555.766 , -3558.5515, ...,  1084.1982,  1062.6016,
           1036.5276],
         ...,
         [-3507.9219, -3491.199 , -3491.9875, ..., -4500.003 , -4489.7275,
          -4501.537 ],
         [-3520.2815, -3508.0112, -3504.5913, ..., -4513.788 , -4514.832 ,
          -4520.433 ],
         [-3527.554 , -3512.358 , -3506.0007, ..., -4521.106 , -4525.564 ,
          -4525.5522]], dtype=float32)

• Coordinates: (2)
  • x
    (x)
    float64
    -2.7e+06 -2.7e+06 ... -2e+06 -2e+06

    long_name :

    x

    actual_range :

    [-2700000. -2000000.]

    axis :

    X

    array([-2700000., -2699500., -2699000., ..., -2001000., -2000500., -2000000.])

  • y
    (y)
    float64
    1e+06 1e+06 ... 2e+06 2e+06

    long_name :

    y

    actual_range :

    [1000000. 2000000.]

    axis :

    Y

    array([1000000., 1000500., 1001000., ..., 1999000., 1999500., 2000000.])

• Indexes: (2)
  • x
    PandasIndex

    PandasIndex(Index([-2700000.0, -2699500.0, -2699000.0, -2698500.0, -2698000.0, -2697500.0,
           -2697000.0, -2696500.0, -2696000.0, -2695500.0,
           ...
           -2004500.0, -2004000.0, -2003500.0, -2003000.0, -2002500.0, -2002000.0,
           -2001500.0, -2001000.0, -2000500.0, -2000000.0],
          dtype='float64', name='x', length=1401))

  • y
    PandasIndex

    PandasIndex(Index([1000000.0, 1000500.0, 1001000.0, 1001500.0, 1002000.0, 1002500.0,
           1003000.0, 1003500.0, 1004000.0, 1004500.0,
           ...
           1995500.0, 1996000.0, 1996500.0, 1997000.0, 1997500.0, 1998000.0,
           1998500.0, 1999000.0, 1999500.0, 2000000.0],
          dtype='float64', name='y', length=2001))

• Attributes: (2)

  long_name :

  z

  actual_range :

  [-4550.37695312 2708.81054688]

Create a projection from the region and a figure height

proj_xy = utils.set_proj(region, fig_height=15)[0]

Use standard PyGMT commands to plot a figure

fig = pygmt.Figure()

fig.grdimage(
    grid=bed,
    cmap="globe",
    projection=proj_xy,
    region=region,
    frame=True,
    nan_transparent=True,
)

# display the figure
fig.show()

# display colorbar 2/3 as wide as figure

fig.colorbar(
    cmap=True,
    position=f"jBC+w{utils.get_fig_width()*(2/3)}c/.5c+jTC+h+o0c/.6c+e",
    frame="xaf+lBedmachine bed (m)",
)

fig.show()

You can also create the projection by giving a figure width instead of height

proj_xy = utils.set_proj(region, fig_width=15)[0]

# use standard PyGMT commands to plot figure
fig = pygmt.Figure()

# create a custom colarmap
pygmt.makecpt(
    cmap="globe",
    series="-4500/2500/250",  # 250m increments between -4.5 and +2.5 km.
)

fig.grdimage(
    grid=bed,
    cmap=True,
    projection=proj_xy,
    region=region,
    frame=True,
    nan_transparent=True,
)

# display the figure
fig.show()

You can switch between standard PyGMT commands and PolarToolkit commands. Here, on the same figure instance, we’ll add: * a colorbar with a histogram * an inset location map * the Antarctic coastline and groundingline

maps.add_coast(fig, pen="1p,red")

maps.add_inset(fig, inset_width=0.4, inset_box_pen="2p,green")

maps.add_colorbar(
    fig,
    cbar_label="Bedmachine bed elevation",
    cbar_unit="m",
    hist=True,
    cpt_lims=(-4500, 2500),
    grid=bed,
    region=region,
    hist_bin_width=250,  # set this to the cmap interval to match hist bins to cmap bins
    # hist_bin_num=20, # use this instead to set the number of bins
)

fig.show()