2.1 Packages

#install the RNCEP, lubridate and tidyverse packages
if(!require("RNCEP")) install.packages("RNCEP")
if(!require("lubridate")) install.packages("lubridate")
if(!require("tidyverse")) install.packages("tidyverse")
if(!require("sf")) install.packages("sf")

#load the packages
library(RNCEP)
library(lubridate) #date and time manipulation
library(tidyverse) #data manipulation and visualization
library(RColorBrewer) #color schemes
library(sf) #to import a spatial object and to work with geom_sf in ggplot2

2.2 Data download

We will download the air temperature of the 850haPa pressure level for the year 2016. The variables and pressure levels can be found in the details of the function ?NCEP.gather. The reanalysis2 argument allows us to download both version I and version II, being by default FALSE, that is, we access reanalysis I. In all the requests we will obtain data of every 6 hours (00:00, 06:00, 12:00 and 18:00). This supposes a total of 1464 values for the year 2016.

#define the necessary arguments
month_range <- c(1,12)     #period of months
year_range <- c(2016,2016) #period of years

lat_range <- c(30,60)      #latitude range
lon_range <- c(-30,50)     #longitude range
 

data <- NCEP.gather("air",    #name of the variable
                    850, #pressure level 850hPa
                    month_range,year_range,
                    lat_range,lon_range,
                    return.units = TRUE,
                    reanalysis2=TRUE)

## [1] Units of variable 'air' are degK
## [1] Units of variable 'air' are degK

#dimensions                    
dim(data) 

## [1]   13   33 1464

#we find lon, lat and time with dimnames()
#date and time
date_time <- dimnames(data)[[3]]
date_time <- ymd_h(date_time)
head(date_time)

## [1] "2016-01-01 00:00:00 UTC" "2016-01-01 06:00:00 UTC"
## [3] "2016-01-01 12:00:00 UTC" "2016-01-01 18:00:00 UTC"
## [5] "2016-01-02 00:00:00 UTC" "2016-01-02 06:00:00 UTC"

#longitude and latitude
lat <- dimnames(data)[[1]]
lon <- dimnames(data)[[2]]
head(lon);head(lat)

## [1] "-30"   "-27.5" "-25"   "-22.5" "-20"   "-17.5"

## [1] "60"   "57.5" "55"   "52.5" "50"   "47.5"

2.3 Monthly average

We see that the downloaded data is an array of three dimensions with [lat, lon, time]. As above mentioned, we extracted latitude, longitude and time. The temperature is given in Kelvin. The objective in the next section will be to show two maps comparing January and July.

#create our grouping variable
group <- month(date_time) 

#estimate the average temperature by month 
data_month <- aperm(
  apply(
    data, #our data
    c(1,2), #apply to each time series 1:row, 2:column a the mean( ) function
    by, #group by
    group, #months
    function(x)ifelse(all(is.na(x)),NA,mean(x))),
  c(2,3,1)) #reorder to get an array like the original

dim(data_month) #850haPa temperature per month January to December

## [1] 13 33 12

2.4 Visualization

Once we got here, we can visualize the 850hPa temperature of January and July with ggplot2. In this example, I use geom_sf( ) from the library sf, which makes the work easier to visualize spatial objects in ggplot (in the near future I will make a post about sf and ggplot). In the dimension of latitude and longitude we saw that it only indicates a value for each row and column. But we need the coordinates of all the cells in the matrix. To create all combinations between two variables we can use the expand.grid( ) function.

#first we create all the combinations of lon-lat
lonlat <- expand.grid(lon=lon,lat=lat)

#as lonlat was a row/column name, it is character, that's why we convert it into numeric
lonlat <- apply(lonlat,2,as.numeric)

#lon and lat are not in the order as we expect
#row=lon; column=lat
data_month <- aperm(data_month,c(2,1,3))

#subtract 273.15K to convert K to ºC.
df <- data.frame(lonlat,
                 Ta01=as.vector(data_month[,,1])-273.15,
                 Ta07=as.vector(data_month[,,7])-273.15)

Before we can make the map with ggplot2, we have to adapt the table. The shapefile with the countries limits can be downloaded here.

#convert the wide table into a long one
df <- gather(df,month,Ta,Ta01:Ta07)%>%
             mutate(month=factor(month,unique(month),c("Jan","Jul")))

#import the countries limits
limit <- st_read("CNTR_RG_03M_2014.shp")

## Reading layer `CNTR_RG_03M_2014' from data source 
##   `E:\GitHub\blog_update_2021\content\en\post\2018-09-15-access-to-climate-reanalysis-data-from-r\CNTR_RG_03M_2014.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 256 features and 3 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -180 ymin: -90 xmax: 180 ymax: 83.66068
## Geodetic CRS:  ETRS89

#color scheme
colbr <- brewer.pal(11,"RdBu")

ggplot(df)+
      geom_tile(aes(lon,lat,fill=Ta))+ #temperature data
      geom_sf(data=limit,fill=NA,size=.5)+ #limits 
        scale_fill_gradientn(colours=rev(colbr))+
          coord_sf(ylim=c(30,60),xlim=c(-30,50))+
          scale_x_continuous(breaks=seq(-30,50,10),expand=c(0,0))+
          scale_y_continuous(breaks=seq(30,60,5),expand=c(0,0))+
          labs(x="",y="",fill="Ta 850hPa (ºC)")+
           facet_grid(month~.)+ #plot panels by month
             theme_bw()

3.1 Installation

if(!require("reticulate")) install.packages("reticulate")
if(!require("ncdf4")) install.packages("ncdf4") #to manage netCDF format

#load packages
library(reticulate)
library(ncdf4)

Once we have installed anaconda and the package reticulate, we can install the library python ecmwfapi. We can carry out the installation, or through the Windows CMD using the command conda install -c conda-forge ecmwf-api-client, or with the R function py_install( ) from the reticulate package. The same function allows us to install any python library from R.

#install the python ECMWF API
py_install("ecmwf-api-client")

3.2 Connection and download with the ECMWF API

In order to access the API, it is required to create a file with the user’s information.

The “.ecmwfapirc” file must contain the following information:

{
    "url"   : "https://api.ecmwf.int/v1",
    "key"   : "XXXXXXXXXXXXXXXXXXXXXX",
    "email" : "john.smith@example.com"
}

The key can be obtained with the user account here.

The file can be created with the Windows notebook.

1. We create a document “ecmwfapirc.txt”.
2. Rename this file to “.ecmwfapirc.”

The last point disappears automatically. Then we save this file in “C:/USERNAME/.ecmwfapirc” or “C:/USERNAME/Documents/.ecmwfapirc”.

#import the python library ecmwfapi
ecmwf <- import('ecmwfapi')

#for this step there must exist the file .ecmwfapirc
server = ecmwf$ECMWFDataServer() #start the connection

One we get here, how do we create a query? The easiest thing is to go to the website of ECMWF, where we choose the database, in this case ERA-Interim surface, to create a script with all the necessary data. More details about the syntax can be found here. When we proceed on the website, we only have to click on “View MARS Request”. This step takes us to the script in python.

With the syntax of the script from the MARS Request, we can create the query in R.

#we create the query
query <-r_to_py(list(
  class='ei',
  dataset= "interim", #dataset
  date= "2017-01-01/to/2017-12-31", #time period
  expver= "1",
  grid= "0.125/0.125", #resolution
  levtype="sfc",
  param= "167.128", # air temperature (2m)
  area="45/-10/30/5", #N/W/S/E
  step= "0",
  stream="oper",
  time="00:00:00/06:00:00/12:00:00/18:00:00", #hours
  type="an",
  format= "netcdf", #format
  target='ta2017.nc' #file name
))

#query to get the ncdf
server$retrieve(query)

The result is a netCDF file that we can process with the library ncdf4.

3.3 Processing ncdf

In the next section, the objective will be the extraction of a time serie from the closest coordinate to a given one. We will use the coordinates of Madrid (40.418889, -3.691944).

#load packages
library(sf)
library(ncdf4)
library(tidyverse)

#open the connection with the ncdf file
nc <- nc_open("ta2017.nc")

#extract lon and lat
lat <- ncvar_get(nc,'latitude')
lon <- ncvar_get(nc,'longitude')
dim(lat);dim(lon)

#extract the time
t <- ncvar_get(nc, "time")

#time unit: hours since 1900-01-01
ncatt_get(nc,'time')

#convert the hours into date + hour
#as_datetime() function of the lubridate package needs seconds
timestamp <- as_datetime(c(t*60*60),origin="1900-01-01")

#import the data
data <- ncvar_get(nc,"t2m")

#close the conection with the ncdf file
nc_close(nc)

In this next section we use the sf package, which is replacing the well known sp and rgdal packages.

#create all the combinations of lon-lat
lonlat <- expand.grid(lon=lon,lat=lat)

#we must convert the coordinates in a spatial object sf
#we also indicate the coordinate system in EPSG code
coord <- st_as_sf(lonlat,coords=c("lon","lat"))%>%
                    st_set_crs(4326)

#we do the same with our coordinate of Madrid
psj <- st_point(c(-3.691944,40.418889))%>%
                   st_sfc()%>%
                     st_set_crs(4326)

#plot all points
plot(st_geometry(coord))
plot(psj,add=TRUE,pch = 3, col = 'red')

In the next steps we calculate the distance of our reference point to all the grid points. Then we look for the one with less distance.

#add the distance to the points
coord <- mutate(coord,dist=st_distance(coord,psj))

#create a distance matrix with the same dimensions as our data
dist_mat <- matrix(coord$dist,dim(data)[-3])

#the arrayInd function is useful to obtain the row and column indexes
mat_index <- as.vector(arrayInd(which.min(dist_mat), dim(dist_mat)))

#we extract the time serie and change the unit from K to ºC
#we convert the time in date + hour
df <- data.frame(ta=data[mat_index[1],mat_index[2],],time=timestamp)%>%
        mutate(ta=ta-273.15,time=ymd_hms(time))

Finally, we visualize our time series.

ggplot(df,
       aes(time,ta))+
    geom_line()+
    labs(y="Temperature (ºC)",
             x="")+
    theme_bw()