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
from iconspy.demo_helpers import setup_plot_area_quickstart, setup_plot_area_mar
import cmocean.cm as cmo

Load and prepare the example data

We will load from netcdf files a 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.

The call to convert_tgrid_data can be slow on large grids. For R2B4 (~100 km) grids the conversion is near instantaneous; however for an R2B9 (~5 km) grid this can take around 30 seconds.

[2]:

from iconspy.tests.conftest import get_ds_tgrid_lr, get_ds_fxgrid_lr

[3]:

ds_tgrid = get_ds_tgrid_lr()
ds_fx = get_ds_fxgrid_lr()

# 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",
    }
)

Overview

When constructing a region or a section a typical workflow will be to define a set of TargetStations which we want to connect. We convert these into ModelStations which lie on the model grid, before joining them together to form Sections or Regions.

Target and Model stations

The target station

When creating sections or regions, we typically start with a set of points we want to connect. Such points are represented in ispy as a TargetStation. A TargetStation object is initialised with a name, the target points longitude followed by its latitude, and finally a boolean stating whether we want the point to be on the land boundary or in the ocean.

[4]:

ispy.TargetStation?

Init signature: ispy.TargetStation(name, lon, lat, boundary=True)
Docstring:
Represents a target station

Parameters
----------
name : str
    Name of the target station
lon : float
    Longitude of the target station
lat : float
    Latitude of the target station
boundary : bool, optional
    Whether the station should be sited on a boundary point or not, by default True

Attributes
----------
name : str
    Name of the target station
target_lon : float
    Longitude of the target station
target_lat : float
    Latitude of the target station
boundary : bool
    Whether the station should be sited on a boundary point or not

Example
-------
>>> import iconspy as ispy
>>> target_station = ispy.TargetStation("My Station", 10.0, 20.0)
>>> target_station.plot()
File:           /work/mh0256/m301014/iconspy/iconspy/core.py
Type:           type
Subclasses:

[5]:

target_ne_corner = ispy.TargetStation("NE Corner", -83.72, -0.928, boundary=True)
target_nw_corner = ispy.TargetStation("NW Corner", -95.794, -4.976, boundary=False)
target_sw_corner = ispy.TargetStation("SW Corner", -92.592, -23.219, boundary=False)
target_se_corner = ispy.TargetStation("SE Corner", -70.285, -18.491, boundary=True)

fig, ax = setup_plot_area_quickstart(ds_IsD)

target_ne_corner.plot(ax=ax, gridlines=False)
target_se_corner.plot(ax=ax, gridlines=False)
target_sw_corner.plot(ax=ax, gridlines=False)
target_nw_corner.plot(ax=ax, gridlines=False)

The model station

A TargetStation will not typically lie exactly on a model grid point. We can convert these stations to model stations as such:

[6]:

# Convert the target stations to model stations
model_ne_corner = target_ne_corner.to_model_station(ds_IsD)
model_nw_corner = target_nw_corner.to_model_station(ds_IsD)
model_sw_corner = target_sw_corner.to_model_station(ds_IsD)
model_se_corner = target_se_corner.to_model_station(ds_IsD)

fig, ax = setup_plot_area_quickstart(ds_IsD)

# Plot the model stations (circles)
model_ne_corner.plot(ax=ax, gridlines=False)
model_nw_corner.plot(ax=ax, gridlines=False)
model_sw_corner.plot(ax=ax, gridlines=False)
model_se_corner.plot(ax=ax, gridlines=False)

# Compare with the target station (crosses)
target_ne_corner.plot(ax=ax, gridlines=False)
target_nw_corner.plot(ax=ax, gridlines=False)
target_sw_corner.plot(ax=ax, gridlines=False)
target_se_corner.plot(ax=ax, gridlines=False)

fig.legend()

[6]:

<matplotlib.legend.Legend at 0x7ffb44713e00>

The ModelStation uses a “ball tree” to find the nearest vertex on the model grid to it. If we indicate we want a boundary vertex in our TargetStation, this will be a point on land. For example, our North East Corner TargetStation sits over the ocean (blue cross); however, we have specified we want a boundary vertex so the ModelStation we identify is at the coast (blue circle).

Its worth noting that the ball tree understand the Earth’s sphericity, so you don’t need to worry about using longitudes of [-180, 180) or [0, 360). You will always get the closest “valid” point.

If you want to see the index of the vertex of the ModelStation you can do so as follows:

[7]:

model_ne_corner.vertex

[7]:

Reconstruct sections from a region

Having produced a region we can also extract the sections that make it up.

Sometimes when we form a region from sections, the raw sections cover some of the same edges. When we construct a region these branches are removed. By reconstructing our sections from the region, these branches will also be removed.

[29]:

# Produce a dictionary of reconstructed sections from the region
reconstructed_sections = example_region.extract_sections_from_region(ds_IsD)
print(reconstructed_sections)

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

OrderedDict({'Southern Edge (great circle)': ReconstructedSection(Southern Edge (great circle), SE Corner, SW Corner, great circle), 'Southern Edge (shortest)': ReconstructedSection(Southern Edge (shortest), SW Corner, SE Corner, shortest)})

Note how the reconstructed section is two edges shorter due to the removal of the edges in common with the other section

[28]:

len(reconstructed_sections["Southern Edge (shortest)"].edge_path), len(southern_edge_shortest.edge_path)

[28]:

(9, 11)