This post aims to shed light on how to download on python environment various climatic datasets from the Climate Data Store (CDS). This includes ERA5, CMIP5, CRU+GPCC bias corrected data and others.
Installation of CDS API¶
In [1]:
#pip install cdsapi
The following packages are needed for running the next instructions.
In [2]:
import cdsapi
import xarray as xr
import pandas as pd
import numpy as np
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.geoaxes import GeoAxes
from mpl_toolkits.axes_grid1 import AxesGrid
import matplotlib.pyplot as plt
import seaborn as sns
import tarfile
%matplotlib inline
Download monthly average temperature 1979-2020 from ERA5¶
ERA 5 is the recent climate reanalysis produced by ECMWF, providing hourly data on many atmospheric, land-surface and sea-state parameters together with estimates of uncertainty. ERA5 combines vast amounts of historical observations into global estimates using advanced modelling and data assimilation systems. It provides various weather parmeters from 1979 to present. ERA 5 is widely used as alternative in areas with low or limited number of weather observations. More information can be found from: https://climate.copernicus.eu/climate-reanalysis
In this example, the datasets will be downloaded over West-Africa (20W, 20E,0,25N) from 1979-2020. Depending of the internet connection the operation below can take some few minutes.
The request code below to download the ERA5 data is automatically generated from CDS website based on the datasets choosed, variables, spatial and temporal scale. Registration to access freely the data is required (https://cds.climate.copernicus.eu/user/register).
In [3]:
import cdsapi
c = cdsapi.Client()
c.retrieve(
'reanalysis-era5-single-levels-monthly-means',
{
'format': 'netcdf',
'product_type': 'monthly_averaged_reanalysis',
'variable': '2m_temperature',
'year': [
'1979', '1980', '1981',
'1982', '1983', '1984',
'1985', '1986', '1987',
'1988', '1989', '1990',
'1991', '1992', '1993',
'1994', '1995', '1996',
'1997', '1998', '1999',
'2000', '2001', '2002',
'2003', '2004', '2005',
'2006', '2007', '2008',
'2009', '2010', '2011',
'2012', '2013', '2014',
'2015', '2016', '2017',
'2018', '2019', '2020',
],
'month': [
'01', '02', '03',
'04', '05', '06',
'07', '08', '09',
'10', '11', '12',
],
'time': '00:00',
'area': [
25, -20, 0,
25,
],
},
'Monthly_averaged_temp_ERA5_1979_2020.nc')
Out[3]:
Download CMIP5 monthly outputs precipitation flux¶
In [5]:
import cdsapi
c = cdsapi.Client()
c.retrieve(
'projections-cmip5-monthly-single-levels',
{
'experiment': 'historical',
'variable': 'mean_precipitation_flux',
'model': 'gfdl_cm2p1',
'ensemble_member': 'r10i1p1',
'period': '199601-200012',
'format': 'tgz',
},
'preflux_cmip5_gfdl_cm2p1_1996_2000.tar.gz')
tf = tarfile.open("preflux_cmip5_gfdl_cm2p1_1996_2000.tar.gz")
tf.extractall()
Hourly rainfall flux from Near surface meteorological variables from 1979 to 2018 derived from bias-corrected reanalysis datasets¶
In [6]:
import cdsapi
c = cdsapi.Client()
c.retrieve(
'derived-near-surface-meteorological-variables',
{
'variable': 'rainfall_flux',
'reference_dataset': 'cru_and_gpcc',
'year': '2016',
'month': '08',
'format': 'tgz',
},
'hourly_rainfall_flux_cru_gpcc_2016_08.tar.gz')
tf = tarfile.open("hourly_rainfall_flux_cru_gpcc_2016_08.tar.gz")
tf.extractall()
Import datasets¶
In [5]:
# import datasets
df_temp= xr.open_dataset('Monthly_averaged_temp_ERA5_1979_2020.nc')
df_prcp= xr.open_dataset('pr_Amon_GFDL-CM2p1_historical_r10i1p1_199601-200012.nc')
In [6]:
df_temp
Out[6]:
In [7]:
df_prcp.pr
Out[7]:
In [8]:
# Converts Temperature data into Degree celsius
# Converts Convect
df_temp_deg=df_temp.t2m-273.15
df_temp_deg.attrs = df_temp.t2m.attrs
df_temp_deg["units"] = "deg C"
df_temp_deg
# Converts Precipitation flux data into mm/month
df_pre=df_prcp.pr*86400*30
Monthly variation¶
In [9]:
#Monthly mean
monmean_temp=df_temp_deg.groupby('time.month').mean().mean('expver')
monmean_prcp=df_pre.groupby('time.month').mean()
In [10]:
monmean_temp
Out[10]:
In [11]:
monmean_prcp
Out[11]:
In [12]:
projection = ccrs.PlateCarree()
axes_class = (GeoAxes,dict(map_projection=projection))
longitude=monmean_temp.longitude.values #longitude
latitude=monmean_temp.latitude.values #latitude
month=monmean_temp.month.values
months=['JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG','SEP','OCT','NOV','DEC'] #
fig = plt.figure()
fig.set_size_inches(18.5, 10.5)
axgrid = AxesGrid(fig, 111, axes_class=axes_class,
nrows_ncols=(3, 4),
axes_pad=0.3,
cbar_location='right',
cbar_mode='single',
cbar_pad=0.2,
cbar_size='3%',
label_mode='')
for n, ax in enumerate(axgrid):
ax.set_xlim(-20,25)
ax.set_ylim(0,25)
ax.set_xticks(np.linspace(-20,25,4), crs=projection)
ax.set_yticks(np.linspace(0,25, 4), crs=projection)
ax.coastlines()
p = ax.contourf(longitude, latitude, monmean_temp[n],
transform=projection,
cmap='RdBu_r')
ax.text(20,20, str(months[n]),fontsize=18, ha='center')
d=axgrid.cbar_axes[0].colorbar(p)
d.set_label_text('Deg C')
# Make ticklabels on inner axes invisible
axes = np.reshape(axgrid, axgrid.get_geometry())
for ax in axes[:-1, :].flatten():
ax.xaxis.set_tick_params(which='both',
labelbottom=False, labeltop=False)
for ax in axes[:, 1:].flatten():
ax.yaxis.set_tick_params(which='both',
labelbottom=False, labeltop=False)
#plt.suptitle('Monthly Temperatuve over West-Africa',fontsize=20,y=1.04)
plt.show()
Temperature standardised anomaly with base climatology (1979-2008)¶
In [34]:
interan_temp=df_temp.groupby('time.year').mean().mean('expver') #Annual mean
df = interan_temp.t2m.to_dataframe()
dd=df.t2m.groupby('year').mean() # Yearly mean
clim_temp=df.loc[1979:2008,].t2m.groupby('year').mean().mean() # climatology mean
clim_sd=df.loc[1979:2008,].t2m.groupby('year').mean().std()# climatology standard deviation
temp_anom=(dd-clim_temp)/clim_sd # anomaly compuatation
In [35]:
fig = plt.figure(figsize=(10,5))
cols=['red' if x <0 else 'blue' for x in temp_anom]
plt.bar(temp_anom.index[1:41],temp_anom[1:41],color=cols)
plt.xlabel('Years'),plt.ylabel('Anomaly')
plt.title('Temperature standardised anomaly', loc='left', fontsize=15)
plt.title('1979-2020',loc='right', fontsize=15)
plt.show()
Interannual variability of precipitation per climatic Zones¶
West African domain is divided into three climatic zones: Guinea (4°N–8°N), Savanna (8°N–11°N), and Sahel (11°N–16°N).
In [36]:
#Annual mean
interan_prcp=df_pre.groupby('time.year').sum()
# extract datasets coordinate extent (longitude,latitude)
lats=interan_prcp.lat.values
lonts=interan_prcp.lon.values
In [37]:
df_pre['time']=pd.date_range('1/1/1996', '31/12/2000', freq='M') #fix time series dimensions
In [38]:
#Mean across climatic zone
lat_guinne=lats[(lats>4) & (lats<8)] #Guinea zone
lat_savanna=lats[(lats>8) & (lats<11)] # Savanna zone
lat_sahel=lats[(lats>11) & (lats<16)] # Sahel Zone
In [39]:
guinee=df_pre.sel(lat=lat_guinne).mean('lat').mean('lon').to_dataframe().reset_index()
savanna=df_pre.sel(lat=lat_savanna).mean('lat').mean('lon').to_dataframe().reset_index()
sahel=df_pre.sel(lat=lat_sahel).mean('lat').mean('lon').to_dataframe().reset_index()
In [49]:
#Ploting
fig, ax= plt.subplots(1, 3, figsize=(20,7),sharey=True)
ax[0].plot(guinee.time,guinee.pr,color='green',marker='p',markersize=10)
ax[0].legend(['Guinee'],loc='upper right', frameon=True)
ax[0].set_xlabel('Years')
ax[0].set_ylabel('mm/momth',fontsize=15)
ax[1].plot(savanna.time,savanna.pr,marker='p',color='orange',markersize=10)
ax[1].legend(['Savanna'],loc='upper right', frameon=True)
ax[2].plot(sahel.time,sahel.pr,marker='p',color='blue',markersize=10)
ax[2].legend(['Sahel'],loc='upper right')
plt.axis(aspect='equal')
plt.suptitle('Precipitation per climatic zone',y=1.04,fontsize=20)
plt.tight_layout()
plt.show()
Daily variation of rainfall in August across climatic zones in West-Africa¶
In [29]:
CRU_GPCC_RAIN=xr.open_dataset('Rainf_WFDE5_CRU+GPCC_201608_v1.0.nc')
CRU_GPCC_RAIN
Out[29]:
In [30]:
# Latitude and Longitude extent of the study area
lat_slice = slice(5,30)
lon_slice = slice(-15,30)
# Extract the datasets over West-Africa
Rain_new=CRU_GPCC_RAIN['Rainf'].sel(lat=lat_slice,lon=lon_slice)
dat=pd.to_datetime(Rain_new.time.values) #retrieved 'time' dimensions in datetime format
ts=dat.strftime('%Y-%m-%d') #convert from datetime to time dimension
Rain_new['time']=ts # force 'time' dimension with the daily time series
Rain_new
Out[30]:
In [31]:
lats=Rain_new.lat.values #get latitude extent of the datasets
lat_guinne=lats[(lats>4) & (lats<8)] #retrieved Guinea zone latitude
lat_savanna=lats[(lats>8) & (lats<11)] #retrieved Savanna zone latitude
lat_sahel=lats[(lats>11) & (lats<16)] #retrieved Guinea zone latitude
In [32]:
#Mean across climatic zone
guinee=Rain_new.sel(lat=lat_guinne).mean('lat').mean('lon').to_dataframe().reset_index()
savanna=Rain_new.sel(lat=lat_savanna).mean('lat').mean('lon').to_dataframe().reset_index()
sahel=Rain_new.sel(lat=lat_sahel).mean('lat').mean('lon').to_dataframe().reset_index()
#daily mean of precipiatation in mm/day
interan_guinee=guinee.groupby('time').sum()*86400
interan_savan=savanna.groupby('time').sum()*86400
interan_sahel=sahel.groupby('time').sum()*86400
In [33]:
#Ploting
fig = plt.figure(figsize=(10,5))
ax = plt.axes()
ax.plot(interan_guinee, '--',marker='o')
ax.plot(interan_savan, '--',marker='o')
ax.plot(interan_sahel, '--',marker='o')
ax.set_xticklabels(range(1,32))
ax.set_xlabel('Years')
ax.set_ylabel('mm')
ax.legend(['Guinee','Savanna','Sahel'],loc='upper right', frameon=True)
ax.set_title('Precipitation per climatic zone - daily variation', loc='left',Fontsize=15)
ax.set_title('Time: 2016',loc='right', fontsize=10)
plt.tight_layout()
plt.show()
Thank you very much
ReplyDeleteVery nice bro.
ReplyDeleteGreat Job!!
ReplyDeleteThanks for the blog post buddy! Keep them coming... teknisk analys program
ReplyDeleteThank you so much. I would world to discuss more with you about how you build your site as well.
Delete