Reading MODIS data with xarray#

Import python packages#

import pandas as pd
import s3fs
import xarray as xr

Connect to bucket (anonymous login for public data only)#

fs = s3fs.S3FileSystem(anon=True,
      client_kwargs={
         'endpoint_url': 'https://climate.uiogeo-apps.sigma2.no/'
      })

fs.ls('ESGF/obs4MIPs/MODIS/MODIS6.1terra')[:10]

['ESGF/obs4MIPs/MODIS/MODIS6.1terra/MakeAllYear.sh',
 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/MakeAllYear_Parallel.sh',
 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_clt_Column_2000_daily.nc',
 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_clt_Column_2001_daily.nc',
 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_clt_Column_2002_daily.nc',
 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_clt_Column_2003_daily.nc',
 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_clt_Column_2004_daily.nc',
 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_clt_Column_2005_daily.nc',
 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_clt_Column_2006_daily.nc',
 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_clt_Column_2007_daily.nc']

s3path = 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/*od550aer*.nc'
remote_files = fs.glob(s3path)
print(remote_files)

['ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_od550aer_Column_2000_daily.nc', 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_od550aer_Column_2001_daily.nc', 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_od550aer_Column_2002_daily.nc', 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_od550aer_Column_2003_daily.nc', 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_od550aer_Column_2004_daily.nc', 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_od550aer_Column_2005_daily.nc', 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_od550aer_Column_2006_daily.nc', 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_od550aer_Column_2007_daily.nc', 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_od550aer_Column_2008_daily.nc', 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_od550aer_Column_2009_daily.nc', 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_od550aer_Column_2010_daily.nc', 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_od550aer_Column_2011_daily.nc', 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_od550aer_Column_2012_daily.nc', 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_od550aer_Column_2013_daily.nc', 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_od550aer_Column_2014_daily.nc', 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_od550aer_Column_2015_daily.nc', 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_od550aer_Column_2016_daily.nc', 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_od550aer_Column_2017_daily.nc', 'ESGF/obs4MIPs/MODIS/MODIS6.1terra/aerocom3_MODIS6.1terra_od550aer_Column_2018_daily.nc']

Access data files#

# Iterate through remote_files to create a fileset
fileset = [fs.open(file) for file in remote_files]

Data reading with xarray#

# Read file with xarray
dset = xr.open_mfdataset(fileset, combine='by_coords')

dset

<xarray.Dataset>
Dimensions:    (latitude: 180, longitude: 360, time: 6940)
Coordinates:
  * latitude   (latitude) float32 89.5 88.5 87.5 86.5 ... -87.5 -88.5 -89.5
  * longitude  (longitude) float32 -179.5 -178.5 -177.5 ... 177.5 178.5 179.5
  * time       (time) datetime64[ns] 2000-01-01 2000-01-02 ... 2018-12-31
Data variables:
    od550aer   (time, latitude, longitude) float64 dask.array<chunksize=(366, 180, 360), meta=np.ndarray>
Attributes:
    HDFEOSVersion:                     HDFEOS_V2.17
    identifier_product_doi:            10.5067/MODIS/MOD08_D3.061
    identifier_product_doi_authority:  http://dx.doi.org
    history:                           Mon Dec  2 09:46:57 2019: ncks -4 -L 5...
    NCO:                               netCDF Operators version 4.7.8 (Homepa...

Plot a single time#

!pip install cmaps

WARNING: The directory '/home/jovyan/.cache/pip' or its parent directory is not owned or is not writable by the current user. The cache has been disabled. Check the permissions and owner of that directory. If executing pip with sudo, you may want sudo's -H flag.
Requirement already satisfied: cmaps in /opt/conda/lib/python3.8/site-packages (1.0.3)
Requirement already satisfied: numpy in /opt/conda/lib/python3.8/site-packages (from cmaps) (1.20.1)
Requirement already satisfied: matplotlib in /opt/conda/lib/python3.8/site-packages (from cmaps) (3.3.4)
Requirement already satisfied: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.3 in /opt/conda/lib/python3.8/site-packages (from matplotlib->cmaps) (2.4.7)
Requirement already satisfied: cycler>=0.10 in /opt/conda/lib/python3.8/site-packages (from matplotlib->cmaps) (0.10.0)
Requirement already satisfied: python-dateutil>=2.1 in /opt/conda/lib/python3.8/site-packages (from matplotlib->cmaps) (2.7.5)
Requirement already satisfied: kiwisolver>=1.0.1 in /opt/conda/lib/python3.8/site-packages (from matplotlib->cmaps) (1.3.1)
Requirement already satisfied: pillow>=6.2.0 in /opt/conda/lib/python3.8/site-packages (from matplotlib->cmaps) (8.1.2)
Requirement already satisfied: six in /opt/conda/lib/python3.8/site-packages (from cycler>=0.10->matplotlib->cmaps) (1.15.0)

import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cmaps

fig=plt.figure(figsize=(20,10))

# We're using cartopy and are plotting in Orthographic projection 
# (see documentation on cartopy)
ax = plt.subplot(1, 1, 1, projection=ccrs.PlateCarree(central_longitude=20.0, globe=None))
ax.coastlines(resolution='110m')

# custom colormap
lcmap = cmaps.BlueYellowRed

# We need to project our data to the new Mercator projection and for this we use `transform`.
# we set the original data projection in transform (here PlateCarree)
# we only plot values greather than 0
dset['od550aer'].sel(time='2001-01-01').plot(ax=ax, transform=ccrs.PlateCarree(), cmap=lcmap, vmin=0., vmax=.5)
ax.set_title('MODIS - 2001-01-01\n ', fontsize=20)

Text(0.5, 1.0, 'MODIS - 2001-01-01\n ')
../../_images/read-MODIS_15_1.png

Interactive plot#

from ipywidgets import interact, interactive, fixed, interact_manual, widgets
from datetime import datetime

def plot_map(date):
    fig=plt.figure(figsize=(20,10))

    ax = plt.subplot(1, 1, 1, projection=ccrs.PlateCarree(central_longitude=20.0, globe=None))
    ax.coastlines(resolution='110m')

    # custom colormap
    lcmap = cmaps.BlueYellowRed

    # We need to project our data to the new Mercator projection and for this we use `transform`.
    # we set the original data projection in transform (here PlateCarree)
    # we only plot values greather than 0
    dset['od550aer'].sel(time=date).plot(ax=ax, transform=ccrs.PlateCarree(), cmap=lcmap, vmin=0., vmax=.5)
    ax.set_title('MODIS - {}\n '.format(date), fontsize=20)

start_date = datetime(2000, 1, 1)
end_date = datetime(2018, 12, 31)

dates = pd.date_range(start_date, end_date, freq='D')

options = [(date.strftime('%Y-%m-%d'), date) for date in dates]
index = (0, len(options)-1)

date_slider = widgets.SelectionSlider(
    options=options,
    orientation='horizontal',
    layout={'width': '800px'}
)

interact(plot_map, date=date_slider);