[1]:

%load_ext autoreload
%autoreload 2

Plot data on a map

poligrain shall make your live easier when plotting data as points, lines and grids. In particular plotting lines from CML paths using a colormap will be easy now.

Since we enforce a fixed naming convention of site cooridnates in the xarray.Datasets that we use, we do not have to fiddle with these during plotting. We can just take a Dataset or a DataArray as input for our plotting functions.

[2]:

import cartopy
import matplotlib.pyplot as plt

import poligrain as plg

Load OpenMRG example data

This is a small subset of the full OpenMRG dataset from https://zenodo.org/records/7107689 for which CML data was processed and all data are provided at a 5-minute resolution.

[3]:

(
    ds_rad,
    ds_cmls,
    ds_gauges_municp,
    ds_gauge_smhi,
) = plg.example_data.load_openmrg(data_dir="example_data", subset="5min_2h")

File already exists at example_data/openmrg_cml_5min_2h.nc
Not downloading!
File already exists at example_data/openmrg_rad_5min_2h.nc
Not downloading!
File already exists at example_data/openmrg_municp_gauge_5min_2h.nc
Not downloading!
File already exists at example_data/openmrg_smhi_gauge_5min_2h.nc
Not downloading!

Plot CML paths

If we have CML data in a xarray.Dataset with our naming convention we can simply pass the Dataset to the plotting function to plot all CML paths onto a map. The lines are draw in one go using LineCollection of matplotlib, and not using a for-loop to iterate over the CMLs, so the plotting function is very fast also for large CML networks.

[4]:

plg.plot_map.plot_lines(ds_cmls, use_lon_lat=True);
../_images/notebooks_Plotting_8_0.png

We can plot directly from the xarray.Dataset using the custom Accessor available via .plg

[5]:

fig, ax = plt.subplots(figsize=(4, 4))
ds_cmls.plg.plot_cmls(line_width=0.5, ax=ax, line_color="k", use_lon_lat=True);
../_images/notebooks_Plotting_10_0.png

If we use an xarray.DataArray, as a subset of the CML Dataset that we did load above, we can easily color the lines based on the values of each CML.

Note that there must not be more than one value per CML, i.e. you have to subset to one sublink_id (if you have two, which is the common case) and to one time step.

Below, this is shown for the rain rate per CML from a processed example data. Here only one sublink is provide, hence, we do not have to subset.

[6]:

da_R = ds_cmls.R.isel(time=11)

fig, ax = plt.subplots(figsize=(5, 4))
lines = da_R.plg.plot_cmls(pad_width=1, vmin=0, vmax=1, ax=ax, use_lon_lat=True)
plt.colorbar(lines, label="rainfall rate in mm/h");
../_images/notebooks_Plotting_12_0.png

This way, we can also easily add color based on CML frequency.

[7]:

fig, ax = plt.subplots(figsize=(5, 4))
da_f_GHz = ds_cmls.frequency / 1e3
lines = da_f_GHz.plg.plot_cmls(vmin=5, vmax=40, cmap="viridis", ax=ax, use_lon_lat=True)
plt.colorbar(lines, label="Frequency in GHz");
../_images/notebooks_Plotting_14_0.png

Add a map background

[8]:

lines = plg.plot_map.plot_lines(
    ds_cmls.frequency / 1e3,
    vmin=5,
    vmax=40,
    cmap="magma",
    use_lon_lat=True,
    background_map="OSM",
)
plt.colorbar(lines, label="Frequency in GHz");
../_images/notebooks_Plotting_16_0.png 
[9]:

lines = plg.plot_map.plot_lines(
    ds_cmls.frequency / 1e3,
    vmin=5,
    vmax=40,
    cmap="magma",
    use_lon_lat=True,
    background_map="NE",
    projection=cartopy.crs.EuroPP(),
)

plt.colorbar(lines, label="Frequency in GHz")

# plt.gca().set_extent([5, 20, 52, 62])
plt.gca().gridlines()

[9]:

<cartopy.mpl.gridliner.Gridliner at 0x127ea80b0>
../_images/notebooks_Plotting_17_1.png

Plot SML paths

to be added…

Plot gauges

[10]:

fig, ax = plt.subplots(figsize=(5, 4))
plg.plot_map.plot_plg(
    da_gauges=ds_gauges_municp.rainfall_amount.sum(dim="time"),
    use_lon_lat=True,
    vmin=0,
    vmax=8,
    cmap="turbo",
    marker_size=100,
    ax=ax,
    colorbar_label="gauge rainfall sum (mm)",
);
../_images/notebooks_Plotting_21_0.png

Plot radar data

to be added…

Plot all sensors at the same time

poligrain provides a function that plots grid, line and point data with a shared color map onto a lon-lat or x-y (in projected coordinates) grid.

The goal of this function is to

  • make it easy to do these standard plots used for comparison

  • provide good defaults parameters so the plots are immediately helpful when visually inspecting results

[11]:

# we still need to rename variables here, until we update the example dataset...
ds_rad = ds_rad.rename({"longitudes": "lon", "latitudes": "lat"})
ds_rad = ds_rad.set_coords(["lon", "lat"])

[12]:

plg.plot_map.plot_plg(
    da_grid=ds_rad.rainfall_amount.sum(dim="time").where(
        ds_rad.rainfall_amount.sum(dim="time")
    ),
    da_gauges=ds_gauges_municp.rainfall_amount.sum(dim="time"),
    use_lon_lat=True,
    da_cmls=ds_cmls.R.sum(dim="time"),
    vmin=0,
    vmax=10,
    colorbar_label="rainfall sum (mm)",
);
../_images/notebooks_Plotting_27_0.png

We can also just plot one or tow of the datasets, e.g. only radar and gauge or only radar and CML

[13]:

plg.plot_map.plot_plg(
    da_grid=ds_rad.rainfall_amount.sum(dim="time"),
    da_gauges=ds_gauges_municp.rainfall_amount.sum(dim="time"),
    use_lon_lat=True,
    vmin=0,
    vmax=10,
    colorbar_label="rainfall sum (mm)",
);
../_images/notebooks_Plotting_29_0.png 
[14]:

plg.plot_map.plot_plg(
    da_grid=ds_rad.rainfall_amount.sum(dim="time"),
    da_cmls=ds_cmls.R.sum(dim="time"),
    use_lon_lat=True,
    vmin=0,
    vmax=10,
    colorbar_label="rainfall sum (mm)",
);
../_images/notebooks_Plotting_30_0.png

We can also adjust the styling of the plotting of the points and lines separately by passing kwargs as a dictionary which is then passed to scatter from matplotlib (for points) and to plot_lines of poligrain (for lines). To see all possible kwargs, check the documentation of these plotting functions. Here we just shown an example.

[15]:

plg.plot_map.plot_plg(
    da_grid=ds_rad.rainfall_amount.sum(dim="time"),
    da_gauges=ds_gauges_municp.rainfall_amount.sum(dim="time"),
    use_lon_lat=True,
    da_cmls=ds_cmls.R.sum(dim="time"),
    vmin=0,
    vmax=5,
    cmap="Grays",
    colorbar_label="radar rainfall sum (mm)",
    kwargs_cmls_plot={
        "edge_color": "b",
        "edge_width": 1,
        "vmax": 10,
        "cmap": "magma_r",
    },
    kwargs_gauges_plot={
        "marker": "*",
        "s": 100,
        "cmap": "Greens",
        "vmin": 4,
        "vmax": 6,
    },
);
../_images/notebooks_Plotting_32_0.png

Or you can just plot the lines and points with a single color to get an overview of the available sensors in relation to the radar grid.

We do this by passing CMLs and gauges as xarray.Dataset instead of xarray.DataArray signalizing that we are not interested in specific values, but in the general structure of the all sensors in the respective Dataset.

[16]:

plg.plot_map.plot_plg(
    da_grid=ds_rad.rainfall_amount.sum(dim="time"),
    da_gauges=ds_gauges_municp,
    da_cmls=ds_cmls,
    use_lon_lat=True,
    line_color="k",
    point_color="r",
    vmin=0,
    vmax=5,
);
../_images/notebooks_Plotting_34_0.png

Or just plot gauges and CMLs without radar data.

[17]:

plg.plot_map.plot_plg(
    da_gauges=ds_gauges_municp,
    da_cmls=ds_cmls,
    use_lon_lat=True,
    line_color="k",
    point_color="r",
    vmin=0,
    vmax=5,
);
../_images/notebooks_Plotting_36_0.png

The same plot can also be done with projected coordinates. But of course, first the projected coordinates have to be there.

So, let’s project the data and add the coordinates according to the poligrain naming convention.

[18]:

# UTM32N: https://epsg.io/32632
ref_str = "EPSG:32632"

(
    ds_gauges_municp.coords["x"],
    ds_gauges_municp.coords["y"],
) = plg.spatial.project_point_coordinates(
    ds_gauges_municp.lon, ds_gauges_municp.lat, ref_str
)

(
    ds_cmls.coords["site_0_x"],
    ds_cmls.coords["site_0_y"],
) = plg.spatial.project_point_coordinates(
    ds_cmls.site_0_lon, ds_cmls.site_0_lat, ref_str
)
(
    ds_cmls.coords["site_1_x"],
    ds_cmls.coords["site_1_y"],
) = plg.spatial.project_point_coordinates(
    ds_cmls.site_1_lon, ds_cmls.site_1_lat, ref_str
)

(
    ds_rad.coords["x_grid"],
    ds_rad.coords["y_grid"],
) = plg.spatial.project_point_coordinates(ds_rad.lon, ds_rad.lat, ref_str)

Then do the same plot as above but in projected coordinates.

[19]:

plg.plot_map.plot_plg(
    da_grid=ds_rad.rainfall_amount.sum(dim="time"),
    da_gauges=ds_gauges_municp.rainfall_amount.sum(dim="time"),
    use_lon_lat=False,
    da_cmls=ds_cmls.R.sum(dim="time"),
    vmin=0,
    vmax=5,
    cmap="Grays",
    colorbar_label="radar rainfall sum (mm)",
    kwargs_cmls_plot={
        "edge_color": "b",
        "edge_width": 1,
        "vmax": 10,
        "cmap": "magma_r",
    },
    kwargs_gauges_plot={
        "marker": "*",
        "s": 100,
        "cmap": "Greens",
        "vmin": 4,
        "vmax": 6,
    },
);
../_images/notebooks_Plotting_40_0.png

Add background map

[20]:

plg.plot_map.plot_plg(
    da_grid=ds_rad.rainfall_amount.sum(dim="time").where(
        ds_rad.rainfall_amount.sum(dim="time")
    ),
    da_gauges=ds_gauges_municp.rainfall_amount.sum(dim="time"),
    use_lon_lat=True,
    da_cmls=ds_cmls.R.sum(dim="time"),
    vmin=0,
    vmax=10,
    colorbar_label="rainfall sum (mm)",
    alpha=0.5,
    background_map="NE",
    projection=cartopy.crs.EuroPP(),
);
../_images/notebooks_Plotting_42_0.png 
[21]:

plg.plot_map.plot_plg(
    da_gauges=ds_gauges_municp.rainfall_amount.sum(dim="time"),
    use_lon_lat=True,
    da_cmls=ds_cmls.R.sum(dim="time"),
    vmin=0,
    vmax=10,
    colorbar_label="rainfall sum (mm)",
    alpha=0.5,
    background_map="NE",
    projection=cartopy.crs.EuroPP(),
    # extent=[-10, 30, 40, 60],
);
# plt.gca().set_extent([-10, 30, 40, 60])
../_images/notebooks_Plotting_43_0.png 
[22]:

plg.plot_map.set_up_axes(
    background_map="NE",
    extent=[-10, 30, 40, 60],
    projection=cartopy.crs.EuroPP(),
)

[22]:

<GeoAxes: >
../_images/notebooks_Plotting_44_1.png 
[ ]: