Reading MODIS data with xarray
Contents
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...
xarray.Dataset
- latitude: 180
- longitude: 360
- time: 6940
- latitude(latitude)float3289.5 88.5 87.5 ... -88.5 -89.5
- standard_name :
- latitude
- units :
- degrees_north
array([ 89.5, 88.5, 87.5, 86.5, 85.5, 84.5, 83.5, 82.5, 81.5, 80.5, 79.5, 78.5, 77.5, 76.5, 75.5, 74.5, 73.5, 72.5, 71.5, 70.5, 69.5, 68.5, 67.5, 66.5, 65.5, 64.5, 63.5, 62.5, 61.5, 60.5, 59.5, 58.5, 57.5, 56.5, 55.5, 54.5, 53.5, 52.5, 51.5, 50.5, 49.5, 48.5, 47.5, 46.5, 45.5, 44.5, 43.5, 42.5, 41.5, 40.5, 39.5, 38.5, 37.5, 36.5, 35.5, 34.5, 33.5, 32.5, 31.5, 30.5, 29.5, 28.5, 27.5, 26.5, 25.5, 24.5, 23.5, 22.5, 21.5, 20.5, 19.5, 18.5, 17.5, 16.5, 15.5, 14.5, 13.5, 12.5, 11.5, 10.5, 9.5, 8.5, 7.5, 6.5, 5.5, 4.5, 3.5, 2.5, 1.5, 0.5, -0.5, -1.5, -2.5, -3.5, -4.5, -5.5, -6.5, -7.5, -8.5, -9.5, -10.5, -11.5, -12.5, -13.5, -14.5, -15.5, -16.5, -17.5, -18.5, -19.5, -20.5, -21.5, -22.5, -23.5, -24.5, -25.5, -26.5, -27.5, -28.5, -29.5, -30.5, -31.5, -32.5, -33.5, -34.5, -35.5, -36.5, -37.5, -38.5, -39.5, -40.5, -41.5, -42.5, -43.5, -44.5, -45.5, -46.5, -47.5, -48.5, -49.5, -50.5, -51.5, -52.5, -53.5, -54.5, -55.5, -56.5, -57.5, -58.5, -59.5, -60.5, -61.5, -62.5, -63.5, -64.5, -65.5, -66.5, -67.5, -68.5, -69.5, -70.5, -71.5, -72.5, -73.5, -74.5, -75.5, -76.5, -77.5, -78.5, -79.5, -80.5, -81.5, -82.5, -83.5, -84.5, -85.5, -86.5, -87.5, -88.5, -89.5], dtype=float32)
- longitude(longitude)float32-179.5 -178.5 ... 178.5 179.5
- standard_name :
- longitude
- units :
- degrees_east
array([-179.5, -178.5, -177.5, ..., 177.5, 178.5, 179.5], dtype=float32)
- time(time)datetime64[ns]2000-01-01 ... 2018-12-31
- long_name :
- time
array(['2000-01-01T00:00:00.000000000', '2000-01-02T00:00:00.000000000', '2000-01-03T00:00:00.000000000', ..., '2018-12-29T00:00:00.000000000', '2018-12-30T00:00:00.000000000', '2018-12-31T00:00:00.000000000'], dtype='datetime64[ns]')
- od550aer(time, latitude, longitude)float64dask.array<chunksize=(366, 180, 360), meta=np.ndarray>
- Aggregation_Data_Set :
- None
- Derived_From_Level_2_Data_Set :
- Optical_Depth_Land_And_Ocean
- Included_Level_2_Nighttime_Data :
- False
- Level_2_Pixel_Values_Read_As :
- Real
- Quality_Assurance_Data_Set :
- None
- Statistic_Type :
- Simple
- long_name :
- Aerosol Optical Thickness at 0.55 microns for both Ocean (best) and Land (corrected): Mean
- units :
- none
- valid_range :
- [-100 5000]
- standard_name :
- atmosphere_optical_thickness_due_to_ambient_aerosol_particles
Array Chunk Bytes 3.60 GB 189.73 MB Shape (6940, 180, 360) (366, 180, 360) Count 57 Tasks 19 Chunks Type float64 numpy.ndarray
- 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 -F -O -v od550aer all.nc aerocom_MODIS6.1terra_od550aer_Column_2000_daily.nc Mon Dec 2 09:46:57 2019: ncatted -O -a long_name,time,o,c,time all.nc Mon Dec 2 09:46:57 2019: ncatted -O -a calendar,time,o,c,standard all.nc Mon Dec 2 09:46:57 2019: ncatted -O -a units,time,o,c,days since 2000-1-1 0:0:0 all.nc Mon Dec 2 09:46:49 2019: ncap2 -o all.nc -O -s time=array(0f,1.,$time) all.nc Mon Dec 2 09:46:34 2019: ncecat -O -u time -n 366,3,1 link_2000_001.nc all.nc Mon Dec 2 09:25:11 2019: ncks -O -o prototype_2000.nc -x -v Aerosol_Optical_Depth_Land_Ocean_Mean,Aerosol_Optical_Depth_Small_Ocean_Mean,Aerosol_Optical_Depth_Small_Ocean_Mean,Aerosol_Optical_Depth_Small_Ocean_Mean prototype_2000.nc Mon Dec 2 09:25:11 2019: ncap2 -o prototype_2000.nc -O -s od660lt1aer(:,:)=-9999. prototype_2000.nc Mon Dec 2 09:25:11 2019: ncap2 -o prototype_2000.nc -O -s od660lt1aer=Aerosol_Optical_Depth_Small_Ocean_Mean(2,:,:) prototype_2000.nc Mon Dec 2 09:25:11 2019: ncap2 -o prototype_2000.nc -O -s od550lt1aer(:,:)=-9999. prototype_2000.nc Mon Dec 2 09:25:10 2019: ncap2 -o prototype_2000.nc -O -s od550lt1aer=Aerosol_Optical_Depth_Small_Ocean_Mean(1,:,:) prototype_2000.nc Mon Dec 2 09:25:10 2019: ncap2 -o prototype_2000.nc -O -s od470lt1aer(:,:)=-9999. prototype_2000.nc Mon Dec 2 09:25:10 2019: ncap2 -o prototype_2000.nc -O -s od470lt1aer=Aerosol_Optical_Depth_Small_Ocean_Mean(0,:,:) prototype_2000.nc Mon Dec 2 09:25:10 2019: ncatted -O -a standard_name,od550aer,o,c,atmosphere_optical_thickness_due_to_ambient_aerosol_particles prototype_2000.nc Mon Dec 2 09:25:09 2019: ncap2 -o prototype_2000.nc -O -s od550aer(:,:)=-9999. prototype_2000.nc Mon Dec 2 09:25:09 2019: ncrename -O -v Aerosol_Optical_Depth_Land_Ocean_Mean,od550aer prototype_2000.nc Mon Dec 2 09:25:09 2019: ncatted -O -a units,longitude,o,c,degrees_east prototype_2000.nc Mon Dec 2 09:25:09 2019: ncatted -O -a standard_name,longitude,o,c,longitude prototype_2000.nc Mon Dec 2 09:25:09 2019: ncap2 -o prototype_2000.nc -O -s longitude=array(-179.5f,1.,$longitude) prototype_2000.nc Mon Dec 2 09:25:09 2019: ncatted -O -a units,latitude,o,c,degrees_north prototype_2000.nc Mon Dec 2 09:25:08 2019: ncatted -O -a standard_name,latitude,o,c,latitude prototype_2000.nc Mon Dec 2 09:25:08 2019: ncap2 -o prototype_2000.nc -O -s latitude=array(89.5f,-1.,$latitude) prototype_2000.nc Mon Dec 2 09:25:08 2019: ncrename -O -d YDim:mod08,latitude -d XDim:mod08,longitude prototype_2000.nc Mon Dec 2 09:25:08 2019: ncpdq -O -U prototype_2000.nc prototype_2000.nc Mon Dec 2 09:25:08 2019: ncks -O -v Aerosol_Optical_Depth_Land_Ocean_Mean,Aerosol_Optical_Depth_Small_Ocean_Mean,Aerosol_Optical_Depth_Small_Ocean_Mean,Aerosol_Optical_Depth_Small_Ocean_Mean ../download/MOD08_D3.A2000055.061.2017276160246.hdf prototype_2000.nc
- NCO :
- netCDF Operators version 4.7.8 (Homepage = http://nco.sf.net, Code = http://github.com/nco/nco)
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 ')
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);