In [1]:
# Ensure I don't use any local plugins. Set it to a readable folder with no Python files to avoid warnings.
%env CIS_PLUGIN_HOME=/Users/watson-parris/Pictures


env: CIS_PLUGIN_HOME=/Users/watson-parris/Pictures

CIS Introduction

CIS has it's own version of the Iris Cube. But it's designed to work with any observational data. The CIS data structure is just called UngriddedData:

First unzip your example data to a folder you can easily find it.


In [2]:
from cis import read_data, read_data_list, get_variables

get_variables('../resources/WorkshopData2016/Aeronet/920801_150530_Brussels.lev20')


Out[2]:
{'340-440Angstrom',
 '380-500Angstrom',
 '440-675Angstrom',
 '440-675Angstrom(Polar)',
 '440-870Angstrom',
 '500-870Angstrom',
 'AOT_1020',
 'AOT_1640',
 'AOT_340',
 'AOT_380',
 'AOT_412',
 'AOT_440',
 'AOT_443',
 'AOT_490',
 'AOT_500',
 'AOT_531',
 'AOT_532',
 'AOT_551',
 'AOT_555',
 'AOT_667',
 'AOT_675',
 'AOT_870',
 'Date(dd-mm-yy)',
 'Julian_Day',
 'Last_Processing_Date(ddmmyyyy)',
 'Solar_Zenith_Angle',
 'Time(hh:mm:ss)',
 'TripletVar_1020',
 'TripletVar_1640',
 'TripletVar_340',
 'TripletVar_380',
 'TripletVar_412',
 'TripletVar_440',
 'TripletVar_443',
 'TripletVar_490',
 'TripletVar_500',
 'TripletVar_531',
 'TripletVar_532',
 'TripletVar_551',
 'TripletVar_555',
 'TripletVar_667',
 'TripletVar_675',
 'TripletVar_870',
 'Water(cm)',
 'WaterError'}

In [3]:
aeronet_aot_500 = read_data("../resources/WorkshopData2016/Aeronet/920801_150530_Brussels.lev20", "AOT_500")
print(aeronet_aot_500)


Ungridded data: AOT_500 / (1) 
     Shape = (15573,)

     Total number of points = 15573
     Number of non-masked points = 15573
     Long name = AOT_500
     Standard name = None
     Units = 1
     Missing value = -999.0
     Range = (0.01864, 1.4709920000000001)
     History = 
     Coordinates: 
       Longitude
          Long name = 
          Standard name = longitude
          Units = degrees_east
          Missing value = None
          Range = (4.3499999999999996, 4.3499999999999996)
          History = 
       Latitude
          Long name = 
          Standard name = latitude
          Units = degrees_north
          Missing value = None
          Range = (50.783000000000001, 50.783000000000001)
          History = 
       Altitude
          Long name = 
          Standard name = altitude
          Units = meters
          Missing value = None
          Range = (120.0, 120.0)
          History = 
       DateTime
          Long name = 
          Standard name = time
          Units = days since 1600-01-01 00:00:00
          Missing value = None
          Range = (2006-07-13 08:03:47, 2014-05-25 08:55:29)
          History = 


In [4]:
aeronet_aot_500.name()


Out[4]:
'AOT_500'

Some example datasets


In [5]:
%matplotlib inline

Ungridded time series data


In [6]:
aeronet_aot_500.plot()


Out[6]:
<matplotlib.axes._subplots.AxesSubplot at 0x113f131d0>

In [7]:
ax = aeronet_aot_500.plot(color='red')
ax.set_yscale('log')



In [8]:
aeronet_aot = read_data_list("../resources/WorkshopData2016/Aeronet/920801_150530_Brussels.lev20", 
                             ['AOT_500', 'AOT_675'])
ax = aeronet_aot.plot()



In [9]:
ax.set_title('Brussels Aeronet AOT')
ax.set_xlabel('Date')


Out[9]:
<matplotlib.text.Text at 0x1138b6ef0>

In [10]:
from datetime import datetime 
ax.set_xlim(datetime(2007,5,5), datetime(2007,8,26))


Out[10]:
(732801.0, 732914.0)

In [11]:
aeronet_aot.plot(how='comparativescatter')
# Note that this will only work if we have two datasets in our list


Out[11]:
<matplotlib.axes._subplots.AxesSubplot at 0x104a87c88>

Subsetting

CIS is able to subset datasets across any of the given coordinates


In [12]:
aeronet_aot_2007 = aeronet_aot.subset(t=[datetime(2007,1,1), datetime(2007,12,31)])
aeronet_aot_2007


Out[12]:
[<cis 'UngriddedData' of Ungridded data: AOT_500 / (1) 
>,
<cis 'UngriddedData' of Ungridded data: AOT_675 / (1) 
>]

In [13]:
aeronet_aot_2007.plot(how='comparativescatter')


Out[13]:
<matplotlib.axes._subplots.AxesSubplot at 0x115022dd8>

Model time series


In [14]:
model_aod = read_data("../resources/WorkshopData2016/od550aer.nc", "od550aer")

In [15]:
print(model_aod)


od550aer / (1)                      (time: 1464; latitude: 96; longitude: 192)
     Dimension coordinates:
          time                           x               -              -
          latitude                       -               x              -
          longitude                      -               -              x
     Attributes:
          CDI: Climate Data Interface version 1.6.9 (http://mpimet.mpg.de/cdi)
          CDO: Climate Data Operators version 1.6.9 (http://mpimet.mpg.de/cdo)
          Conventions: CF-1.0
          advection: Lin & Rood
          date_time: 20140613 051140
          echam_version: 6.1.00
          grid_type: gaussian
          history: Fri Dec 04 17:15:40 2015: cdo -r copy od550aer.nc ../../temp.nc
Tue Jul...
          host_name: p101
          institution: Max-Planck-Institute for Meteorology
          jsbach_version: 2.01
          operating_system: AIX 6.1 Power6
          physics: Modified ECMWF physics
          radiation: Modified ECMWF radiation
          title: AEROCOM
          truncation: 63
          user_name: Andreas Veira (m300259)

In [16]:
import iris.analysis
maod_global_mean, = model_aod.collapsed(['longitude', 'latitude'], iris.analysis.MEAN)


WARNING:root:Creating guessed bounds as none exist in file
WARNING:root:Creating guessed bounds as none exist in file
WARNING:root:Creating guessed bounds as none exist in file
/Users/watson-parris/anaconda/envs/python_workshop/lib/python3.5/site-packages/iris/analysis/cartography.py:376: UserWarning: Using DEFAULT_SPHERICAL_EARTH_RADIUS.
  warnings.warn("Using DEFAULT_SPHERICAL_EARTH_RADIUS.")

In [17]:
print(maod_global_mean)


od550aer / (1)                      (time: 1464)
     Dimension coordinates:
          time                           x
     Scalar coordinates:
          latitude: 0.0 degrees, bound=(-89.4969872937, 89.4969872937) degrees
          longitude: 179.0625 degrees, bound=(-0.9375, 359.0625) degrees
     Attributes:
          CDI: Climate Data Interface version 1.6.9 (http://mpimet.mpg.de/cdi)
          CDO: Climate Data Operators version 1.6.9 (http://mpimet.mpg.de/cdo)
          Conventions: CF-1.0
          advection: Lin & Rood
          date_time: 20140613 051140
          echam_version: 6.1.00
          grid_type: gaussian
          history: Fri Dec 04 17:15:40 2015: cdo -r copy od550aer.nc ../../temp.nc
Tue Jul...
          host_name: p101
          institution: Max-Planck-Institute for Meteorology
          jsbach_version: 2.01
          operating_system: AIX 6.1 Power6
          physics: Modified ECMWF physics
          radiation: Modified ECMWF radiation
          title: AEROCOM
          truncation: 63
          user_name: Andreas Veira (m300259)
     Cell methods:
          mean: longitude, latitude

In [18]:
ax = maod_global_mean.plot(itemwidth=2)



In [19]:
aeronet_aot_500.plot(ax=ax)


Out[19]:
<matplotlib.axes._subplots.AxesSubplot at 0x118c129b0>

Aircraft data


In [20]:
number_concentration = read_data('../resources/WorkshopData2016/ARCPAC_2008', 
                                 'NUMBER_CONCENTRATION')
print(number_concentration)


WARNING:root:Identified 1 point(s) which were missing values for some or all coordinates - these points have been removed from the data.
Ungridded data: NUMBER_CONCENTRATION / (#/cm3) 
     Shape = (1268,)

     Total number of points = 1268
     Number of non-masked points = 981
     Long name = 
     Standard name = None
     Units = #/cm3
     Missing value = -9999
     Range = (17.0, 2191.0)
     History = 
     Misc attributes: 
       Missing_Value = -9999
     Coordinates: 
       TIME
          Long name = 
          Standard name = time
          Units = days since 1600-01-01 00:00:00
          Missing value = -9999
          Range = (2008-04-11 20:01:30, 2008-04-16 03:03:30)
          History = 
          Misc attributes: 
            Missing_Value = -9999
       LATITUDE
          Long name = 
          Standard name = latitude
          Units = degrees
          Missing value = -9999
          Range = (63.303199999999997, 75.128600000000006)
          History = 
          Misc attributes: 
            Missing_Value = -9999
       LONGITUDE
          Long name = 
          Standard name = longitude
          Units = degrees
          Missing value = -9999
          Range = (-165.1241, -139.23159999999999)
          History = 
          Misc attributes: 
            Missing_Value = -9999
       ALTITUDE
          Long name = 
          Standard name = altitude
          Units = m
          Missing value = -9999
          Range = (77.200000000000003, 6579.6000000000004)
          History = 
          Misc attributes: 
            Missing_Value = -9999


In [21]:
ax = number_concentration.plot()



In [22]:
ax.bluemarble()

Satellite data


In [23]:
aerosol_cci = read_data('../resources/WorkshopData2016/AerosolCCI', 
                        'AOD550')
aerosol_cci.plot()


Out[23]:
<cartopy.mpl.geoaxes.GeoAxesSubplot at 0x1137de9b0>

In [24]:
aerosol_cci_one_day = read_data('../resources/WorkshopData2016/AerosolCCI/20080415*.nc', 
                                'AOD550')
ax = aerosol_cci_one_day.plot()



In [25]:
aerosol_cci_one_day.plot(projection='Orthographic')


Out[25]:
<cartopy.mpl.geoaxes.GeoAxesSubplot at 0x118a7a240>

In [26]:
ax=aerosol_cci_one_day.plot(projection='InterruptedGoodeHomolosine')
ax.bluemarble()


Aggregation

Given a set of UngriddedData...

... we can perform an aggregation over a specified grid...

... to create a new GriddedData object (which is essentiall an Iris Cube)


In [27]:
gridded_aerosol_cci_one_day = aerosol_cci_one_day.aggregate(x=[-180,180,10], y=[-90,90,5])

In [28]:
gridded_aerosol_cci_one_day[0].plot()


Out[28]:
<cartopy.mpl.geoaxes.GeoAxesSubplot at 0x116f9e5f8>

Exercises

1. Read in AOD550 and AOD670 from the 5 days of satellite data


In [ ]:

2. Subset this data down to the region covered by the aircraft data


In [ ]:

3. Try plotting AOD550 against AOD670 from the subsetted satellite data using a comparative scatter plot


In [ ]:

Collocation

Model onto Aeronet

This is an gridded onto un-gridded collocation and can be done using either linear interpolation or nearest neighbour.

This is very quick and in general CIS can even handle hybrid height coordinates:


In [29]:
# Lets take a closer look at the model data
print(model_aod)


od550aer / (1)                      (time: 1464; latitude: 96; longitude: 192)
     Dimension coordinates:
          time                           x               -              -
          latitude                       -               x              -
          longitude                      -               -              x
     Attributes:
          CDI: Climate Data Interface version 1.6.9 (http://mpimet.mpg.de/cdi)
          CDO: Climate Data Operators version 1.6.9 (http://mpimet.mpg.de/cdo)
          Conventions: CF-1.0
          advection: Lin & Rood
          date_time: 20140613 051140
          echam_version: 6.1.00
          grid_type: gaussian
          history: Fri Dec 04 17:15:40 2015: cdo -r copy od550aer.nc ../../temp.nc
Tue Jul...
          host_name: p101
          institution: Max-Planck-Institute for Meteorology
          jsbach_version: 2.01
          operating_system: AIX 6.1 Power6
          physics: Modified ECMWF physics
          radiation: Modified ECMWF radiation
          title: AEROCOM
          truncation: 63
          user_name: Andreas Veira (m300259)

In [30]:
from cis.time_util import PartialDateTime
# First subset the aeronet data:
aeronet_aot_2008 = aeronet_aot_500.subset(t=PartialDateTime(2008))

Note that we don’t actually have to do this subsetting, but that otherwise CIS will interpolate the nearest values, which in this case we don’t really want.


In [31]:
# Now do the collocation:
model_aod_onto_aeronet = model_aod.collocated_onto(aeronet_aot_2008)


/Users/watson-parris/anaconda/envs/python_workshop/lib/python3.5/site-packages/cis/collocation/gridded_interpolation.py:181: MaskedArrayFutureWarning: setting an item on a masked array which has a shared mask will not copy the mask and also change the original mask array in the future.
Check the NumPy 1.11 release notes for more information.
  expanded_result[~self.missing_mask] = result

In [32]:
print(model_aod_onto_aeronet[0])


Ungridded data: od550aer / (1) 
     Shape = (1885,)

     Total number of points = 1885
     Number of non-masked points = 1885
     Long name = Optical thickness - total 550nm
     Standard name = None
     Units = 1
     Missing value = nan
     Range = (0.02192891966422272, 1.6932252635449034)
     History = 2016-12-07T11:53:21Z Collocated onto sampling from: [] 
using CIS version 1.5.1 
variables: od550aer 
with files: ['../resources/WorkshopData2016/od550aer.nc'] 
using collocator: <cis.collocation.col_implementations.GriddedUngriddedCollocator object at 0x11041bf60> 
kernel: lin
     Coordinates: 
       Longitude
          Long name = 
          Standard name = longitude
          Units = degrees_east
          Missing value = None
          Range = (4.3500000000000227, 4.3500000000000227)
          History = 
       Latitude
          Long name = 
          Standard name = latitude
          Units = degrees_north
          Missing value = None
          Range = (50.783000000000001, 50.783000000000001)
          History = 
       Altitude
          Long name = 
          Standard name = altitude
          Units = meters
          Missing value = None
          Range = (120.0, 120.0)
          History = 
       DateTime
          Long name = 
          Standard name = time
          Units = days since 1600-01-01 00:00:00
          Missing value = None
          Range = (2008-03-03 14:25:42, 2008-12-30 13:46:04)
          History = 

Note the updated history


In [33]:
from cis.plotting.plot import multilayer_plot, taylor_plot
ax = multilayer_plot([model_aod_onto_aeronet[0], aeronet_aot_2008], 
                     layer_opts=[dict(label='Model'), 
                                 dict(label='Aeronet')], xaxis='time',
                    itemwidth=1)