{ "metadata": { }, "nbformat": 4, "nbformat_minor": 5, "cells": [ { "id": "metadata", "cell_type": "markdown", "source": "
\n\n# Pangeo Notebook in Galaxy - Introduction to Xarray\n\nby [Anne Fouilloux](https://training.galaxyproject.org/hall-of-fame/annefou/)\n\nCC-BY licensed content from the [Galaxy Training Network](https://training.galaxyproject.org/)\n\n**Objectives**\n\n- What is Pangeo notebook?\n- How to start Pangeo Notebook in Galaxy?\n- What are the main software components of the Pangeo ecosystem?\n- What is Xarray?\n- How to manipulate Xarrays?\n- How to print metadata information?\n- How to make a selection?\n- How to visualize?\n- How to filter?\n- How to make reduction operations (mean, max, min)?\n- How to resample my data?\n- Where to go next?\n\n**Objectives**\n\n- Learn to get metadata information using Xarray Galaxy Tools\n- Learn to select data\n- Learn to visualize data\n- Learn to filter, make reduction operations (mean, max, min)\n- Learn to resample my data\n- Learn to cite and contribute to Pangeo\n\n**Time Estimation: 1H**\n
\n", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-0", "source": "

Introduction

\n

In this tutorial, we will learn about Xarray, one of the most used Python library from the Pangeo ecosystem.

\n

We will be using data from Copernicus Atmosphere Monitoring Service\nand more precisely PM2.5 (Particle Matter < 2.5 μm) 4 days forecast from December, 22 2021. Parallel data analysis with Pangeo is not covered in this tutorial.

\n
\n
💬 Comment: Remark
\n

This tutorial uses data on a regular latitude-longitude grid. More complex and irregular grids are not discussed in this tutorial. In addition,\nthis tutorial is not meant to cover all the different possibilities offered by Xarrays but shows functionalities we find useful for day to day\nanalysis.

\n
\n
\n
Agenda
\n

In this tutorial, we will cover:

\n
    \n
  1. Introduction
    2. \n
  2. Analysis
      \n
    1. Import Python packages
      2. \n
    \n
    3. \n
\n
\n

Analysis

\n

Import Python packages

\n

Some packages may need to be installed first. For example cmcrameri is missing, so we need to install it by entering the following command in a new cell of your Jupyter Notebook:

\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-1", "source": [ "pip install cmcrameri" ], "cell_type": "code", "execution_count": null, "outputs": [ ], "metadata": { "attributes": { "classes": [ "python" ], "id": "" } } }, { "id": "cell-2", "source": "

Then we need to import all the necessary packages in our Jupyter Notebook.

\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-3", "source": [ "import numpy as np\n", "import xarray as xr\n", "import cartopy.crs as ccrs\n", "import matplotlib.pyplot as plt\n", "import cmcrameri.cm as cmc\n", "import pandas as pd" ], "cell_type": "code", "execution_count": null, "outputs": [ ], "metadata": { "attributes": { "classes": [ "python" ], "id": "" } } }, { "id": "cell-4", "source": "

Open and read metadata

\n

The netCDF dataset can now be opened with Xarray:

\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-5", "source": [ "dset = xr.open_dataset(\"data/CAMS-PM2_5-20211222.netcdf\")" ], "cell_type": "code", "execution_count": null, "outputs": [ ], "metadata": { "attributes": { "classes": [ "python" ], "id": "" } } }, { "id": "cell-6", "source": "

Once opened, we can get metadata using print statement.

\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-7", "source": [ "print(dset)" ], "cell_type": "code", "execution_count": null, "outputs": [ ], "metadata": { "attributes": { "classes": [ "python" ], "id": "" } } }, { "id": "cell-8", "source": "

Below is what you should get if everything goes fine.

\n
\n
🖥 Output
\n
 <xarray.Dataset>\n Dimensions:     (longitude: 700, latitude: 400, level: 1, time: 97)\n Coordinates:\n   * longitude   (longitude) float32 -24.95 -24.85 -24.75 ... 44.75 44.85 44.95\n   * latitude    (latitude) float32 69.95 69.85 69.75 69.65 ... 30.25 30.15 30.05\n   * level       (level) float32 0.0\n   * time        (time) timedelta64[ns] 00:00:00 01:00:00 ... 4 days 00:00:00\nData variables:\n  pm2p5_conc  (time, level, latitude, longitude) float32 0.4202 ... 7.501\nAttributes:\n  title:        PM25 Air Pollutant FORECAST at the Surface\n  institution:  Data produced by Meteo France\n  source:       Data from ENSEMBLE model\n  history:      Model ENSEMBLE FORECAST\n  FORECAST:     Europe, 20211222+[0H_96H]\n  summary:      ENSEMBLE model hourly FORECAST of PM25 concentration at the...\n  project:      MACC-RAQ (http://macc-raq.gmes-atmosphere.eu)\n
\n
\n
\n
💡 Tip: Command not found
\n

If you get an error with the previous command, first check the location of the input file CAMS-PM2_5-20211222.netcdf:\nit needs to be in the same directory as your Jupyter Notebook.

\n
\n

We can identify 4 different sections:

\n
    \n
  1. Dimensions: name of dimensions and corresponding number of elements;
    2. \n
  2. Coordinates: contains coordinate arrays (longitude, latitude, level and time) with their values.
    3. \n
  3. Data variables: contains all the variables available in the dataset. Here, we only have one variable. For each variable, we get information on its shape and values.
    4. \n
  4. Attributes: at this level, we get all the attributes of the dataset.
    5. \n
\n

We can also get metadata information for each coordinate and data variables using “.” followed by the coordinate or data variable name.

\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-9", "source": [ "print(dset.time)" ], "cell_type": "code", "execution_count": null, "outputs": [ ], "metadata": { "attributes": { "classes": [ "python" ], "id": "" } } }, { "id": "cell-10", "source": "
\n
❓ Question: CAM PM2.5 Dataset
\n

What is the name of the variable for Particle matter < 2.5 μm and its physical units?

\n
👁 View solution\n
👁 Solution
\n
    \n
  1. To get metadata information from pm2p5_conc Data variable, we use its variable name and print it. Printing it will only print metadata, not the values.\n
      \n
    • Variable name: mass_concentration_of_pm2p5_ambient_aerosol_in_air
      • \n
    • Units: µg/m3
      • \n
    \n
    2. \n
\n
\n
⌨️ Input: Python
\n
print(dset.pm2p5_conc)\n
\n
\n
\n
🖥 Output
\n
   <xarray.DataArray 'pm2p5_conc' (time: 97, level: 1, latitude: 400, longitude: 700)>\n[27160000 values with dtype=float32]\n   Coordinates:\n* longitude  (longitude) float32 335.0 335.1 335.2 335.4 ... 44.75 44.85 44.95\n* latitude   (latitude) float32 69.95 69.85 69.75 69.65 ... 30.25 30.15 30.05\n* level      (level) float32 0.0\n* time       (time) timedelta64[ns] 00:00:00 01:00:00 ... 4 days 00:00:00\n   Attributes:\n       species:        PM2.5 Aerosol\n       units:          µg/m3\n       value:          hourly values\n       standard_name:  mass_concentration_of_pm2p5_ambient_aerosol_in_air\n
\n
\n
\n\n
\n
💬 Comment: Different ways to access Data variables
\n

To access a variable or coordinate, we can use “.” or specify its name as a string between squared brackets “[” “]”. For example:

\n
print(dset['pm2p5_conc'])\n# or alternatively\nprint(dset.pm2p5_conc)\n
\n

When we print a variable or coordinate, we do not get all the individual values but a DataArray that contains a lot of very useful metadata such as coordinates (if they have some), all the attributes such as the name, the physical units, etc.

\n
\n

Select / Subset from coordinates

\n

We often want to select elements from the coordinates for instance to subset a geographical area or select specific times or a specific time range.

\n

There are two different ways to make a selection:

\n\n

Select elements from coordinate by index

\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-11", "source": [ "print(dset.isel(time=0))" ], "cell_type": "code", "execution_count": null, "outputs": [ ], "metadata": { "attributes": { "classes": [ "python" ], "id": "" } } }, { "id": "cell-12", "source": "

You should see that the coordinate time “disappeared” from the Dimensions and now the variable pm2p5_conc is a 3D field with longitude, latitude and level.

\n

Select elements from coordinates by value

\n

When selecting elements by the value of the coordinate, we need to use the same datatype. For instance, to select an element from\ntime, we need to use timedelta64. The code below will give the same result as isel(time=0).

\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-13", "source": [ "print(dset.sel(time=np.timedelta64(0)))" ], "cell_type": "code", "execution_count": null, "outputs": [ ], "metadata": { "attributes": { "classes": [ "python" ], "id": "" } } }, { "id": "cell-14", "source": "

The output will be very similar to what we did previously when selecting from coordinates by index.

\n
\n
❓ Question: Select a single time for PM2.5
\n

How to select the forecast for December, 24th 2021 at 12:00 UTC?

\n
👁 View solution\n
👁 Solution
\n

Data starts on December, 22nd 2021 at 00:00 UTC so we need to add 2 days and 12 hours to select the correct time index.

\n
\n
⌨️ Input: Python
\n
print(dset.sel(time=(np.timedelta64(2,'D')+ np.timedelta64(12,'h'))))\n
\n
\n
\n
🖥 Output
\n
<xarray.Dataset>\nDimensions:     (longitude: 700, latitude: 400, level: 1)\nCoordinates:\n* longitude   (longitude) float32 -24.95 -24.85 -24.75 ... 44.75 44.85 44.95\n* latitude    (latitude) float32 69.95 69.85 69.75 69.65 ... 30.25 30.15 30.05\n* level       (level) float32 0.0\ntime        timedelta64[ns] 2 days 12:00:00\nData variables:\n    pm2p5_conc  (level, latitude, longitude) float32 0.4499 0.4421 ... 10.71\nAttributes:\n    title:        PM25 Air Pollutant FORECAST at the Surface\n    institution:  Data produced by Meteo France\nsource:       Data from ENSEMBLE model\nhistory:      Model ENSEMBLE FORECAST\n    FORECAST:     Europe, 20211222+[0H_96H]\n    summary:      ENSEMBLE model hourly FORECAST of PM25 concentration at the...\n    project:      MACC-RAQ (http://macc-raq.gmes-atmosphere.eu)\n
\n
\n
\n\n

Plotting

\n

To plot a map, you need to select a variable with data on geographical coordinates (latitude, longitude). In addition, coordinates need to be sorted, and preferably in increasing order. This is not the case for the coordinate “longitude” which is given between 360 and 0.

\n

Let’s shift the longitudes by 180 degrees so that they come in the range of -180 to 180.

\n

Shift longitudes

\n

We print the longitudes before and after shifting them so we can see what is happening.

\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-15", "source": [ "print(dset.longitude)" ], "cell_type": "code", "execution_count": null, "outputs": [ ], "metadata": { "attributes": { "classes": [ "python" ], "id": "" } } }, { "id": "cell-16", "source": "

The longitude values are between 335.05 and 44.95 degrees.

\n

Let’s now shift the longitudes to get values between -180, 180 degrees.

\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-17", "source": [ "dset.coords['longitude'] = (dset['longitude'] + 180) % 360 - 180\n", "print(dset.longitude)" ], "cell_type": "code", "execution_count": null, "outputs": [ ], "metadata": { "attributes": { "classes": [ "python" ], "id": "" } } }, { "id": "cell-18", "source": "

Indeed, the longitudes have been shifted and now the values are between -24.95 and 44.95.

\n

Visualize on a map PM2.5 for December, 24th 2021 at 12:00 UTC

\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-19", "source": [ "dset.sel(time=(np.timedelta64(2,'D')+ np.timedelta64(12,'h'))).pm2p5_conc.plot()" ], "cell_type": "code", "execution_count": null, "outputs": [ ], "metadata": { "attributes": { "classes": [ "python" ], "id": "" } } }, { "id": "cell-20", "source": "

We will get a figure like the one below:

\n

\"CAMS

\n
\n
level\" style=\"font-size: 150%\">💬 Comment: What about level
\n

Note that in the previous plot, we did not need to select level because there is one value only. However, if we had more than one level, we would need to add a selection on the level before plotting

\n
\n

Customize your plot

\n

There are many ways to customize your plots and we will only detail what we think is important for creating publication ready figures:

\n\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-21", "source": [ "fig = plt.figure(1, figsize=[15,10])\n", "\n", "# We're using cartopy to project our data.\n", "# (see documentation on cartopy)\n", "ax = plt.subplot(1, 1, 1, projection=ccrs.Mercator())\n", "ax.coastlines(resolution='10m')\n", "\n", "# We need to project our data to the new projection and for this we use `transform`.\n", "# we set the original data projection in transform (here PlateCarree)\n", "dset.sel(time=(np.timedelta64(2,'D') + np.timedelta64(12,'h')))['pm2p5_conc'].plot(ax=ax,\n", " transform=ccrs.PlateCarree(),\n", " vmin = 0, vmax = 35,\n", " cmap=cmc.roma_r)\n", "# One way to customize your title\n", "plt.title(\"Copernicus Atmosphere Monitoring Service PM2.5, 2 day forecasts\\n 24th December 2021 at 12:00 UTC\", fontsize=18)\n", "plt.savefig(\"CAMS-PM2_5-fc-20211224.png\")" ], "cell_type": "code", "execution_count": null, "outputs": [ ], "metadata": { "attributes": { "classes": [ "python" ], "id": "" } } }, { "id": "cell-22", "source": "
\n
\n

And you should get the following plot:

\n

\"Customized

\n

Multi-plots

\n

Now, we will plot several times on the same figure in different sub-plots; we will not plot all the times (too many) but the first 24 forecasted values.

\n

Firstly, we need to create a list of times and convert it to pandas datetime in order to make it easier to format times when plotting:

\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-23", "source": [ "list_times = np.datetime64('2021-12-22') + dset.time.sel(time=slice(np.timedelta64(0),np.timedelta64(1,'D')))\n", "print(pd.to_datetime(list_times).strftime(\"%d %b %H:%S UTC\"))" ], "cell_type": "code", "execution_count": null, "outputs": [ ], "metadata": { "attributes": { "classes": [ "python" ], "id": "" } } }, { "id": "cell-24", "source": "

Secondly, we need to use the same plotting method as earlier, but we pass additional parameters:

\n\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-25", "source": [ "fig = plt.figure(1, figsize=[10,10])\n", "\n", "# We're using cartopy to project our data.\n", "# (see documentation on cartopy)\n", "proj_plot = ccrs.Mercator()\n", "\n", "# We need to project our data to the new projection and for this we use `transform`.\n", "# we set the original data projection in transform (here PlateCarree)\n", "p = dset.sel(time=slice(np.timedelta64(1,'h'),np.timedelta64(1,'D')))['pm2p5_conc'].plot(transform=ccrs.PlateCarree(),\n", " vmin = 0, vmax = 35,\n", " subplot_kws={\"projection\": proj_plot},\n", " col='time', col_wrap=4,\n", " robust=True,\n", " aspect=dset.dims[\"longitude\"] / dset.dims[\"latitude\"], # for a sensible figsize\n", " cmap=cmc.roma_r)\n", "# We have to set the map's options on our axes\n", "for ax,i in zip(p.axes.flat, (np.datetime64('2021-12-22') + dset.time.sel(time=slice(np.timedelta64(0),np.timedelta64(1,'D')))).values):\n", " ax.coastlines('10m')\n", " ax.set_title(\"CAMS PM2.5 \" + pd.to_datetime(i).strftime(\"%d %b %H:%S UTC\"), fontsize=12)\n", "# Save your figure\n", "plt.savefig(\"CAMS-PM2_5-fc-multi.png\")" ], "cell_type": "code", "execution_count": null, "outputs": [ ], "metadata": { "attributes": { "classes": [ "python" ], "id": "" } } }, { "id": "cell-26", "source": "

In the second part of our plot, we are going to customize each subplot (this is why we loop for each of them and get their axes) by adding:

\n\n

\"Customized

\n
\n
❓ Question: PM2.5 over Italy
\n

Using a Multi-plot between Rome and Naples, can you tell us if the forecasted PM2.5 will increase or decrease during the first 24 hours?

\n
👁 View solution\n
👁 Solution
\n

We will select a sub-area: 11. East to 15.0 East and 40. N to 43. N. PM2.5 will increase and reach values close to 35 μm.m-3.\nWe will use slice to select the area and we slice latitudes with latitude=slice(47.3, 36.5) and not latitude=slice(36.5, 47.3).\nThe reason is that when using slice, you need to specify values using the same order as in the coordinates. Latitudes are specified in\ndecreasing order for CAMS.

\n
\n
⌨️ Input: Python
\n
fig = plt.figure(1, figsize=[10,10])\n\n# We're using cartopy to project our data.\n# (see documentation on cartopy)\nproj_plot = ccrs.Mercator()\n\n# We need to project our data to the new projection and for this we use <code>transform</code>.\n# we set the original data projection in transform (here PlateCarree)\np = dset.sel(time=slice(np.timedelta64(1,'h'),np.timedelta64(1,'D'))).sel(latitude=slice(43., 40.),\n                                                                          longitude=slice(11.,15.))['pm2p5_conc'].plot(transform=ccrs.PlateCarree(),\n                                                                                     vmin = 0, vmax = 35,\n                                                                                     subplot_kws={\"projection\": proj_plot},\n                                                                                     col='time', col_wrap=4,\n                                                                                     robust=True,\n                                                                                     aspect=dset.dims[\"longitude\"] / dset.dims[\"latitude\"],  # for a sensible figsize\n                                                                                     cmap=cmc.roma_r)\n# We have to set the map's options on all axes\nfor ax,i in zip(p.axes.flat,  (np.datetime64('2021-12-22') + dset.time.sel(time=slice(np.timedelta64(0),np.timedelta64(1,'D')))).values):\n    ax.coastlines('10m')\n    ax.set_title(\"CAMS PM2.5 \" + pd.to_datetime(i).strftime(\"%d %b %H:%S UTC\"), fontsize=12)\n# Save your figure\nplt.savefig(\"CAMS-PM2_5-fc-multi-Italy.png\")\n
\n
\n

\"PM2.5

\n
\n\n

How to use the where method

\n

Sometimes we may want to make more complex selections with criteria on the values of a given variable and not only on its coordinates. For this purpose, we use the where method. For instance, we may want to only keep PM2.5 if values are greater than 25 μm.m-3 (or any threshold you would like to choose).

\n

Mask values that do not meet a criteria with Where

\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-27", "source": [ "print(dset.where(dset['pm2p5_conc'] > 25))" ], "cell_type": "code", "execution_count": null, "outputs": [ ], "metadata": { "attributes": { "classes": [ "python" ], "id": "" } } }, { "id": "cell-28", "source": "
\n
💬 Comment: What happened?
\n

Each element of the dataset where the criteria within the where statement is not met, e.g. when PM2.5 <= 25, will be set to nan.\nYou may not see any changes when printing the dataset but if you look carefuly at pm2p5_conc values, you will see many nan.

\n
\n

Let’s plot one time to better see what happened:

\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-29", "source": [ "######################\n", "# Plotting with mask #\n", "######################\n", "\n", "fig = plt.figure(1, figsize=[15,10])\n", "\n", "# We're using cartopy to project our data.\n", "# (see documentation on cartopy)\n", "ax = plt.subplot(1, 1, 1, projection=ccrs.Mercator())\n", "ax.coastlines(resolution='10m')\n", "\n", "# We need to project our data to the new projection and for this we use `transform`.\n", "# we set the original data projection in transform (here PlateCarree)\n", "dset.where(dset['pm2p5_conc'] > 25).isel(time=0)['pm2p5_conc'].plot(ax=ax,\n", " transform=ccrs.PlateCarree(),\n", " vmin = 0, vmax = 35,\n", " cmap=cmc.roma_r)\n", "# One way to customize your title\n", "plt.title(\"Copernicus Atmosphere Monitoring Service PM2.5, 2 day forecasts\\n 24th December 2021 at 12:00 UTC\\n only values > 25\", fontsize=18)\n", "plt.savefig(\"CAMS-PM2_5-fc-20211224-25.png\")" ], "cell_type": "code", "execution_count": null, "outputs": [ ], "metadata": { "attributes": { "classes": [ "python" ], "id": "" } } }, { "id": "cell-30", "source": "

\"PM2.5

\n

We can then make the same multi-plot as earlier (over Italy) but with a where statement to mask values lower than 25 μm.m-3:

\n

Multi-plot over Italy using a mask

\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-31", "source": [ "fig = plt.figure(1, figsize=[10,10])\n", "\n", "# We're using cartopy to project our data.\n", "# (see documentation on cartopy)\n", "proj_plot = ccrs.Mercator()\n", "\n", "# We need to project our data to the new projection and for this we use `transform`.\n", "# we set the original data projection in transform (here PlateCarree)\n", "p = dset.where(dset['pm2p5_conc'] > 25).sel(time=slice(np.timedelta64(1,'h'),np.timedelta64(1,'D'))).sel(latitude=slice(43., 40.),\n", " longitude=slice(11.,15.))['pm2p5_conc'].plot(transform=ccrs.PlateCarree(),\n", " vmin = 0, vmax = 35,\n", " subplot_kws={\"projection\": proj_plot},\n", " col='time', col_wrap=4,\n", " robust=True,\n", " aspect=dset.dims[\"longitude\"] / dset.dims[\"latitude\"], # for a sensible figsize\n", " cmap=cmc.roma_r)\n", "# We have to set the map's options on all four axes\n", "for ax,i in zip(p.axes.flat, (np.datetime64('2021-12-22') + dset.time.sel(time=slice(np.timedelta64(0),np.timedelta64(1,'D')))).values):\n", " ax.coastlines('10m')\n", " ax.set_title(\"PM2.5 > 25 μm.m-3\" + pd.to_datetime(i).strftime(\"%d %b %H:%S UTC\"), fontsize=12)\n", "# Save your figure\n", "plt.savefig(\"CAMS-PM2_5-fc-multi-Italy-25.png\")" ], "cell_type": "code", "execution_count": null, "outputs": [ ], "metadata": { "attributes": { "classes": [ "python" ], "id": "" } } }, { "id": "cell-32", "source": "

\"Multi-plot

\n

Reduction operations

\n

We often want to compute the mean of all our datasets, or along a dimension (for instance time). If you do not pass any argument to the operation then it is done over all dimensions.

\n

Details on the mean method

\n

When we do not specify any parameters, we get a single value.

\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-33", "source": [ "print(dset.sel(latitude=slice(43., 40.), longitude=slice(11.,15.)).mean())" ], "cell_type": "code", "execution_count": null, "outputs": [ ], "metadata": { "attributes": { "classes": [ "python" ], "id": "" } } }, { "id": "cell-34", "source": "
\n
🖥 Output
\n
<xarray.Dataset>\nDimensions:     ()\n Data variables:\n    pm2p5_conc  float32 9.118\n
\n
\n
\n
❓ Question: Maximum PM2.5 over Italy
\n

What is the maximum forecasted PM2.5 value over the Rome-Naples region?

\n
👁 View solution\n
👁 Solution
\n

We select the same sub-area: 11. East to 15.0 East and 40. N to 43. N and compute the maximum with max. The maximum PM2.5 value is 59.13694382 μm.m-3 (that is rounded to 59.14).

\n
\n
⌨️ Input: Python
\n
dset.sel(latitude=slice(43., 40.), longitude=slice(11.,15.)).max()\n
\n
\n
\n
🖥 Output
\n
xarray.Dataset\nDimensions:\nCoordinates: (0)\nData variables:\npm2p5_conc\n()\nfloat64\n59.14\narray(59.13694382)\nAttributes: (0)\n
\n
\n
\n\n
\n
❓ Question: Find when the maximum PM2.5 is forecasted
\n

When is the maximum PM2.5 value forecasted?

\n
👁 View solution\n
👁 Solution
\n

We will select a sub-area: 11. East to 15.0 East and 40. N to 43. N and average over the entire selected area and search where the maximum PM2.5 value of 59.13694382 μm.m-3 is found. The maximum PM2.5 value occurs on 2021-12-22 at 20:00 UTC.

\n
\n
⌨️ Input: Python
\n
dset_tmean = dset.sel(latitude=slice(43., 40.), longitude=slice(11.,15.)).max(dim=('latitude', 'longitude'))\ndset_tmean_max = dset_tmean.where(dset_tmean['pm2p5_conc'] == 59.13694382, drop=True)\nprint(dset_tmean_max)\n
\n
\n
\n
🖥 Output
\n
<xarray.Dataset>\nDimensions:     (time: 1, level: 1)\nCoordinates:\n* level       (level) float32 0.0\n* time        (time) timedelta64[ns] 20:00:00\nData variables:\n    pm2p5_conc  (time, level) float32 59.14\n
\n
\n
\n\n
\n
💬 Comment: Pixel size when averaging
\n

We average over a relatively small area so we do not make a weighted average. Use weighted averages when averaging over the entire globe or over a large area where the pixel sizes may vary (depending on the latitude).

\n
\n

Details on the resample method

\n

1 day Resampling

\n

The resampling frequency is lower than our original data, so we would need to apply a global operation on the data we group together such as mean, min, max:

\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-35", "source": [ "print(dset.resample(time='1D').mean())" ], "cell_type": "code", "execution_count": null, "outputs": [ ], "metadata": { "attributes": { "classes": [ "python" ], "id": "" } } }, { "id": "cell-36", "source": "
\n
🖥 Output
\n
<xarray.Dataset>\nDimensions:     (time: 5, longitude: 700, latitude: 400, level: 1)\nCoordinates:\n* time        (time) timedelta64[ns] 0 days 1 days 2 days 3 days 4 days\n* longitude   (longitude) float32 -24.95 -24.85 -24.75 ... 44.75 44.85 44.95\n* latitude    (latitude) float32 69.95 69.85 69.75 69.65 ... 30.25 30.15 30.05\n* level       (level) float32 0.0\nData variables:\n    pm2p5_conc  (time, level, latitude, longitude) float32 0.4298 ... 7.501\n
\n
\n

30 minute resampling

\n

When the resampling frequency is higher than the original data, we need to indicate how to fill the gaps, for instance, interpolate and indicate which interpolation method to apply or select nearest values, etc.:

\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "id": "cell-37", "source": [ "print(dset.resample(time='30min').interpolate('linear'))" ], "cell_type": "code", "execution_count": null, "outputs": [ ], "metadata": { "attributes": { "classes": [ "python" ], "id": "" } } }, { "id": "cell-38", "source": "
\n
💬 Comment: Be careful when sub-sampling!
\n

Increasing the frequency of your data e.g. artificially creating data may not be scientifically relevant. Please use it carefully! Interpolating is not always scientifically relevant and sometimes you may prefer to choose a different method, like taking the nearest value for instance:

\n
\n
⌨️ Input: Python
\n
 dset.resample(time='30min').nearest()\n
\n
\n
\n
\n
❓ Question: PM2.5 over Italy in the next 4 days
\n

Using a Multi-plot between Rome and Naples, and making averages per day, can you tell us if forecasted PM2.5 will increase or decrease?

\n
👁 View solution\n
👁 Solution
\n

PM2.5 over Italy is overall decreasing over the next 4 forecasted days.

\n
\n
⌨️ Input: Python
\n
fig = plt.figure(1, figsize=[10,10])\n\n# We're using cartopy to project our data.\n# (see documentation on cartopy)\nproj_plot = ccrs.Mercator()\n\nsub_dset = dset.sel(latitude=slice(43., 40.), longitude=slice(11.,15.)).resample(time='1D').mean()\n# We need to project our data to the new projection and for this we use <code>transform</code>.\n# we set the original data projection in transform (here PlateCarree)\np = sub_dset['pm2p5_conc'].plot(transform=ccrs.PlateCarree(),\n                                vmin = 0, vmax = 35,\n                                subplot_kws={\"projection\": proj_plot},\n                                col='time', col_wrap=5,\n                                robust=True,\n                                aspect=dset.dims[\"longitude\"] / dset.dims[\"latitude\"],  # for a sensible figsize\n                                cmap=cmc.roma_r)\n# We have to set the map's options on all axes\nfor ax,i in zip(p.axes.flat,  (np.datetime64('2021-12-22') + dset.time.sel(time=slice(np.timedelta64(0),np.timedelta64(1,'D')))).values):\n    ax.coastlines('10m')\n    ax.set_title(\"CAMS PM2.5 \" + pd.to_datetime(i).strftime(\"%d %b %H:%S UTC\"), fontsize=12)\n# Save your figure\nplt.savefig(\"CAMS-PM2_5-fc-multi-Italy-mean-per-day.png\")\n
\n
\n

\"Daily

\n
\n\n
\n
Grouby versus resample\" style=\"font-size: 150%\">💬 Comment: Grouby versus resample
\n

Use groupby instead of resample when you wish to group over a dimension that is not time. groupby is very similar to resample but can be applied to any coordinates and not only to time.

\n
\n

Conclusion

\n

Well done! Pangeo is a fantastic community with many more resources for learning and/or contributing! Please, if you use any Python packages from the Pangeo ecosystem, do not forget to cite Pangeo Abernathey et al. 2017, Abernathey et al. 2021, Gentemann et al. 2021 and Sambasivan et al. 2021!

\n

Have a look at the Pangeo Tutorial Gallery to pick up your next Pangeo training material!

\n", "cell_type": "markdown", "metadata": { "editable": false, "collapsed": false } }, { "cell_type": "markdown", "id": "final-ending-cell", "metadata": { "editable": false, "collapsed": false }, "source": [ "# Key Points\n\n", "- Pangeo ecosystem enables big data analysis in geosciences\n", "- Xarray is an important Python package for big data analysis in geosciences\n", "- Xarray can be used to read, select, mask and plot netCDF data\n", "- Xarray can also be used to perform global operations such as mean, max, min or resampling data\n", "\n# Congratulations on successfully completing this tutorial!\n\n", "Please [fill out the feedback on the GTN website](https://training.galaxyproject.org/training-material/topics/climate/tutorials/pangeo-notebook/tutorial.html#feedback) and check there for further resources!\n" ] } ] }