Quickstart

The following notebook introduces iconspy’s datatypes and some of their associated methods.

Having followed through this tutorial you will be able to construct and visusalise iconspy sections.

[1]:

from pathlib import Path
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import xarray as xr
import iconspy as ispy
import cmocean.cm as cmo

Load and prepare the example data

We will load from netcdf files as tgrid which describes the model grid, and an fxgrid which contains bathymetry information.

We then have to put them in a format that iconspy can understand. For the tgrid this means calling ispy.convert_tgrid_data, and for the the fxgrid we must make sure the dimensions have the correct names.

[2]:

shared_data_path = Path("/pool/data/ICON/oes/input/r0006/")
R02B04_path = shared_data_path / "icon_grid_0036_R02B04_O"
grid_path = R02B04_path / "R2B4_ocean-grid.nc"
fx_path = R02B04_path / "R2B4L40_fx.nc"

ds_tgrid = xr.open_dataset(grid_path)  # horizontal grid information
ds_fx = xr.open_dataset(fx_path)  # Contains bathymetry etc.

# Put datasets into the iconspy format
ds_IsD = ispy.convert_tgrid_data(ds_tgrid)

ds_fx = ds_fx.rename(
    {
        "ncells": "cell",
        "ncells_2": "edge",
        "ncells_3": "vertex",
    }
)

Define the stations we want to connect

To define a section we need two or more points to connect. Typically we have a rough idea of where we want these points to be, and we would like to choose the model grid points closest to these coordinates. Sometimes we have a desire for one of the points to be on land.

In the below we will define a section near Fram Strait, stretching from a point on Greenland to a point offshore.

[3]:

def setup_plot_area(ds_IsD):
    Slat, Nlat = -20, 20
    Wlon, Elon = -60, 0

    edges_in_region = ds_IsD["edge"].where(
        (ds_IsD["elon"] > Wlon) * (ds_IsD["elon"] < Elon) * (ds_IsD["elat"] > Slat) * (ds_IsD["elat"] < Nlat)
        , drop=True).astype("int32")


    lons = ds_IsD["vlon"].sel(vertex=ds_IsD["edge_vertices"])
    lats = ds_IsD["vlat"].sel(vertex=ds_IsD["edge_vertices"])

    fig, ax = plt.subplots(subplot_kw={"projection": ccrs.PlateCarree()})

    for edge in edges_in_region:
        ax.plot(
            lons.isel(edge=edge),
            lats.isel(edge=edge),
            color="black",
            alpha=0.2,
        )

    ax.scatter(
        ds_IsD["vlon"],
        ds_IsD["vlat"],
        s=4,
        transform=ccrs.PlateCarree(),
        color="tab:green"
    )

    ax.set_xlim(Wlon, Elon)
    ax.set_ylim(Slat, Nlat)
    ax.set_aspect("equal")

    ax.coastlines()
    ax.grid(False)
    return fig, ax

[4]:

# Put in the approximate coordinates of the points we want
# Specify that we want the George Town point to be on the boundary
target_GeorgeTown = ispy.TargetStation("George Town", -59, 6, boundary=True)
target_offshore = ispy.TargetStation("Offshore", -35, 15, boundary=False)
target_Natal = ispy.TargetStation("Natal", -35, -5, boundary=True)

# Visualise the points
fig, ax = setup_plot_area(ds_IsD)
target_GeorgeTown.plot(ax=ax, gridlines=False)
target_offshore.plot(ax=ax, gridlines=False)
target_Natal.plot(ax=ax, gridlines=False)
gl = ax.gridlines()
gl.xlines, gl.ylines = False, False
ax.grid(visible=False)
ax.legend(loc=2)

[4]:

<matplotlib.legend.Legend at 0x7f3b77288250>

The map shows the approximate location of the points we want to join. We can now find the model grid points nearest to them and plot them.

[5]:

# Convert the target stations to model stations
model_GeorgeTown = target_GeorgeTown.to_model_station(ds_IsD)
model_offshore = target_offshore.to_model_station(ds_IsD)
model_Natal = target_Natal.to_model_station(ds_IsD)

fig, ax = setup_plot_area(ds_IsD)

# Plot the model stations (circles)
model_GeorgeTown.plot(ax=ax, gridlines=False)
model_offshore.plot(ax=ax, gridlines=False)
model_Natal.plot(ax=ax, gridlines=False)

# Compare with the target station (crosses)
target_GeorgeTown.plot(ax=ax, gridlines=False)
target_offshore.plot(ax=ax, gridlines=False)
target_Natal.plot(ax=ax, gridlines=False)

fig.legend()

[5]:

<matplotlib.legend.Legend at 0x7f3b771272d0>

Reconstruct sections from a region

Having produced a region we can also extract the sections that make it up. This can be useful as ordinary sections don’t have edge orientations.

[13]:

reconstructed_sections = example_region.extract_sections_from_region(ds_IsD)

[14]:

fig, ax = setup_plot_area(ds_IsD)
for section in reconstructed_sections:
    reconstructed_sections[section].plot(ax=ax, gridlines=False)

In the above we can see a funny straight line in the blue and green section. This is a bug that will be adressed in a future version. It arises from the repeating of the first vertex along the section at the end of the section. The edge paths remain unaffected though.

Quickstart

Load and prepare the example data

Define the stations we want to connect

Connect stations to form sections

Connect sections to form a region

Reconstruct sections from a region