Model evaluation using pyaerocom¶
After this introduction, you will know how to:
- Search for model and observation data in the provided database
- Read model and observation data
- Explore the data (
GriddedDataandUngriddedData) - Colocate data and assess model performance (
ColocatedData) - Convert from pyaerocom to other common libraries such as
pandasorxarray - Store intermediate results as NetCDF that may be used for further analysis.
Pyaerocom API¶
Get started¶
In [1]:
import pyaerocom as pya
import matplotlib.pyplot as plt
pya.change_verbosity('critical', log=pya.const.print_log) # don't output warnings
pya.__version__
Out[1]:
Please make sure to use version 0.8.1.dev1. If you have an earlier version, follow these instructions in order to update your installation.
Change basic data directory to location of pyaerocom-testdata¶
In [2]:
import socket
if socket.gethostname() == 'pc4971':
print('I am on Jonas PC')
DATA_BASEDIR = '/home/jonasg/MyPyaerocom/pyaerocom-testdata'
else:
print('I assume I am on the Abisko Jupyter hub')
DATA_BASEDIR = '/home/notebook/shared-ns1000k/inputs/pyaerocom-testdata/'
In [3]:
pya.const.BASEDIR = DATA_BASEDIR
Searching for available data in pyaerocom-testdata (method browse_database)¶
In [4]:
pya.browse_database('*TM5*')
You can use the Data ID or the Data directory to read this dataset (next step)
Reading of model data (ReadGridded class)¶
In [5]:
DATA_ID = 'TM5_AP3-CTRL2919'
In [6]:
reader = pya.io.ReadGridded(DATA_ID)
print(reader)
In case you are not familiar with the variable names:¶
In [7]:
pya.get_variable('od550aer')
Out[7]:
Read 2010 Aerosol optical depth (at 550 nm) from TM5 model¶
In [8]:
model_data = reader.read_var('od550aer', start=2010)
type(model_data)
Out[8]:
Working with the GriddedData object¶
In [9]:
model_data.ts_type # temporal resolution
Out[9]:
In [10]:
model_data.quickplot_map('June 2010');
Branching off: Converting GriddedData to iris.Cube or xarray.DataArray¶
Actually, the GriddedData object is based on the iris.Cube object which can be accessed via the cube attr:
In [11]:
cube = model_data.cube
In [12]:
type(cube)
Out[12]:
In [13]:
cube
Out[13]:
In [14]:
try: # pyaerocom >= 0.8.1
data_arr = model_data.to_xarray()
except: # pyaerocom < 0.8.1
print('Coming soon in pyaerocom v0.8.1')
import xarray as xarr
data_arr = xarr.DataArray.from_iris(model_data.cube)
data_arr
Out[14]:
Other methods of the GriddedData object that may be useful¶
- model_data.interpolate
- model_data.regrid
- model_data.crop
- model_data.resample_time
- model_data.sel
Extract data at a certain coordinate and plot timeseries¶
In [15]:
LON_ABISKO = 18.8312 #° E
LAT_ABISKO = 68.3495 #° N
subset_abisko = model_data.sel(longitude=LON_ABISKO, latitude=LAT_ABISKO)
try: # pyaerocom >= 0.8.1
subset_abisko.to_xarray().plot();
except: # pyaerocom < 0.8.1
print('Coming soon in pyaerocom v0.8.1')
import xarray as xarr
data_arr = xarr.DataArray.from_iris(subset_abisko.cube).plot()
Reading of ungridded data (ReadUngridded class)¶
- Ungridded data usually comprises timeseries data from different locations on earth and sampled at different times.
- The data is often provided per station in form of text files (e.g. CSV).
Read Aeronet Sun version 3 level 2 data¶
In [16]:
pya.browse_database('Aeronet*V3*Lev2*')
In [17]:
OBS_ID = 'AeronetSunV3Lev2.daily'
In [18]:
obs_reader = pya.io.ReadUngridded(OBS_ID)
print(obs_reader)
Read in all stations (files) that have AOD data¶
In [19]:
obs_data = obs_reader.read(vars_to_retrieve='od550aer')
obs_data
Out[19]:
The UngriddedData object¶
Getting an overview¶
In [20]:
ax = obs_data.plot_station_coordinates(markersize=16)
ax = obs_data.plot_station_coordinates(var_name='od550aer',
start=2008, stop=2012,
markersize=12, color='lime', ax=ax)
Metadata -> available per file that was read in¶
In [21]:
len(obs_data.metadata)
Out[21]:
e.g. metadata of first file:
In [22]:
obs_data.metadata[0]
Out[22]:
Filter data¶
In [23]:
obs_arctic = obs_data.apply_filters(latitude=[70, 90])
obs_arctic.plot_station_coordinates(markersize=40, marker='x', color='r');
Extract data from individual stations (e.g. Lille)¶
In [24]:
lille_data = obs_data.to_station_data('*Lille*')
type(lille_data)
Out[24]:
Quick intro into StationData¶
- Essentially an extended python dictionary with some extras (e.g. you can use
.and[]to access and assign attributes) - contains timeseries data as
pandas.Series - contains metadata
In [25]:
lille_data.keys()
Out[25]:
In [26]:
lille_data.longitude
Out[26]:
In [27]:
lille_data.latitude
Out[27]:
The actual timeseries data can be accessed as pandas.Series simply by using the variable name¶
In [28]:
lille_data.od550aer # THIS is exactly equivalent to the command lille_data['od550aer']
Out[28]:
In [29]:
ax = lille_data.plot_timeseries(var_name='od550aer', marker='x', linestyle='none')
ax = lille_data.plot_timeseries(var_name='od550aer', ts_type='monthly', linestyle='-', lw=3, ax=ax)
Colocation of model and obsdata¶
In [30]:
col_data = pya.colocation.colocate_gridded_ungridded(model_data,
obs_data,
ts_type='monthly',
start=2010)
col_data
Out[30]:
In [31]:
col_data.plot_scatter(marker='o', color='blue', alpha=0.1);
The xarray.DataArray object can be simply accessed via¶
In [32]:
col_data.data
Out[32]:
Save ColocatedData object to NetCDF¶
In [33]:
col_data.savename_aerocom
Out[33]:
In [34]:
col_data.to_netcdf('.')
... and from there, Sara, Diego, Anne et al. have provided you with the perfect introduction.¶
High level colocation and creation of colocated NetCDF files¶
The example above did essentially the following steps:
- Find a model and variable of interest
- Find a corresponding observation network that measures that variable
- Load both model and obsdata for this variable into
GriddedDataandUngriddedData, respectively - Colocate the model to the observation locations and create
ColocatedDataobject - Save the colocated data object as NetCDF file
All these steps can be done with a one-liner in pyaerocom using the Colocator class
Colocate the same model and observation network but now use the Angstrom exponent (ang4487aer) instead of the AOD:
In [35]:
colocator = pya.Colocator(model_id=DATA_ID,
obs_id=OBS_ID,
obs_vars='ang4487aer',
ts_type='monthly',
start=2010,
basedir_coldata='.',
reanalyse_existing=True,
save_coldata=True)
colocator.run()
The Colocator objects stores all ColocatedData objects that were created in it's data attribute, which is a nested dictionary, organised via model_id and var_name:
In [36]:
colocator.data
Out[36]:
In [37]:
colocator.data[DATA_ID]['ang4487aer'].plot_scatter(marker='o', alpha=0.1, color='g');
Do the same for another model (e.g. CAM5.3-Oslo)¶
In [38]:
pya.browse_database('*CAM*')
In [40]:
colocator.run(model_id=ANOTHER_MODEL_ID)
colocator.data.keys()
Out[40]:
In [41]:
colocator.data[ANOTHER_MODEL_ID]['ang4487aer'].plot_scatter(marker='o', alpha=0.1, color='g');
Usage example: Compute AOD bias and plot bias map¶
In [42]:
_arr = col_data.data
mean_bias = ((_arr[1] - _arr[0])/_arr[0]).mean('time') * 100 #%
mean_bias
Out[42]:
In [43]:
ax = pya.plot.mapping.init_map()
_sc = ax.scatter(mean_bias.longitude, mean_bias.latitude, marker='o', c=mean_bias.data, cmap='bwr',
vmin=-100, vmax=100, s=80)
cb = plt.colorbar(_sc)
cb.ax.set_ylabel('Normalised mean bias [%]');
Now if you colocate a lot of models against multiple different observations, you might end up here:¶
https://aerocom-evaluation.met.no/overall.php?project=aerocom&exp=PIII-optics2019