Zonal and Anomaly variation over West-Africa
This post shed light on how to produce Hovmöller diagrams in R software and compute anomaly on basics weather parameters.
Load library
library("raster")
library("rasterVis")
library("maptools")
library("rgdal")
library("gridExtra")Datasets
Temperature and Precipitation are from Climatic Research Unit (CRU) Time-Series (TS) version 4.04 and can be downloaded from the link 1 and link 2 respectively.
Import datasets
# Create a RasterStack object
precip <- stack('cru_ts4.04.1901.2019.pre.dat.nc')## Loading required namespace: ncdf4tmp <- stack('cru_ts4.04.1901.2019.tmp.dat.nc')
precip## class : RasterStack
## dimensions : 360, 720, 259200, 1428 (nrow, ncol, ncell, nlayers)
## resolution : 0.5, 0.5 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : X1901.01.16, X1901.02.15, X1901.03.16, X1901.04.16, X1901.05.16, X1901.06.16, X1901.07.16, X1901.08.16, X1901.09.16, X1901.10.16, X1901.11.16, X1901.12.16, X1902.01.16, X1902.02.15, X1902.03.16, ...tmp## class : RasterStack
## dimensions : 360, 720, 259200, 1428 (nrow, ncol, ncell, nlayers)
## resolution : 0.5, 0.5 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : X1901.01.16, X1901.02.15, X1901.03.16, X1901.04.16, X1901.05.16, X1901.06.16, X1901.07.16, X1901.08.16, X1901.09.16, X1901.10.16, X1901.11.16, X1901.12.16, X1902.01.16, X1902.02.15, X1902.03.16, ...Set time dimension
# Create time series from the dataset and set Z layer
dates <- seq(as.Date('1901-01-16'), as.Date('2019-12-16'), 'month')
tmp_new <- setZ(tmp,dates)
precip_new <- setZ(precip,dates)
tmp_new## class : RasterStack
## dimensions : 360, 720, 259200, 1428 (nrow, ncol, ncell, nlayers)
## resolution : 0.5, 0.5 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : X1901.01.16, X1901.02.15, X1901.03.16, X1901.04.16, X1901.05.16, X1901.06.16, X1901.07.16, X1901.08.16, X1901.09.16, X1901.10.16, X1901.11.16, X1901.12.16, X1902.01.16, X1902.02.15, X1902.03.16, ...
## time : 1901-01-16 - 2019-12-16 (range)precip_new## class : RasterStack
## dimensions : 360, 720, 259200, 1428 (nrow, ncol, ncell, nlayers)
## resolution : 0.5, 0.5 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : X1901.01.16, X1901.02.15, X1901.03.16, X1901.04.16, X1901.05.16, X1901.06.16, X1901.07.16, X1901.08.16, X1901.09.16, X1901.10.16, X1901.11.16, X1901.12.16, X1902.01.16, X1902.02.15, X1902.03.16, ...
## time : 1901-01-16 - 2019-12-16 (range)Geographic crop of datasets across West-Africa
WA_temp <- crop(tmp_new,extent(-20,30,0,30))
WA_rainfall <- crop(precip_new,extent(-20,30,0,30))Monthly mean from 1901 to 2019
# Retrieve month
months <- as.integer(format(dates, "%m"))
#Apply mean fuction on the monthly dataset
mon_mean_temp<- stackApply(WA_temp, months, fun=mean)
mon_mean_rainfall<- stackApply(WA_rainfall, months, fun=mean)
#Set time dimension
mon_mean_temp<- setZ(mon_mean_temp, month.abb)
mon_mean_rainfall<- setZ(mon_mean_rainfall, month.abb)
names(mon_mean_temp)<-month.abb
names(mon_mean_rainfall)<-month.abbMonthly average of across the dataset
Define color palette
#define color palette for each type of dataset
panel_temp=c("#FFFFC8","#FFF8BD" ,"#FDEDAB" ,"#FBE093" ,"#F9D277", "#F7C252", "#F6B00B","#F49D00" ,
"#F28800","#EF7000","#EB5500", "#DF3A00" ,"#CA2700" ,"#B21400","#980217","#7D0025")
panel_prcp=c("#F9FFAF" ,"#E0F6AB" ,"#C7ECA7", "#AFE2A3", "#96D7A1", "#7FCD9F"
,"#68C19D" ,"#51B69C", "#3BAA9A", "#259E98", "#109195" ,"#028491"
,"#09778D" ,"#186A87" ,"#245D81" ,"#2E4F79")Import Africa shapefile
download.file('http://www.maplibrary.org/library/stacks/Africa/Africa_SHP.zip',destfile = 'Africa_SHP.zip')
unzip('Africa_SHP.zip')
# Define projection system
proj <- CRS('+proj=longlat +ellps=WGS84')## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum Unknown based on WGS84 ellipsoid in CRS definition# import shapefile
Africa_shp<- readOGR('Africa.shp')## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\...\Africa.shp", layer: "Africa"
## with 762 features
## It has 3 fieldsplot(Africa_shp)
Display monthly average of precipitation
## Precipitation
levelplot(mon_mean_rainfall,col.regions=panel_prcp,
main=list('Monthly Average of Precipitation : 1901 - 2019'),
colorkey = list(title = "[mm]",
title.gpar = list(cex = 0.8,
font = 3,
col = '#2E4F79')))+
layer(sp.lines(Africa_shp, lwd=0.8, col='darkgray'))
Display monthly average of temperature
levelplot(mon_mean_temp,col.regions =panel_temp,
main=list('Monthly Average of Temperature : 1901 - 2019'),
colorkey = list(title = "[Deg C]",
title.gpar = list(cex = 1,
font = 3,
col = '#7D0025')))+
layer(sp.lines(Africa_shp, lwd=0.8, col='darkgray'))
Hovmöller diagrams - Precipitation
# Compute weights over latitude dimension
plat<-hovmoller(mon_mean_rainfall, xlab='Latitude',ylab='Months',dirXY=y,
main='Latitudinal variation [mm]',
panel=panel.2dsmoother, args=list(span=0.3)
,col.regions = panel_prcp)
# Compute weights over Longitude dimension
plon<-hovmoller(mon_mean_rainfall, xlab='Longitude',ylab='Months',dirXY=x,
main='Longitudinal variation [mm]',
col.regions = panel_prcp,
panel=panel.2dsmoother, args=list(span=0.3))
grid.arrange(plat, plon, ncol=2)
Hovmöller diagrams - Temperature
# Compute weights over latitude dimension
tlat<-hovmoller(mon_mean_temp, xlab='Latitude',ylab='Months',
main='Latitudinal variation [°C]',
panel=panel.2dsmoother, args=list(span=0.3)
,col.regions = panel_temp)
# Compute weights over Longitude dimension
tlon<-hovmoller(mon_mean_temp, xlab='Longitude',ylab='Months',dirXY=x,
main='Longitudinal variation[°C]',
col.regions = panel_temp,
panel=panel.2dsmoother, args=list(span=0.3))
grid.arrange(tlat, tlon, ncol=2)
Anomaly Hovmöller diagrams (1976-2019)
The base years range considered is from 1976-2005.
Function anomaly
## Extract month value from object
month <- function(x)as.integer(format(getZ(x),'%m'))
## Compute anomaly across the dataset with the function anomaly
anomaly <- function(input,start,midend,end){
# Define time series for the base year and the new set of files
new_set <- seq(as.Date(start), as.Date(end), 'month')
baseyear <- seq(as.Date(start), as.Date(midend), 'month')
new_subset=subset(input,which(dates==start):which(dates==end))
new_subset <- setZ(new_subset,new_set)
names(new_subset)<- as.character(new_set)
base <-subset(input,which(dates==start):which(dates==midend))
base <- setZ(base,baseyear)
## Monthly means across the base year
mm <-stackApply(base, month(base),mean)
## Monthly standard deviation across the base year
msd <-stackApply(new_subset, month(new_subset),sd)
## anomaly
tt=(new_subset - mm)/msd
tt=setZ(tt, new_set)
return(tt)
}Perform anomaly on the datasets
start='1976-01-16'
midend='2005-12-16'
end='2019-12-16'
WA_temp_anom=anomaly(WA_temp,start,midend,end)
WA_rain_anom=anomaly(WA_rainfall,start,midend,end)Temperature anomaly
anomtemp_lat=hovmoller(WA_temp_anom, xlab='Latitude',
main='Latidunal anomaly',
col.regions = panel_temp,
panel=panel.2dsmoother, args=list(span=0.3))
anomtemp_lon=hovmoller(WA_temp_anom, xlab='Longitude',dirXY=x,
main='Longitudinal anomaly',
col.regions = panel_temp,
panel=panel.2dsmoother, args=list(span=0.1))
grid.arrange(anomtemp_lat,anomtemp_lon, ncol=2)
Precipitation anomaly
anomrain_lat<-hovmoller(WA_rain_anom, xlab='Latitude',dirXY=y,
main='Latitudinal Anomaly',
panel=panel.2dsmoother, args=list(span=0.3)
,col.regions = panel_prcp)
anomrain_lon<-hovmoller(WA_rain_anom, xlab='Longitude',dirXY=x,
main='Longitudinal Anomaly',
col.regions = panel_prcp,
panel=panel.2dsmoother, args=list(span=0.3))
grid.arrange(anomrain_lat,anomrain_lon, ncol=2)