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]:

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

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.

[3]:

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:

[4]:

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:

[9]:

# 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()

[9]:

<matplotlib.legend.Legend at 0x7ffea9846bd0>

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:

[5]:

model_ne_corner.vertex

[5]:

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.