Some tips with xarray and pandas

Contents

Some tips with xarray and pandas#

  • We have massively different levels here

  • Try to make some aims for technical skills you can learn!

    • If you are beginning with python → learn the basics

    • If you are good at basic python → learn new packages and tricks

    • If you know all the packages → improve your skills with producing your own software, organising your code etc.

    • If you don’t know git and github → get better at this!

  • Learn from each other!

  • Every time I taught on this course, I learned something!

  • AI is a real game changer and we are figuring out together how to use it well!

    • Especially great to give us feedback on this.

Please feel free to come with suggestions and extra input as we go!

CHANGE: Etherpad for the course

Note to Sara: do menti

Jupyter notebooks:#

Jupyter is an interactive computing environment that allows users to write, run, and document code in a single interface, typically using notebooks that combine live code, text, equations, and visualizations. It supports multiple programming languages (like Python, R, and Julia) and is widely used in data science, research, and education for reproducible analysis and exploration.

a=1
import numpy as np
a
1
  • cells can be markdown or code

  • run cell by shift+enter, cmd+enter, ctrl+enter

  • Kernel: You can restart, interrupt etc.

  • Tricks:

    • mark multiple places at the same time

    • cmd or ctr + D: select same snippet.

    • Escape to move out of cell

      • Once out of cell mode, you can copy (c), paste (v), insert cell above (a), incert cell below (b) etc.

  • Contextual help: cmd+i or ctrl+i

.

.

.

.

.

What are pandas and xarray?#

  • Pandas → like a spreadsheet 2D data with columns and rows

  • xarray → like pandas, but in N dimensions

  • Use the functionality these packages gives you! Will help you avoid mistakes. Try to get as good as possible :)

Pandas#

image.png (Source: https://www.geeksforgeeks.org/python-pandas-dataframe/)

Xarray#

image.png (Source: https://docs.xarray.dev/)

.

.

.

.

.

1. Import python packages#

import xarray as xr
xr.set_options(display_style='html')
import intake
import cftime
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
import pandas as pd
import datetime
import seaborn as sns
%matplotlib inline

2. Pandas basics:#

Can read a range of data formats, e.g.

pd.read_csv('filename.csv')
pd.read_excel('filename')
pd.read_sql('filename')

See overview: https://pandas.pydata.org/pandas-docs/stable/user_guide/io.html

#fn = '~/shared-craas2-ns9988k-ns9252k/escience2025/data/group3/eScience2024/zeppelin-ebc-2015-2019-1.csv'
fn = '~/shared-clinfra-ns12127k-ns9252k/eScience-courses/escience2024/data/data_group4/zeppelin-ebc-2015-2019-1.csv'

Read in as a DataFrame#

import pandas as pd
df = pd.read_csv(fn, index_col=0, parse_dates=True )
df
CVI_status BC_cvi z BC_tot flag visibility
time
2015-11-13 00:00:00 off 0.00 0.30583 0.04 bad 34596.2
2015-11-13 00:01:00 off 0.00 0.43400 0.02 bad 34481.7
2015-11-13 00:02:00 off 0.04 0.53317 0.01 bad 34572.3
2015-11-13 00:03:00 off 0.05 0.60250 0.02 bad 34061.8
2015-11-13 00:04:00 off -0.01 1.16433 0.02 bad 34461.8
... ... ... ... ... ... ...
2019-11-15 10:30:00 off 0.02 1.14300 0.02 good 36983.4
2019-11-15 10:31:00 off 0.00 1.34233 0.01 good 38964.4
2019-11-15 10:32:00 off 0.05 0.68667 -0.00 good 45981.7
2019-11-15 10:33:00 off -0.02 1.00450 -0.00 good 48805.2
2019-11-15 10:34:00 off 0.04 0.86300 0.01 good 47112.9

2107355 rows × 6 columns

Select a column:#

One column is called a Series

df['BC_cvi']
time
2015-11-13 00:00:00    0.00
2015-11-13 00:01:00    0.00
2015-11-13 00:02:00    0.04
2015-11-13 00:03:00    0.05
2015-11-13 00:04:00   -0.01
                       ... 
2019-11-15 10:30:00    0.02
2019-11-15 10:31:00    0.00
2019-11-15 10:32:00    0.05
2019-11-15 10:33:00   -0.02
2019-11-15 10:34:00    0.04
Name: BC_cvi, Length: 2107355, dtype: float64

Select data in general:#

Several ways to do this:

df.loc[<from_index>: <to_index>, <from_column>:<to_column>]
df.loc['2016-01-01':'2016-01-02', 'CVI_status':'BC_cvi']#['BC_cvi']
CVI_status BC_cvi
time
2016-01-01 00:00:00 off 0.13
2016-01-01 00:01:00 off 0.11
2016-01-01 00:02:00 off 0.03
2016-01-01 00:03:00 off 0.03
2016-01-01 00:04:00 off 0.00
... ... ...
2016-01-02 23:55:00 off 0.02
2016-01-02 23:56:00 off 0.00
2016-01-02 23:57:00 off 0.12
2016-01-02 23:58:00 off 0.00
2016-01-02 23:59:00 off 0.07

2880 rows × 2 columns

The same, but using indexes:

df.iloc[0:5, 0:2]#['BC_cvi']
CVI_status BC_cvi
time
2015-11-13 00:00:00 off 0.00
2015-11-13 00:01:00 off 0.00
2015-11-13 00:02:00 off 0.04
2015-11-13 00:03:00 off 0.05
2015-11-13 00:04:00 off -0.01

Plot something simple:#

df.loc['2016-01-01':'2016-02-01','BC_cvi'].plot()
<Axes: xlabel='time'>
../../../../_images/ec49e342857987b59759af32be8f4a410155a836f0a4596b41dd0ba77adfe71d.png

Save your dataframe:#

df.to_csv('my_new_dataset.csv')

3. Xarray and netcdf#

3.1 Reading in the data from a netcdf file:#

Netcdf format:#

  • NetCDF (Network Common Data Form) is a widely used, platform-independent file format designed for array-oriented scientific data.

  • It is self-describing, meaning metadata (e.g., units, descriptions) is stored alongside the data.

  • Supports multidimensional data structures, such as time-varying grids (e.g., time, latitude, longitude, level).

  • Commonly used in earth and atmospheric sciences for storing model outputs, reanalyses, and satellite data. Developed at UCAR.

  • Organizes content into:

    • Dimensions (e.g., time, lat, lon)

    • Variables (e.g., temperature, pressure)

    • Attributes (e.g., units, standard names, conventions)

  • Enables efficient storage and retrieval, even for very large datasets.

  • Supported by many scientific tools and libraries including xarray, netCDF4, NCO, and CDO.

  • Ideal for interoperability, reproducibility, and long-term data archiving.

NetCDF + xarray = <3

Xarray is a Python library designed to work with multi-dimensional arrays and datasets, particularly those used in earth sciences, climate science, and atmospheric science. It builds upon and extends the functionality of NumPy, Pandas, and NetCDF, providing a high-level interface for working with labeled, multi-dimensional data.

.

.

.

.

.

Open dataset#

  • Use xarray python package to analyze netCDF dataset

  • open_dataset allows to get all the metadata without loading data into memory.

  • with xarray, we only load into memory what is needed.

path='filename.nc'
ds = xr.open_dataset(path)
Opening multiple files:#

list_of_files = [
    'file1.nc',
    'file2.nc'
]
xr.open_mfdataset(list_of_files, concat_dim='time',combine='by_coords')

Example:#

(I’ve gone hunting in your data folders)

f1 = '~/shared-clinfra-ns12127k-ns9252k/eScience-courses/escience2025/data/group6/CMIP6/abrupt-4xCO2/rlut_Amon_CAMS-CSM1-0_abrupt-4xCO2_r1i1p1f1_gn.nc'
#f2 = '~/shared-craas2-ns9988k-ns9252k/escience2025/data/group6/CMIP6/abrupt-4xCO2/rsut_Amon_CAMS-CSM1-0_abrupt-4xCO2_r1i1p1f1_gn.nc'
f2 = '~/shared-clinfra-ns12127k-ns9252k/eScience-courses/escience2025/data/group6/CMIP6/abrupt-4xCO2/rsut_Amon_CAMS-CSM1-0_abrupt-4xCO2_r1i1p1f1_gn.nc'
Open one file:#
xr.open_dataset(f1)
/opt/conda/envs/pangeo/lib/python3.12/site-packages/dask/config.py:786: FutureWarning: Dask configuration key 'distributed.p2p.disk' has been deprecated; please use 'distributed.p2p.storage.disk' instead
  warnings.warn(
<xarray.Dataset> Size: 369MB
Dimensions:  (time: 1800, lat: 160, lon: 320)
Coordinates:
  * time     (time) object 14kB 0001-01-01 00:00:00 ... 0150-12-01 00:00:00
  * lat      (lat) float64 1kB -89.14 -88.03 -86.91 -85.79 ... 86.91 88.03 89.14
  * lon      (lon) float64 3kB 0.0 1.125 2.25 3.375 ... 355.5 356.6 357.8 358.9
Data variables:
    rlut     (time, lat, lon) float32 369MB ...
Attributes: (12/50)
    CDI:                    Climate Data Interface version 2.0.5 (https://mpi...
    Conventions:            CF-1.7 CMIP-6.2
    source:                 CAMS_CSM 1.0 (2016): \naerosol: none\natmos: ECHA...
    institution:            Chinese Academy of Meteorological Sciences, Beiji...
    activity_id:            CMIP
    branch_method:          Standard
    ...                     ...
    tracking_id:            hdl:21.14100/73320c27-d385-493c-9674-1470acacb8ed
    variable_id:            rlut
    variant_label:          r1i1p1f1
    license:                CMIP6 model data produced by Lawrence Livermore P...
    cmor_version:           3.4.0
    CDO:                    Climate Data Operators version 2.0.5 (https://mpi...
Open multiple files:#
ds1 = xr.open_mfdataset([f1, f2])

3.1 Check how your dataset looks#

Notice:

  • Coordinates

  • Data variables (each with assigned coordinates)

  • Attributes

ds1
<xarray.Dataset> Size: 737MB
Dimensions:  (time: 1800, lat: 160, lon: 320)
Coordinates:
  * time     (time) object 14kB 0001-01-01 00:00:00 ... 0150-12-01 00:00:00
  * lat      (lat) float64 1kB -89.14 -88.03 -86.91 -85.79 ... 86.91 88.03 89.14
  * lon      (lon) float64 3kB 0.0 1.125 2.25 3.375 ... 355.5 356.6 357.8 358.9
Data variables:
    rlut     (time, lat, lon) float32 369MB dask.array<chunksize=(1, 160, 320), meta=np.ndarray>
    rsut     (time, lat, lon) float32 369MB dask.array<chunksize=(1, 160, 320), meta=np.ndarray>
Attributes: (12/50)
    CDI:                    Climate Data Interface version 2.0.5 (https://mpi...
    Conventions:            CF-1.7 CMIP-6.2
    source:                 CAMS_CSM 1.0 (2016): \naerosol: none\natmos: ECHA...
    institution:            Chinese Academy of Meteorological Sciences, Beiji...
    activity_id:            CMIP
    branch_method:          Standard
    ...                     ...
    tracking_id:            hdl:21.14100/73320c27-d385-493c-9674-1470acacb8ed
    variable_id:            rlut
    variant_label:          r1i1p1f1
    license:                CMIP6 model data produced by Lawrence Livermore P...
    cmor_version:           3.4.0
    CDO:                    Climate Data Operators version 2.0.5 (https://mpi...

.

.

.

.

.

3.2 Read in CMIP6 data: We will skip this next part, but you can check it later to read data:#

cat_url = "https://storage.googleapis.com/cmip6/pangeo-cmip6.json"
#cat_url = '/mnt/clinfra-ns12127k/data/catalogs/noresm.json'#craas2-ns9988k/data/catalogs/cmip6.json'
col = intake.open_esm_datastore(cat_url)
col

pangeo-cmip6 catalog with 7674 dataset(s) from 514818 asset(s):

unique
activity_id 18
institution_id 36
source_id 88
experiment_id 170
member_id 657
table_id 37
variable_id 700
grid_label 10
zstore 514818
dcpp_init_year 61
version 736
derived_variable_id 0

Search corresponding data#

Please check here for info about CMIP and variables :)

Particularly useful is maybe the variable search which you find here: https://clipc-services.ceda.ac.uk/dreq/mipVars.html

cat = col.search(source_id = ['CESM2'], 
                 experiment_id=['historical'], 
                 table_id=['Amon','fx','AERmon'], 
                 variable_id=['tas','hurs', 'areacella','mmrso4' ], 
                 member_id=['r1i1p1f1'],
                )
cat.df
activity_id institution_id source_id experiment_id member_id table_id variable_id grid_label zstore dcpp_init_year version
0 CMIP NCAR CESM2 historical r1i1p1f1 fx areacella gn gs://cmip6/CMIP6/CMIP/NCAR/CESM2/historical/r1... <NA> 20190308
1 CMIP NCAR CESM2 historical r1i1p1f1 AERmon mmrso4 gn gs://cmip6/CMIP6/CMIP/NCAR/CESM2/historical/r1... <NA> 20190308
2 CMIP NCAR CESM2 historical r1i1p1f1 Amon hurs gn gs://cmip6/CMIP6/CMIP/NCAR/CESM2/historical/r1... <NA> 20190308
3 CMIP NCAR CESM2 historical r1i1p1f1 Amon tas gn gs://cmip6/CMIP6/CMIP/NCAR/CESM2/historical/r1... <NA> 20190308
cat.esmcat.aggregation_control.groupby_attrs = ['activity_id','experiment_id', 'source_id','table_id','grid_label']
cat.esmcat.aggregation_control.groupby_attrs
['activity_id', 'experiment_id', 'source_id', 'table_id', 'grid_label']

Create dictionary from the list of datasets we found#

  • This step may take several minutes so be patient!

dset_dict = cat.to_dataset_dict(zarr_kwargs={'use_cftime':True})
--> The keys in the returned dictionary of datasets are constructed as follows:
	'activity_id.experiment_id.source_id.table_id.grid_label'
100.00% [3/3 00:29<00:00]
/opt/conda/envs/pangeo/lib/python3.12/site-packages/intake_esm/source.py:308: FutureWarning: In a future version of xarray the default value for compat will change from compat='no_conflicts' to compat='override'. This is likely to lead to different results when combining overlapping variables with the same name. To opt in to new defaults and get rid of these warnings now use `set_options(use_new_combine_kwarg_defaults=True) or set compat explicitly.
  self._ds = xr.combine_by_coords(
/opt/conda/envs/pangeo/lib/python3.12/site-packages/intake_esm/source.py:308: FutureWarning: In a future version of xarray the default value for compat will change from compat='no_conflicts' to compat='override'. This is likely to lead to different results when combining overlapping variables with the same name. To opt in to new defaults and get rid of these warnings now use `set_options(use_new_combine_kwarg_defaults=True) or set compat explicitly.
  self._ds = xr.combine_by_coords(
/opt/conda/envs/pangeo/lib/python3.12/site-packages/intake_esm/source.py:308: FutureWarning: In a future version of xarray the default value for compat will change from compat='no_conflicts' to compat='override'. This is likely to lead to different results when combining overlapping variables with the same name. To opt in to new defaults and get rid of these warnings now use `set_options(use_new_combine_kwarg_defaults=True) or set compat explicitly.
  self._ds = xr.combine_by_coords(
list(dset_dict.keys())
['CMIP.historical.CESM2.fx.gn',
 'CMIP.historical.CESM2.AERmon.gn',
 'CMIP.historical.CESM2.Amon.gn']
Since I have already checked that these datasets are on the same grid, we can merge them:
ds_list =[]
for k in dset_dict.keys():
    ds = dset_dict[k]
    for v in ['lon_bnds', 'lat_bnds', 'time_bnds']:
        if v in ds:
            ds= ds.drop_vars(v)
    ds_list.append(ds)     
ds = xr.merge(ds_list,compat='override')

.

.

.

.

.

4. Check how your dataset looks#

Notice:

  • Coordinates

  • Data variables (each with assigned coordinates)

  • Attributes

ds
<xarray.Dataset> Size: 15GB
Dimensions:         (member_id: 1, dcpp_init_year: 1, lat: 192, lon: 288,
                     lev: 32, time: 1980, nbnd: 2)
Coordinates:
  * member_id       (member_id) object 8B 'r1i1p1f1'
  * dcpp_init_year  (dcpp_init_year) object 8B None
  * lat             (lat) float64 2kB -90.0 -89.06 -88.12 ... 88.12 89.06 90.0
  * lon             (lon) float64 2kB 0.0 1.25 2.5 3.75 ... 356.2 357.5 358.8
  * lev             (lev) float64 256B -3.643 -7.595 -14.36 ... -976.3 -992.6
  * time            (time) object 16kB 1850-01-15 12:00:00 ... 2014-12-15 12:...
    a_bnds          (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
    b_bnds          (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
    lev_bnds        (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
Dimensions without coordinates: nbnd
Data variables:
    areacella       (member_id, dcpp_init_year, lat, lon) float32 221kB dask.array<chunksize=(1, 1, 192, 288), meta=np.ndarray>
    a               (lev) float64 256B dask.array<chunksize=(32,), meta=np.ndarray>
    b               (lev) float64 256B dask.array<chunksize=(32,), meta=np.ndarray>
    mmrso4          (member_id, dcpp_init_year, time, lev, lat, lon) float32 14GB dask.array<chunksize=(1, 1, 10, 32, 192, 288), meta=np.ndarray>
    p0              float32 4B ...
    ps              (time, lat, lon) float32 438MB dask.array<chunksize=(600, 192, 288), meta=np.ndarray>
    hurs            (member_id, dcpp_init_year, time, lat, lon) float32 438MB dask.array<chunksize=(1, 1, 600, 192, 288), meta=np.ndarray>
    tas             (member_id, dcpp_init_year, time, lat, lon) float32 438MB dask.array<chunksize=(1, 1, 600, 192, 288), meta=np.ndarray>
Attributes: (12/58)
    Conventions:                      CF-1.7 CMIP-6.2
    activity_id:                      CMIP
    branch_method:                    standard
    branch_time_in_child:             674885.0
    branch_time_in_parent:            219000.0
    case_id:                          15
    ...                               ...
    intake_esm_attrs:variable_id:     areacella
    intake_esm_attrs:grid_label:      gn
    intake_esm_attrs:zstore:          gs://cmip6/CMIP6/CMIP/NCAR/CESM2/histor...
    intake_esm_attrs:version:         20190308
    intake_esm_attrs:_data_format_:   zarr
    intake_esm_dataset_key:           CMIP.historical.CESM2.fx.gn

.

.

.

.

.

5. Sometimes we want to do some nice tweaks before we start:#

5.1 Indexing/selecting data:#

  • Xarray loads data only when it needs to (it’s lazy, someone else can explain) – pandas loads everything right away

  • You might want to early on define the subset of data you want to look at so that you don’t end up loading a lot of extra data.

See here for nice overview#

In order to reduce unecessary calculations and loading of data, think about which part of the data you want, and slice early on.

Slice in time: the sel method#

dss = ds.sel(time = slice('1990-01-01','2010-01-01'))
dss
<xarray.Dataset> Size: 2GB
Dimensions:         (member_id: 1, dcpp_init_year: 1, lat: 192, lon: 288,
                     lev: 32, time: 240, nbnd: 2)
Coordinates:
  * member_id       (member_id) object 8B 'r1i1p1f1'
  * dcpp_init_year  (dcpp_init_year) object 8B None
  * lat             (lat) float64 2kB -90.0 -89.06 -88.12 ... 88.12 89.06 90.0
  * lon             (lon) float64 2kB 0.0 1.25 2.5 3.75 ... 356.2 357.5 358.8
  * lev             (lev) float64 256B -3.643 -7.595 -14.36 ... -976.3 -992.6
  * time            (time) object 2kB 1990-01-15 12:00:00 ... 2009-12-15 12:0...
    a_bnds          (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
    b_bnds          (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
    lev_bnds        (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
Dimensions without coordinates: nbnd
Data variables:
    areacella       (member_id, dcpp_init_year, lat, lon) float32 221kB dask.array<chunksize=(1, 1, 192, 288), meta=np.ndarray>
    a               (lev) float64 256B dask.array<chunksize=(32,), meta=np.ndarray>
    b               (lev) float64 256B dask.array<chunksize=(32,), meta=np.ndarray>
    mmrso4          (member_id, dcpp_init_year, time, lev, lat, lon) float32 2GB dask.array<chunksize=(1, 1, 10, 32, 192, 288), meta=np.ndarray>
    p0              float32 4B ...
    ps              (time, lat, lon) float32 53MB dask.array<chunksize=(120, 192, 288), meta=np.ndarray>
    hurs            (member_id, dcpp_init_year, time, lat, lon) float32 53MB dask.array<chunksize=(1, 1, 120, 192, 288), meta=np.ndarray>
    tas             (member_id, dcpp_init_year, time, lat, lon) float32 53MB dask.array<chunksize=(1, 1, 120, 192, 288), meta=np.ndarray>
Attributes: (12/58)
    Conventions:                      CF-1.7 CMIP-6.2
    activity_id:                      CMIP
    branch_method:                    standard
    branch_time_in_child:             674885.0
    branch_time_in_parent:            219000.0
    case_id:                          15
    ...                               ...
    intake_esm_attrs:variable_id:     areacella
    intake_esm_attrs:grid_label:      gn
    intake_esm_attrs:zstore:          gs://cmip6/CMIP6/CMIP/NCAR/CESM2/histor...
    intake_esm_attrs:version:         20190308
    intake_esm_attrs:_data_format_:   zarr
    intake_esm_dataset_key:           CMIP.historical.CESM2.fx.gn
NB: Always check that the function you are using is doing what you expect!
  • A lot of issues arise from people not doing sanity checks.

Example 2: You might want to select only the arctic e.g.:

dss_arctic = dss.sel(lat = slice(60,None))
dss_arctic
<xarray.Dataset> Size: 310MB
Dimensions:         (member_id: 1, dcpp_init_year: 1, lat: 32, lon: 288,
                     lev: 32, time: 240, nbnd: 2)
Coordinates:
  * member_id       (member_id) object 8B 'r1i1p1f1'
  * dcpp_init_year  (dcpp_init_year) object 8B None
  * lat             (lat) float64 256B 60.79 61.73 62.67 ... 88.12 89.06 90.0
  * lon             (lon) float64 2kB 0.0 1.25 2.5 3.75 ... 356.2 357.5 358.8
  * lev             (lev) float64 256B -3.643 -7.595 -14.36 ... -976.3 -992.6
  * time            (time) object 2kB 1990-01-15 12:00:00 ... 2009-12-15 12:0...
    a_bnds          (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
    b_bnds          (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
    lev_bnds        (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
Dimensions without coordinates: nbnd
Data variables:
    areacella       (member_id, dcpp_init_year, lat, lon) float32 37kB dask.array<chunksize=(1, 1, 32, 288), meta=np.ndarray>
    a               (lev) float64 256B dask.array<chunksize=(32,), meta=np.ndarray>
    b               (lev) float64 256B dask.array<chunksize=(32,), meta=np.ndarray>
    mmrso4          (member_id, dcpp_init_year, time, lev, lat, lon) float32 283MB dask.array<chunksize=(1, 1, 10, 32, 32, 288), meta=np.ndarray>
    p0              float32 4B ...
    ps              (time, lat, lon) float32 9MB dask.array<chunksize=(120, 32, 288), meta=np.ndarray>
    hurs            (member_id, dcpp_init_year, time, lat, lon) float32 9MB dask.array<chunksize=(1, 1, 120, 32, 288), meta=np.ndarray>
    tas             (member_id, dcpp_init_year, time, lat, lon) float32 9MB dask.array<chunksize=(1, 1, 120, 32, 288), meta=np.ndarray>
Attributes: (12/58)
    Conventions:                      CF-1.7 CMIP-6.2
    activity_id:                      CMIP
    branch_method:                    standard
    branch_time_in_child:             674885.0
    branch_time_in_parent:            219000.0
    case_id:                          15
    ...                               ...
    intake_esm_attrs:variable_id:     areacella
    intake_esm_attrs:grid_label:      gn
    intake_esm_attrs:zstore:          gs://cmip6/CMIP6/CMIP/NCAR/CESM2/histor...
    intake_esm_attrs:version:         20190308
    intake_esm_attrs:_data_format_:   zarr
    intake_esm_dataset_key:           CMIP.historical.CESM2.fx.gn

isel, sel: index selecting#

  • Select the surface (which in this case is the last index :)

dss
<xarray.Dataset> Size: 2GB
Dimensions:         (member_id: 1, dcpp_init_year: 1, lat: 192, lon: 288,
                     lev: 32, time: 240, nbnd: 2)
Coordinates:
  * member_id       (member_id) object 8B 'r1i1p1f1'
  * dcpp_init_year  (dcpp_init_year) object 8B None
  * lat             (lat) float64 2kB -90.0 -89.06 -88.12 ... 88.12 89.06 90.0
  * lon             (lon) float64 2kB 0.0 1.25 2.5 3.75 ... 356.2 357.5 358.8
  * lev             (lev) float64 256B -3.643 -7.595 -14.36 ... -976.3 -992.6
  * time            (time) object 2kB 1990-01-15 12:00:00 ... 2009-12-15 12:0...
    a_bnds          (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
    b_bnds          (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
    lev_bnds        (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
Dimensions without coordinates: nbnd
Data variables:
    areacella       (member_id, dcpp_init_year, lat, lon) float32 221kB dask.array<chunksize=(1, 1, 192, 288), meta=np.ndarray>
    a               (lev) float64 256B dask.array<chunksize=(32,), meta=np.ndarray>
    b               (lev) float64 256B dask.array<chunksize=(32,), meta=np.ndarray>
    mmrso4          (member_id, dcpp_init_year, time, lev, lat, lon) float32 2GB dask.array<chunksize=(1, 1, 10, 32, 192, 288), meta=np.ndarray>
    p0              float32 4B ...
    ps              (time, lat, lon) float32 53MB dask.array<chunksize=(120, 192, 288), meta=np.ndarray>
    hurs            (member_id, dcpp_init_year, time, lat, lon) float32 53MB dask.array<chunksize=(1, 1, 120, 192, 288), meta=np.ndarray>
    tas             (member_id, dcpp_init_year, time, lat, lon) float32 53MB dask.array<chunksize=(1, 1, 120, 192, 288), meta=np.ndarray>
Attributes: (12/58)
    Conventions:                      CF-1.7 CMIP-6.2
    activity_id:                      CMIP
    branch_method:                    standard
    branch_time_in_child:             674885.0
    branch_time_in_parent:            219000.0
    case_id:                          15
    ...                               ...
    intake_esm_attrs:variable_id:     areacella
    intake_esm_attrs:grid_label:      gn
    intake_esm_attrs:zstore:          gs://cmip6/CMIP6/CMIP/NCAR/CESM2/histor...
    intake_esm_attrs:version:         20190308
    intake_esm_attrs:_data_format_:   zarr
    intake_esm_dataset_key:           CMIP.historical.CESM2.fx.gn
dss_s = dss.isel(lev=-1)
dss_s
<xarray.Dataset> Size: 213MB
Dimensions:         (member_id: 1, dcpp_init_year: 1, lat: 192, lon: 288,
                     time: 240, nbnd: 2)
Coordinates:
  * member_id       (member_id) object 8B 'r1i1p1f1'
  * dcpp_init_year  (dcpp_init_year) object 8B None
  * lat             (lat) float64 2kB -90.0 -89.06 -88.12 ... 88.12 89.06 90.0
  * lon             (lon) float64 2kB 0.0 1.25 2.5 3.75 ... 356.2 357.5 358.8
  * time            (time) object 2kB 1990-01-15 12:00:00 ... 2009-12-15 12:0...
    a_bnds          (nbnd) float64 16B dask.array<chunksize=(2,), meta=np.ndarray>
    b_bnds          (nbnd) float64 16B dask.array<chunksize=(2,), meta=np.ndarray>
    lev             float64 8B -992.6
    lev_bnds        (nbnd) float64 16B dask.array<chunksize=(2,), meta=np.ndarray>
Dimensions without coordinates: nbnd
Data variables:
    areacella       (member_id, dcpp_init_year, lat, lon) float32 221kB dask.array<chunksize=(1, 1, 192, 288), meta=np.ndarray>
    a               float64 8B dask.array<chunksize=(), meta=np.ndarray>
    b               float64 8B dask.array<chunksize=(), meta=np.ndarray>
    mmrso4          (member_id, dcpp_init_year, time, lat, lon) float32 53MB dask.array<chunksize=(1, 1, 10, 192, 288), meta=np.ndarray>
    p0              float32 4B ...
    ps              (time, lat, lon) float32 53MB dask.array<chunksize=(120, 192, 288), meta=np.ndarray>
    hurs            (member_id, dcpp_init_year, time, lat, lon) float32 53MB dask.array<chunksize=(1, 1, 120, 192, 288), meta=np.ndarray>
    tas             (member_id, dcpp_init_year, time, lat, lon) float32 53MB dask.array<chunksize=(1, 1, 120, 192, 288), meta=np.ndarray>
Attributes: (12/58)
    Conventions:                      CF-1.7 CMIP-6.2
    activity_id:                      CMIP
    branch_method:                    standard
    branch_time_in_child:             674885.0
    branch_time_in_parent:            219000.0
    case_id:                          15
    ...                               ...
    intake_esm_attrs:variable_id:     areacella
    intake_esm_attrs:grid_label:      gn
    intake_esm_attrs:zstore:          gs://cmip6/CMIP6/CMIP/NCAR/CESM2/histor...
    intake_esm_attrs:version:         20190308
    intake_esm_attrs:_data_format_:   zarr
    intake_esm_dataset_key:           CMIP.historical.CESM2.fx.gn

.

.

.

.

.

5.2 Calculates variables and assign attributes!#

Nice for plotting and to keep track of what is in your dataset (especially ‘units’ and ‘standard_name’/’long_name’ will be looked for by xarray.

dss['T_C'] = dss['tas']-273.15
dss['T_C']
<xarray.DataArray 'T_C' (member_id: 1, dcpp_init_year: 1, time: 240, lat: 192,
                         lon: 288)> Size: 53MB
dask.array<sub, shape=(1, 1, 240, 192, 288), dtype=float32, chunksize=(1, 1, 120, 192, 288), chunktype=numpy.ndarray>
Coordinates:
  * member_id       (member_id) object 8B 'r1i1p1f1'
  * dcpp_init_year  (dcpp_init_year) object 8B None
  * time            (time) object 2kB 1990-01-15 12:00:00 ... 2009-12-15 12:0...
  * lat             (lat) float64 2kB -90.0 -89.06 -88.12 ... 88.12 89.06 90.0
  * lon             (lon) float64 2kB 0.0 1.25 2.5 3.75 ... 356.2 357.5 358.8
Attributes: (12/19)
    cell_measures:  area: areacella
    cell_methods:   area: time: mean
    comment:        near-surface (usually, 2 meter) air temperature
    description:    near-surface (usually, 2 meter) air temperature
    frequency:      mon
    id:             tas
    ...             ...
    time_label:     time-mean
    time_title:     Temporal mean
    title:          Near-Surface Air Temperature
    type:           real
    units:          K
    variable_id:    tas
#dss['T_C'] = dss['T_C'].assign_attrs({'units': 'deg C'})
dss['T_C'] = dss['T_C'].assign_attrs({'units': '$^\circ$ C'})
dss['T_C'] = dss['T_C'].assign_attrs({'long_name': 'Temperature at surface'})
<>:2: SyntaxWarning: invalid escape sequence '\c'
<>:2: SyntaxWarning: invalid escape sequence '\c'
/tmp/ipykernel_1197/3247141791.py:2: SyntaxWarning: invalid escape sequence '\c'
  dss['T_C'] = dss['T_C'].assign_attrs({'units': '$^\circ$ C'})
dss['T_C']
<xarray.DataArray 'T_C' (member_id: 1, dcpp_init_year: 1, time: 240, lat: 192,
                         lon: 288)> Size: 53MB
dask.array<sub, shape=(1, 1, 240, 192, 288), dtype=float32, chunksize=(1, 1, 120, 192, 288), chunktype=numpy.ndarray>
Coordinates:
  * member_id       (member_id) object 8B 'r1i1p1f1'
  * dcpp_init_year  (dcpp_init_year) object 8B None
  * time            (time) object 2kB 1990-01-15 12:00:00 ... 2009-12-15 12:0...
  * lat             (lat) float64 2kB -90.0 -89.06 -88.12 ... 88.12 89.06 90.0
  * lon             (lon) float64 2kB 0.0 1.25 2.5 3.75 ... 356.2 357.5 358.8
Attributes: (12/19)
    cell_measures:  area: areacella
    cell_methods:   area: time: mean
    comment:        near-surface (usually, 2 meter) air temperature
    description:    near-surface (usually, 2 meter) air temperature
    frequency:      mon
    id:             tas
    ...             ...
    time_label:     time-mean
    time_title:     Temporal mean
    title:          Near-Surface Air Temperature
    type:           real
    units:          $^\circ$ C
    variable_id:    tas
May always be small things you need to adjust:
dss['time']
<xarray.DataArray 'time' (time: 240)> Size: 2kB
array([cftime.DatetimeNoLeap(1990, 1, 15, 12, 0, 0, 0, has_year_zero=True),
       cftime.DatetimeNoLeap(1990, 2, 14, 0, 0, 0, 0, has_year_zero=True),
       cftime.DatetimeNoLeap(1990, 3, 15, 12, 0, 0, 0, has_year_zero=True),
       ...,
       cftime.DatetimeNoLeap(2009, 10, 15, 12, 0, 0, 0, has_year_zero=True),
       cftime.DatetimeNoLeap(2009, 11, 15, 0, 0, 0, 0, has_year_zero=True),
       cftime.DatetimeNoLeap(2009, 12, 15, 12, 0, 0, 0, has_year_zero=True)],
      shape=(240,), dtype=object)
Coordinates:
  * time     (time) object 2kB 1990-01-15 12:00:00 ... 2009-12-15 12:00:00
Attributes:
    axis:           T
    bounds:         time_bnds
    standard_name:  time
    title:          time
    type:           double
  • This calendar is in cftime and noLeap. Sometimes this causes issues when plotting timeseries, so just for fun we will convert to a standard datetime64 index calendar because it’s anyway monthly.

dss['time'] = dss['time'].to_dataframe().index.to_datetimeindex()
/tmp/ipykernel_1197/4195133999.py:1: FutureWarning: In a future version of xarray to_datetimeindex will default to returning a 'us'-resolution DatetimeIndex instead of a 'ns'-resolution DatetimeIndex. This warning can be silenced by explicitly passing the `time_unit` keyword argument.
  dss['time'] = dss['time'].to_dataframe().index.to_datetimeindex()
/tmp/ipykernel_1197/4195133999.py:1: RuntimeWarning: Converting a CFTimeIndex with dates from a non-standard calendar, 'noleap', to a pandas.DatetimeIndex, which uses dates from the standard calendar.  This may lead to subtle errors in operations that depend on the length of time between dates.
  dss['time'] = dss['time'].to_dataframe().index.to_datetimeindex()
dss['time']
<xarray.DataArray 'time' (time: 240)> Size: 2kB
array(['1990-01-15T12:00:00.000000000', '1990-02-14T00:00:00.000000000',
       '1990-03-15T12:00:00.000000000', ..., '2009-10-15T12:00:00.000000000',
       '2009-11-15T00:00:00.000000000', '2009-12-15T12:00:00.000000000'],
      shape=(240,), dtype='datetime64[ns]')
Coordinates:
  * time     (time) datetime64[ns] 2kB 1990-01-15T12:00:00 ... 2009-12-15T12:...
We get a warning, but it's ok as long as we know what we are doing.

3.3 Convert longitude:#

this data comes in 0–360 degrees, but often -180 to 180 is more convenient. So we can convert:

NOTE:

Maybe you want to put this in a module? Or a package..

dss
<xarray.Dataset> Size: 2GB
Dimensions:         (member_id: 1, dcpp_init_year: 1, lat: 192, lon: 288,
                     lev: 32, time: 240, nbnd: 2)
Coordinates:
  * member_id       (member_id) object 8B 'r1i1p1f1'
  * dcpp_init_year  (dcpp_init_year) object 8B None
  * lat             (lat) float64 2kB -90.0 -89.06 -88.12 ... 88.12 89.06 90.0
  * lon             (lon) float64 2kB 0.0 1.25 2.5 3.75 ... 356.2 357.5 358.8
  * lev             (lev) float64 256B -3.643 -7.595 -14.36 ... -976.3 -992.6
  * time            (time) datetime64[ns] 2kB 1990-01-15T12:00:00 ... 2009-12...
    a_bnds          (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
    b_bnds          (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
    lev_bnds        (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
Dimensions without coordinates: nbnd
Data variables:
    areacella       (member_id, dcpp_init_year, lat, lon) float32 221kB dask.array<chunksize=(1, 1, 192, 288), meta=np.ndarray>
    a               (lev) float64 256B dask.array<chunksize=(32,), meta=np.ndarray>
    b               (lev) float64 256B dask.array<chunksize=(32,), meta=np.ndarray>
    mmrso4          (member_id, dcpp_init_year, time, lev, lat, lon) float32 2GB dask.array<chunksize=(1, 1, 10, 32, 192, 288), meta=np.ndarray>
    p0              float32 4B ...
    ps              (time, lat, lon) float32 53MB dask.array<chunksize=(120, 192, 288), meta=np.ndarray>
    hurs            (member_id, dcpp_init_year, time, lat, lon) float32 53MB dask.array<chunksize=(1, 1, 120, 192, 288), meta=np.ndarray>
    tas             (member_id, dcpp_init_year, time, lat, lon) float32 53MB dask.array<chunksize=(1, 1, 120, 192, 288), meta=np.ndarray>
    T_C             (member_id, dcpp_init_year, time, lat, lon) float32 53MB dask.array<chunksize=(1, 1, 120, 192, 288), meta=np.ndarray>
Attributes: (12/58)
    Conventions:                      CF-1.7 CMIP-6.2
    activity_id:                      CMIP
    branch_method:                    standard
    branch_time_in_child:             674885.0
    branch_time_in_parent:            219000.0
    case_id:                          15
    ...                               ...
    intake_esm_attrs:variable_id:     areacella
    intake_esm_attrs:grid_label:      gn
    intake_esm_attrs:zstore:          gs://cmip6/CMIP6/CMIP/NCAR/CESM2/histor...
    intake_esm_attrs:version:         20190308
    intake_esm_attrs:_data_format_:   zarr
    intake_esm_dataset_key:           CMIP.historical.CESM2.fx.gn
from my_package import convert_lon_to_180, convert_lon_to_180
def convert_lon_to_180(ds):
    ds['lon'] = (ds['lon']+ 180) % 360 - 180
    ds = ds.sortby('lon')
    return ds

def convert_lon_to_360(ds):
    ds['lon'] =ds['lon']% 360
    ds = ds.sortby('lon')
    return lon % 360

(migth want to move this to a module!)

dss = convert_lon_to_180(dss)
dss
<xarray.Dataset> Size: 2GB
Dimensions:         (member_id: 1, dcpp_init_year: 1, lat: 192, lon: 288,
                     lev: 32, time: 240, nbnd: 2)
Coordinates:
  * member_id       (member_id) object 8B 'r1i1p1f1'
  * dcpp_init_year  (dcpp_init_year) object 8B None
  * lat             (lat) float64 2kB -90.0 -89.06 -88.12 ... 88.12 89.06 90.0
  * lon             (lon) float64 2kB -180.0 -178.8 -177.5 ... 176.2 177.5 178.8
  * lev             (lev) float64 256B -3.643 -7.595 -14.36 ... -976.3 -992.6
  * time            (time) datetime64[ns] 2kB 1990-01-15T12:00:00 ... 2009-12...
    a_bnds          (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
    b_bnds          (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
    lev_bnds        (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
Dimensions without coordinates: nbnd
Data variables:
    areacella       (member_id, dcpp_init_year, lat, lon) float32 221kB dask.array<chunksize=(1, 1, 192, 288), meta=np.ndarray>
    a               (lev) float64 256B dask.array<chunksize=(32,), meta=np.ndarray>
    b               (lev) float64 256B dask.array<chunksize=(32,), meta=np.ndarray>
    mmrso4          (member_id, dcpp_init_year, time, lev, lat, lon) float32 2GB dask.array<chunksize=(1, 1, 10, 32, 192, 288), meta=np.ndarray>
    p0              float32 4B ...
    ps              (time, lat, lon) float32 53MB dask.array<chunksize=(120, 192, 288), meta=np.ndarray>
    hurs            (member_id, dcpp_init_year, time, lat, lon) float32 53MB dask.array<chunksize=(1, 1, 120, 192, 288), meta=np.ndarray>
    tas             (member_id, dcpp_init_year, time, lat, lon) float32 53MB dask.array<chunksize=(1, 1, 120, 192, 288), meta=np.ndarray>
    T_C             (member_id, dcpp_init_year, time, lat, lon) float32 53MB dask.array<chunksize=(1, 1, 120, 192, 288), meta=np.ndarray>
Attributes: (12/58)
    Conventions:                      CF-1.7 CMIP-6.2
    activity_id:                      CMIP
    branch_method:                    standard
    branch_time_in_child:             674885.0
    branch_time_in_parent:            219000.0
    case_id:                          15
    ...                               ...
    intake_esm_attrs:variable_id:     areacella
    intake_esm_attrs:grid_label:      gn
    intake_esm_attrs:zstore:          gs://cmip6/CMIP6/CMIP/NCAR/CESM2/histor...
    intake_esm_attrs:version:         20190308
    intake_esm_attrs:_data_format_:   zarr
    intake_esm_dataset_key:           CMIP.historical.CESM2.fx.gn
  • Notice that it lost the units now..

dss['lon'].attrs['units'] = '$^\circ$ East'
<>:1: SyntaxWarning: invalid escape sequence '\c'
<>:1: SyntaxWarning: invalid escape sequence '\c'
/tmp/ipykernel_1197/4061564318.py:1: SyntaxWarning: invalid escape sequence '\c'
  dss['lon'].attrs['units'] = '$^\circ$ East'

Notice how the plotting labels use both the attribute “standard_name” and “units” from the dataset.

.

.

.

.

.

6. The easiest interpolation: select with ‘nearest’ neighboor#

Example: let’s select Zeppelin station: 78.906661, 11.889203

lat_zep = 78.90
lon_zep = 11.89
dss['T_C'].sel(lon=lon_zep, lat=lat_zep, method='nearest').plot()
[<matplotlib.lines.Line2D at 0x7f7eaa505610>]
../../../../_images/05a0fa67e6f72acf27283f3da564a255d9fc45947d1b2ce17d753f3c649227ab.png

7. Super quick averaging etc with xarray#

da_so4 = dss['mmrso4']
da_so4
<xarray.DataArray 'mmrso4' (member_id: 1, dcpp_init_year: 1, time: 240,
                            lev: 32, lat: 192, lon: 288)> Size: 2GB
dask.array<getitem, shape=(1, 1, 240, 32, 192, 288), dtype=float32, chunksize=(1, 1, 10, 32, 192, 288), chunktype=numpy.ndarray>
Coordinates:
  * member_id       (member_id) object 8B 'r1i1p1f1'
  * dcpp_init_year  (dcpp_init_year) object 8B None
  * time            (time) datetime64[ns] 2kB 1990-01-15T12:00:00 ... 2009-12...
  * lev             (lev) float64 256B -3.643 -7.595 -14.36 ... -976.3 -992.6
  * lat             (lat) float64 2kB -90.0 -89.06 -88.12 ... 88.12 89.06 90.0
  * lon             (lon) float64 2kB -180.0 -178.8 -177.5 ... 176.2 177.5 178.8
Attributes: (12/19)
    cell_measures:  area: areacella
    cell_methods:   area: time: mean
    comment:        Dry mass of sulfate (SO4) in aerosol particles as a fract...
    description:    Dry mass of sulfate (SO4) in aerosol particles as a fract...
    frequency:      mon
    id:             mmrso4
    ...             ...
    time_label:     time-mean
    time_title:     Temporal mean
    title:          Aerosol Sulfate Mass Mixing Ratio
    type:           real
    units:          kg kg-1
    variable_id:    mmrso4

Mean:

da_so4_lmean = da_so4.mean(['time','lon'], keep_attrs=True) #.plot()#ylim=[1000,100], yscale='log')
da_so4_lmean
<xarray.DataArray 'mmrso4' (member_id: 1, dcpp_init_year: 1, lev: 32, lat: 192)> Size: 25kB
dask.array<mean_agg-aggregate, shape=(1, 1, 32, 192), dtype=float32, chunksize=(1, 1, 32, 192), chunktype=numpy.ndarray>
Coordinates:
  * member_id       (member_id) object 8B 'r1i1p1f1'
  * dcpp_init_year  (dcpp_init_year) object 8B None
  * lev             (lev) float64 256B -3.643 -7.595 -14.36 ... -976.3 -992.6
  * lat             (lat) float64 2kB -90.0 -89.06 -88.12 ... 88.12 89.06 90.0
Attributes: (12/19)
    cell_measures:  area: areacella
    cell_methods:   area: time: mean
    comment:        Dry mass of sulfate (SO4) in aerosol particles as a fract...
    description:    Dry mass of sulfate (SO4) in aerosol particles as a fract...
    frequency:      mon
    id:             mmrso4
    ...             ...
    time_label:     time-mean
    time_title:     Temporal mean
    title:          Aerosol Sulfate Mass Mixing Ratio
    type:           real
    units:          kg kg-1
    variable_id:    mmrso4

8. Plotting with xarray#

  • Notice it’s slow now, why?

da_so4_lmean.plot()
<matplotlib.collections.QuadMesh at 0x7f7df23e0a40>
../../../../_images/67b9f2a334d2fc7785379f21bed92df3253f7b1510f8de9964ca175b619c4c3b.png
da_so4_lmean['lev'] = np.abs(da_so4_lmean['lev'].values)
da_so4_lmean.plot(
    ylim=[1000,100], 
    yscale='log'
)
<matplotlib.collections.QuadMesh at 0x7f7df0e74e30>
../../../../_images/8598ced3386adf66ab41fa84fb79f3e7c40b48c65672ed77d634177ddc076cfb.png

Standard deviation

dss['T_C'].std(['time']).plot()
<matplotlib.collections.QuadMesh at 0x7f7df046e8a0>
../../../../_images/cff9d963ce75d0979ef9a822c7ae29320c84eb07cda1917a08a94c35e78b01ce.png

Temperature change much stronger over land than ocean and higher at high latitudes…

9. Mask data and groupby: pick out seasons#

month = ds['time.month']
month
<xarray.DataArray 'month' (time: 1980)> Size: 16kB
array([ 1,  2,  3, ..., 10, 11, 12], shape=(1980,))
Coordinates:
  * time     (time) object 16kB 1850-01-15 12:00:00 ... 2014-12-15 12:00:00
dss_JA = (dss
    .where(month.isin([7,8]))
    .mean('time')
         )
dss_JA
<xarray.Dataset> Size: 8MB
Dimensions:         (member_id: 1, dcpp_init_year: 1, lat: 192, lon: 288,
                     lev: 32, nbnd: 2)
Coordinates:
  * member_id       (member_id) object 8B 'r1i1p1f1'
  * dcpp_init_year  (dcpp_init_year) object 8B None
  * lat             (lat) float64 2kB -90.0 -89.06 -88.12 ... 88.12 89.06 90.0
  * lon             (lon) float64 2kB -180.0 -178.8 -177.5 ... 176.2 177.5 178.8
  * lev             (lev) float64 256B -3.643 -7.595 -14.36 ... -976.3 -992.6
    a_bnds          (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
    b_bnds          (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
    lev_bnds        (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
Dimensions without coordinates: nbnd
Data variables:
    areacella       (member_id, dcpp_init_year, lat, lon) float32 221kB dask.array<chunksize=(1, 1, 192, 288), meta=np.ndarray>
    a               (lev) float64 256B dask.array<chunksize=(32,), meta=np.ndarray>
    b               (lev) float64 256B dask.array<chunksize=(32,), meta=np.ndarray>
    mmrso4          (member_id, dcpp_init_year, lev, lat, lon) float32 7MB dask.array<chunksize=(1, 1, 32, 192, 288), meta=np.ndarray>
    p0              float32 4B nan
    ps              (lat, lon) float32 221kB dask.array<chunksize=(192, 288), meta=np.ndarray>
    hurs            (member_id, dcpp_init_year, lat, lon) float32 221kB dask.array<chunksize=(1, 1, 192, 288), meta=np.ndarray>
    tas             (member_id, dcpp_init_year, lat, lon) float32 221kB dask.array<chunksize=(1, 1, 192, 288), meta=np.ndarray>
    T_C             (member_id, dcpp_init_year, lat, lon) float32 221kB dask.array<chunksize=(1, 1, 192, 288), meta=np.ndarray>
Attributes: (12/58)
    Conventions:                      CF-1.7 CMIP-6.2
    activity_id:                      CMIP
    branch_method:                    standard
    branch_time_in_child:             674885.0
    branch_time_in_parent:            219000.0
    case_id:                          15
    ...                               ...
    intake_esm_attrs:variable_id:     areacella
    intake_esm_attrs:grid_label:      gn
    intake_esm_attrs:zstore:          gs://cmip6/CMIP6/CMIP/NCAR/CESM2/histor...
    intake_esm_attrs:version:         20190308
    intake_esm_attrs:_data_format_:   zarr
    intake_esm_dataset_key:           CMIP.historical.CESM2.fx.gn
dss_season = (dss
    .groupby('time.season')
    .mean(keep_attrs=True)
             )
dss_season
<xarray.Dataset> Size: 33MB
Dimensions:         (season: 4, member_id: 1, dcpp_init_year: 1, lev: 32,
                     lat: 192, lon: 288, nbnd: 2)
Coordinates:
  * season          (season) object 32B 'DJF' 'JJA' 'MAM' 'SON'
  * member_id       (member_id) object 8B 'r1i1p1f1'
  * dcpp_init_year  (dcpp_init_year) object 8B None
  * lev             (lev) float64 256B -3.643 -7.595 -14.36 ... -976.3 -992.6
  * lat             (lat) float64 2kB -90.0 -89.06 -88.12 ... 88.12 89.06 90.0
  * lon             (lon) float64 2kB -180.0 -178.8 -177.5 ... 176.2 177.5 178.8
    a_bnds          (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
    b_bnds          (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
    lev_bnds        (lev, nbnd) float64 512B dask.array<chunksize=(32, 2), meta=np.ndarray>
Dimensions without coordinates: nbnd
Data variables:
    mmrso4          (season, member_id, dcpp_init_year, lev, lat, lon) float32 28MB dask.array<chunksize=(4, 1, 1, 32, 192, 288), meta=np.ndarray>
    ps              (season, lat, lon) float32 885kB dask.array<chunksize=(4, 192, 288), meta=np.ndarray>
    hurs            (season, member_id, dcpp_init_year, lat, lon) float32 885kB dask.array<chunksize=(4, 1, 1, 192, 288), meta=np.ndarray>
    tas             (season, member_id, dcpp_init_year, lat, lon) float32 885kB dask.array<chunksize=(4, 1, 1, 192, 288), meta=np.ndarray>
    T_C             (season, member_id, dcpp_init_year, lat, lon) float32 885kB dask.array<chunksize=(4, 1, 1, 192, 288), meta=np.ndarray>
    areacella       (season, member_id, dcpp_init_year, lat, lon) float32 885kB dask.array<chunksize=(4, 1, 1, 192, 288), meta=np.ndarray>
    a               (season, lev) float64 1kB dask.array<chunksize=(4, 32), meta=np.ndarray>
    b               (season, lev) float64 1kB dask.array<chunksize=(4, 32), meta=np.ndarray>
    p0              (season) float32 16B 1e+05 1e+05 1e+05 1e+05
Attributes: (12/58)
    Conventions:                      CF-1.7 CMIP-6.2
    activity_id:                      CMIP
    branch_method:                    standard
    branch_time_in_child:             674885.0
    branch_time_in_parent:            219000.0
    case_id:                          15
    ...                               ...
    intake_esm_attrs:variable_id:     areacella
    intake_esm_attrs:grid_label:      gn
    intake_esm_attrs:zstore:          gs://cmip6/CMIP6/CMIP/NCAR/CESM2/histor...
    intake_esm_attrs:version:         20190308
    intake_esm_attrs:_data_format_:   zarr
    intake_esm_dataset_key:           CMIP.historical.CESM2.fx.gn
da_TC = dss_season['T_C']
%%time
da_TC.plot(col='season')
CPU times: user 1.24 s, sys: 294 ms, total: 1.53 s
Wall time: 3.49 s
<xarray.plot.facetgrid.FacetGrid at 0x7f7df031ba10>
../../../../_images/6d29d441ece69da2f4eafef127502ea8f72a6f102a42050627ba6c7418a6d814.png
Note to Sara: Do quiz...

Tip: Might want to load or compute.#

from dask.diagnostics import ProgressBar
%%time
with ProgressBar():
    da_TC.load()
[########################################] | 100% Completed | 3.32 ss
CPU times: user 910 ms, sys: 256 ms, total: 1.17 s
Wall time: 3.33 s

10. Controle the plot visuals:#

%%time
da_TC.plot(
    col = 'season',
    cmap = 'coolwarm',
    xlim = [-100,100],
    ylim = [-0,90],
    cbar_kwargs = {'label':'Temperature [$^\circ$C]'}
    )
CPU times: user 151 ms, sys: 3.23 ms, total: 155 ms
Wall time: 154 ms
<unknown>:6: SyntaxWarning: invalid escape sequence '\c'
<xarray.plot.facetgrid.FacetGrid at 0x7f7dc79df050>
../../../../_images/5af442fb1450552c6b812fbec580567f1ff5f7e1b445b7d2d9a0ba5755f3c932.png
  • Xarray uses certain attributes as labels when plotting:

da_TC = da_TC.assign_attrs({'long_name': 'Temperature near surface'})
da_TC.mean('season', keep_attrs=True).plot()
<matplotlib.collections.QuadMesh at 0x7f7dc7649cd0>
../../../../_images/4a7a52fbcbfa6971e71530c82a4b63aec18a9a114f0aa617752d6ad1d409168f.png

11. Plotting on maps with cartopy#

See more here:

import cartopy as cy
import cartopy.crs as ccrs
da_plt = dss['mmrso4'].isel(lev=-1).mean('time', keep_attrs=True).squeeze()#('member_id')
da_plt
<xarray.DataArray 'mmrso4' (lat: 192, lon: 288)> Size: 221kB
dask.array<getitem, shape=(192, 288), dtype=float32, chunksize=(192, 288), chunktype=numpy.ndarray>
Coordinates:
  * lat             (lat) float64 2kB -90.0 -89.06 -88.12 ... 88.12 89.06 90.0
  * lon             (lon) float64 2kB -180.0 -178.8 -177.5 ... 176.2 177.5 178.8
    member_id       <U8 32B 'r1i1p1f1'
    dcpp_init_year  object 8B None
    lev             float64 8B -992.6
Attributes: (12/19)
    cell_measures:  area: areacella
    cell_methods:   area: time: mean
    comment:        Dry mass of sulfate (SO4) in aerosol particles as a fract...
    description:    Dry mass of sulfate (SO4) in aerosol particles as a fract...
    frequency:      mon
    id:             mmrso4
    ...             ...
    time_label:     time-mean
    time_title:     Temporal mean
    title:          Aerosol Sulfate Mass Mixing Ratio
    type:           real
    units:          kg kg-1
    variable_id:    mmrso4
from matplotlib.colors import LogNorm
with ProgressBar():
    da_plt.load()
[########################################] | 100% Completed | 15.82 s
da_plt
<xarray.DataArray 'mmrso4' (lat: 192, lon: 288)> Size: 221kB
array([[7.0655860e-12, 7.0656250e-12, 7.0656757e-12, ..., 7.0654576e-12,
        7.0654355e-12, 7.0655361e-12],
       [7.4523243e-12, 7.4701807e-12, 7.4873970e-12, ..., 7.3945355e-12,
        7.4097612e-12, 7.4321209e-12],
       [8.0703352e-12, 8.0987569e-12, 8.1331625e-12, ..., 7.9838211e-12,
        8.0185901e-12, 8.0491074e-12],
       ...,
       [3.2195208e-11, 3.2224824e-11, 3.2185154e-11, ..., 3.2204971e-11,
        3.2198522e-11, 3.2200648e-11],
       [3.1399085e-11, 3.1360515e-11, 3.1390588e-11, ..., 3.1466201e-11,
        3.1436742e-11, 3.1405469e-11],
       [3.0563073e-11, 3.0563111e-11, 3.0562844e-11, ..., 3.0563923e-11,
        3.0565089e-11, 3.0564860e-11]], shape=(192, 288), dtype=float32)
Coordinates:
  * lat             (lat) float64 2kB -90.0 -89.06 -88.12 ... 88.12 89.06 90.0
  * lon             (lon) float64 2kB -180.0 -178.8 -177.5 ... 176.2 177.5 178.8
    member_id       <U8 32B 'r1i1p1f1'
    dcpp_init_year  object 8B None
    lev             float64 8B -992.6
Attributes: (12/19)
    cell_measures:  area: areacella
    cell_methods:   area: time: mean
    comment:        Dry mass of sulfate (SO4) in aerosol particles as a fract...
    description:    Dry mass of sulfate (SO4) in aerosol particles as a fract...
    frequency:      mon
    id:             mmrso4
    ...             ...
    time_label:     time-mean
    time_title:     Temporal mean
    title:          Aerosol Sulfate Mass Mixing Ratio
    type:           real
    units:          kg kg-1
    variable_id:    mmrso4
  • projection tells cartopy how we want the data visualized

  • transform tells it what kind of grid the data is on:

    • PlateCarree is a regular lat, lon grid.

Note to Sara: robust, LogNorm, cbar_kwargs, levels
f, ax = plt.subplots(1, 1, dpi=150, figsize = [6,3],
                    subplot_kw={'projection':ccrs.Robinson()},
                    #subplot_kw={'projection':ccrs.Orthographic()},
                   )

da_plt.plot.pcolormesh(
    cmap = 'Reds',
    ax=ax,
    transform = ccrs.PlateCarree(), 
    robust=True,
    norm = LogNorm(), # to add log scale colorbar
    #cbar_kwargs={
    #    #'label':'Wind Speed [m/s]', 
    #    'orientation':'horizontal',
    #    'shrink':.8
    #},
    levels = 6
)
#ax.set_title('Bottom level: Mean over Time')


ax.coastlines(linewidth=.1)
gl = ax.gridlines(draw_labels=True,#{'bottom':True, 'left':True},
                  xlocs=[-180,-140,0,140, 180],
                  ylocs=[-90,-45,0,45,90],
                  x_inline=False,
                  y_inline=False,
                  color='k',
                  linestyle='dotted'
                 )
../../../../_images/de90912d7caa5110b1574f5f7fe258eaaa04c15681e6b04cd081e88ba5f962aa.png

12. Global average: weighted means!#

dss['T_C'].mean().compute()
<xarray.DataArray 'T_C' ()> Size: 4B
array(5.4478564, dtype=float32)
Attributes: (12/19)
    cell_measures:  area: areacella
    cell_methods:   area: time: mean
    comment:        near-surface (usually, 2 meter) air temperature
    description:    near-surface (usually, 2 meter) air temperature
    frequency:      mon
    id:             tas
    ...             ...
    time_label:     time-mean
    time_title:     Temporal mean
    title:          Near-Surface Air Temperature
    type:           real
    units:          $^\circ$ C
    variable_id:    tas

Why is this wrong?

(dss['T_C']
    .weighted(dss['areacella'])
    .mean()
    .compute()
)
<xarray.DataArray 'T_C' ()> Size: 4B
array(14.715101, dtype=float32)
Attributes: (12/19)
    cell_measures:  area: areacella
    cell_methods:   area: time: mean
    comment:        near-surface (usually, 2 meter) air temperature
    description:    near-surface (usually, 2 meter) air temperature
    frequency:      mon
    id:             tas
    ...             ...
    time_label:     time-mean
    time_title:     Temporal mean
    title:          Near-Surface Air Temperature
    type:           real
    units:          $^\circ$ C
    variable_id:    tas
dss['areacella'].plot()
<matplotlib.collections.QuadMesh at 0x7f7ea8e70f50>
../../../../_images/211c28729cead391cf68869e462390b6a1a9ceac17689b390c058114c98c1b8b.png

13 Filtering, grouping etc#

Pick out station:#

lat_kristineberg = 58.24
lon_kristineberg = 11.44
# pick out surface
ds_surf =dss.isel(lev=-1)
ds_kristineberg = ds_surf.sel(lat=lat_kristineberg, lon = lon_kristineberg, method ='nearest')
ds_kristineberg['mmrso4'].plot()
[<matplotlib.lines.Line2D at 0x7f7ea8fac230>]
../../../../_images/76c629faf1785e3307147b236e84054fb29ff6fa168ee07abca69b4fedf6fe79.png

Pick out the vairables we need.

with ProgressBar():
    ds_kristineberg.load()
[########################################] | 100% Completed | 18.64 ss
vl = ['mmrso4','hurs','tas','T_C']
ds_kristineberg[vl]
<xarray.Dataset> Size: 6kB
Dimensions:         (member_id: 1, dcpp_init_year: 1, time: 240)
Coordinates:
  * member_id       (member_id) object 8B 'r1i1p1f1'
  * dcpp_init_year  (dcpp_init_year) object 8B None
  * time            (time) datetime64[ns] 2kB 1990-01-15T12:00:00 ... 2009-12...
    lat             float64 8B 57.96
    lon             float64 8B 11.25
    lev             float64 8B -992.6
Data variables:
    mmrso4          (member_id, dcpp_init_year, time) float32 960B 9.319e-10 ...
    hurs            (member_id, dcpp_init_year, time) float32 960B 85.65 ... ...
    tas             (member_id, dcpp_init_year, time) float32 960B 277.0 ... ...
    T_C             (member_id, dcpp_init_year, time) float32 960B 3.886 ... ...
Attributes: (12/58)
    Conventions:                      CF-1.7 CMIP-6.2
    activity_id:                      CMIP
    branch_method:                    standard
    branch_time_in_child:             674885.0
    branch_time_in_parent:            219000.0
    case_id:                          15
    ...                               ...
    intake_esm_attrs:variable_id:     areacella
    intake_esm_attrs:grid_label:      gn
    intake_esm_attrs:zstore:          gs://cmip6/CMIP6/CMIP/NCAR/CESM2/histor...
    intake_esm_attrs:version:         20190308
    intake_esm_attrs:_data_format_:   zarr
    intake_esm_dataset_key:           CMIP.historical.CESM2.fx.gn

13.1 Resample#

ds_yearly = ds_kristineberg.resample(time='YE').mean('time')
ds_yearly['T_C'].plot()
ds_kristineberg['T_C'].plot(linewidth=0, marker='.')
[<matplotlib.lines.Line2D at 0x7f7dc68b87a0>]
../../../../_images/ff890364572e97b62cc5689c6233e6bcdf95d8a7792b9f74230cd2d72dfafa4a.png

13.2 Filter data using .where():#

ds_kristineberg['season'] = ds_kristineberg['time.season']
for s in ['MAM','JJA','SON','DJF']:
    ds_kristineberg.where(ds_kristineberg['season']==s)['T_C'].plot.hist(alpha=0.5, bins=20, label=s)
    
plt.legend()
<matplotlib.legend.Legend at 0x7f7dc68869c0>
../../../../_images/a10b43065cc2a71ca0e61a659405976cd71b4b7b98240a424930c62ceb8af054.png

13.3 Group by with .groupby()#

ds_kristineberg['month'] = ds_kristineberg['time.month']
ds_kristineberg
<xarray.Dataset> Size: 12kB
Dimensions:         (member_id: 1, dcpp_init_year: 1, time: 240, nbnd: 2)
Coordinates:
  * member_id       (member_id) object 8B 'r1i1p1f1'
  * dcpp_init_year  (dcpp_init_year) object 8B None
  * time            (time) datetime64[ns] 2kB 1990-01-15T12:00:00 ... 2009-12...
    lat             float64 8B 57.96
    lon             float64 8B 11.25
    a_bnds          (nbnd) float64 16B 0.005032 0.002255
    b_bnds          (nbnd) float64 16B 0.0 0.0
    lev             float64 8B -992.6
    lev_bnds        (nbnd) float64 16B 0.005032 0.002255
Dimensions without coordinates: nbnd
Data variables:
    areacella       (member_id, dcpp_init_year) float32 4B 7.728e+09
    a               float64 8B 0.0
    b               float64 8B 0.9926
    mmrso4          (member_id, dcpp_init_year, time) float32 960B 9.319e-10 ...
    p0              float32 4B 1e+05
    ps              (time) float32 960B 1.009e+05 9.98e+04 ... 1.001e+05
    hurs            (member_id, dcpp_init_year, time) float32 960B 85.65 ... ...
    tas             (member_id, dcpp_init_year, time) float32 960B 277.0 ... ...
    T_C             (member_id, dcpp_init_year, time) float32 960B 3.886 ... ...
    season          (time) <U3 3kB 'DJF' 'DJF' 'MAM' 'MAM' ... 'SON' 'SON' 'DJF'
    month           (time) int64 2kB 1 2 3 4 5 6 7 8 9 ... 4 5 6 7 8 9 10 11 12
Attributes: (12/58)
    Conventions:                      CF-1.7 CMIP-6.2
    activity_id:                      CMIP
    branch_method:                    standard
    branch_time_in_child:             674885.0
    branch_time_in_parent:            219000.0
    case_id:                          15
    ...                               ...
    intake_esm_attrs:variable_id:     areacella
    intake_esm_attrs:grid_label:      gn
    intake_esm_attrs:zstore:          gs://cmip6/CMIP6/CMIP/NCAR/CESM2/histor...
    intake_esm_attrs:version:         20190308
    intake_esm_attrs:_data_format_:   zarr
    intake_esm_dataset_key:           CMIP.historical.CESM2.fx.gn
(ds_kristineberg['T_C']
    .groupby(ds_kristineberg['month'])
    .median()
    .plot()
)
[<matplotlib.lines.Line2D at 0x7f7dc504e5d0>]
../../../../_images/d56d3234254d9254d80fd3000f7871cd2aee3f14b2bf5d3e593f15318c303806.png

14. Convert to pandas & do some random fun stuff:#

Maybe we e.g. want to compare with a station, or just use some of the considerable functionalities available from pandas. It’s easy to convert back and forth between xarray and pandas:

A lot of these functions also exist in xarray!

Obs: might want to convert to pandas#

At this point we basically have a time series and the dataset is not to big, so might want to convert to pandas:

subset_of_variables = ['mmrso4','hurs','tas','T_C','season']
df_kristineberg = ds_kristineberg[subset_of_variables].squeeze().to_dataframe()
df_kristineberg

lets do something unnecesarily complicated :D#

Convert to pandas dataframe#

Because more functionality

df = dss.isel(lev=-1)[vl].to_dataframe()
df_ri = df.reset_index()
df_ri.head()

qcut, cut#

qcut splits the data into quantile ranges

df_ri['hurs_cat'] = pd.qcut(df_ri['hurs'],
                            q=[0.05,0.17, 0.34,0.66, 0.83,0.95],
                            labels=['very low','low','med','high','very high'],
                           )

Cut cuts into categories

df_ri['lat_cat'] = pd.cut(df_ri['lat'], [-90,-60,-30,0,30,60,90], 
                          labels=['S polar','S mid','S tropics', 'N tropic', 'N mid','N polar'])
df_ri
df_ri.sample(n=20000).groupby('lat_cat').mean(numeric_only=True)
sns.boxenplot(x="lat_cat", y="hurs",
              color="b",
              #scale="linear", 
              width_method='linear',
              data=df_ri)#.sample(n=20000))
sns.boxenplot(x="hurs_cat", 
              y="mmrso4",
              color="b",
              width_method='linear',
              data=df_ri.sample(n=20000),
             )
sns.displot(x="mmrso4",
            hue='lat_cat',
            log_scale=True,
            kind='kde',
            data=df_ri.sample(n=20000), 
            multiple="stack"
           )

Convert back to xarray if we need:#

ds_new = df_ri.set_index(['time','lat','lon']).to_xarray()
ds_new

mask by category#

ds_new.where(ds_new['hurs_cat']=='very low').mean(['time','lon'])['mmrso4'].plot(label='very low')#vmin = 0, vmax = 1.5e-8)
ds_new.where(ds_new['hurs_cat']=='low').mean(['time','lon'])['mmrso4'].plot(label='low')#vmin = 0, vmax = 1.5e-8)

ds_new.where(ds_new['hurs_cat']=='high').mean(['time','lon'])['mmrso4'].plot(label='high')#vmin = 0, vmax = 1.5e-8)
ds_new.where(ds_new['hurs_cat']=='very high').mean(['time','lon'], keep_attrs=True)['mmrso4'].plot(label='very high')#vmin = 0, vmax = 1.5e-8)
plt.legend()

Improve and reduce risk of typo:#

for cat in ['very low','low','high','very high']:
    ds_new.where(ds_new['hurs_cat']==cat).mean(['time','lon'])['mmrso4'].plot(label=cat) 
plt.legend()