Introduction

The 2019-20 class of the MSc Climate Change at Maynooth University were set a group assignment to re-evalaute extreme heat temperature records using Ireland as a proxy. This is the script used to conduct data acquisition, data processing and data visualization of the two reanalysis products used which are 20th Century Reanalysis V3: NOAA and ERA5 cds.climate.copernicus.eu with the parameters prmsl and 850 hpa temperatures ensemble mean gathered at 1500GMT.

Process

NOAA

  • 20th Century Reanalysis V3: NOAA Physical Sciences Laboratory Netcdf files of prmsl and 850 hpa tempertures ensemble mean at 1500GMT were downloaded at its native resolution of 0.5° x 0.5°.
  • Files are imported into R and opened, processed and plotted using various libraries.
  • The variable prmsl is converted into msl.
  • The variable temperature is changed from kelvin to degrees Celsius.
  • The variables are then plotted on a map with temperature shown as colors and msl as contours.

ERA5

  • Search results (copernicus.eu), Netcdf files were downloaded from the Copernicus climate data store onto a local directory.
  • The variables msl and 850 hpa tempertures were downloaded
  • Files are imported into R and opened, processed and plotted.
  • As ERA5 is already in msl no processing was needed for this variable.
  • The variable temperature is changed from kelvin to degrees Celsius.
  • The variables are then plotted on a map with temperature shown as colors and msl as contours.

Install packages and Load libraries

library(ncdf4) # package for netcdf manipulation
library(raster) # package for raster manipulation
## Warning: package 'raster' was built under R version 4.1.3
## Loading required package: sp
## Warning: package 'sp' was built under R version 4.1.3
library(rgdal) # package for geospatial analysis
## Warning: package 'rgdal' was built under R version 4.1.3
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-32, (SVN revision 1176)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.1, released 2021/12/27
## Path to GDAL shared files: C:/Users/ckelly/Documents/R/win-library/4.1/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: C:/Users/ckelly/Documents/R/win-library/4.1/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.5-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
library(ggplot2) # package for plotting
## Warning: package 'ggplot2' was built under R version 4.1.3
#install.packages("weathermetrics")
library(weathermetrics)
library(rnaturalearth)       
library(rnaturalearthdata)
#install.packages("rcartocolor")
library(rcartocolor)
## Warning: package 'rcartocolor' was built under R version 4.1.3
library(metR)
library(gridExtra)
library(grid)

Open local directory and import netcdf files

This script loads the netcdf files one by one. The netcdf files have a naming structure as follows T850_YYYY_mm and MSLP_YYYY. The follow script is an example of loading in the year 1796 and month 07

Open local directory

setwd("C:/Users/Met Eireann/Desktop/2OTH_&_ERA5/May2022/1970S")

Load files

setwd("C:/Users/Met Eireann/Desktop/2OTH_&_ERA5/May2022/1970S")
nc_data <- nc_open('T850_197607.nc')
nc_data_msl <- nc_open('MSLP_197607.nc')

Process 20th netcdf file of T850.

Connect to the NetCDF file using the nc_open() function from the ncdf4 package. We save the connection to the file as a variable called nc.

lon <- ncvar_get(nc_data, "lon")
lon[lon > 180] <- lon[lon > 180] - 360
lat <- ncvar_get(nc_data, "lat", verbose = F)
t <- ncvar_get(nc_data, "time")
head(lon) 
## [1] -22 -21 -20 -19 -18 -17
ndvi.array <- ncvar_get(nc_data, "air") # store the data in a 3-dimensional array
dim(ndvi.array) 
## [1] 30 17
fillvalue <- ncatt_get(nc_data, "air", "_FillValue")
fillvalue
## $hasatt
## [1] FALSE
## 
## $value
## [1] 0
nc_close(nc_data) 

ndvi.array[ndvi.array == fillvalue$value] <- NA


ndvi.slice <- ndvi.array

dim(ndvi.array)
## [1] 30 17
r <- raster(t(ndvi.slice), xmn=min(lon), xmx=max(lon), ymn=min(lat), ymx=max(lat), crs=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))

r <- flip(r, direction='y')
head(r)
##         1      2      3      4      5      6      7      8      9     10     11
## 1  277.01 276.91 276.81 276.91 277.31 277.81 278.41 279.11 279.71 280.21 280.61
## 2  277.21 277.11 277.01 277.21 277.51 278.01 278.61 279.21 279.81 280.31 280.91
## 3  277.21 277.11 277.21 277.41 277.71 278.21 278.81 279.41 280.01 280.61 281.21
## 4  277.31 277.21 277.31 277.51 277.91 278.51 279.11 279.71 280.41 281.01 281.61
## 5  277.51 277.51 277.51 277.81 278.21 278.81 279.41 280.11 280.91 281.61 282.11
## 6  277.81 277.71 277.81 278.01 278.41 279.01 279.71 280.51 281.31 282.11 282.71
## 7  278.01 278.01 278.11 278.31 278.61 279.21 280.01 280.81 281.61 282.41 283.11
## 8  278.31 278.31 278.41 278.61 279.01 279.51 280.21 280.91 281.71 282.41 283.11
## 9  278.71 278.81 278.91 279.11 279.41 279.91 280.51 281.11 281.71 282.41 282.91
## 10 279.31 279.51 279.71 280.01 280.31 280.61 281.01 281.41 281.91 282.31 282.71
##        12     13     14     15     16     17     18     19     20
## 1  281.21 281.81 282.41 283.11 283.61 284.11 284.51 284.81 285.01
## 2  281.41 282.11 282.71 283.41 283.91 284.51 284.91 285.21 285.41
## 3  281.81 282.41 283.01 283.71 284.31 284.81 285.21 285.51 285.71
## 4  282.21 282.71 283.41 284.11 284.71 285.31 285.61 285.81 285.81
## 5  282.71 283.21 283.81 284.51 285.11 285.51 285.81 285.91 285.91
## 6  283.21 283.71 284.31 284.81 285.21 285.41 285.31 285.51 285.71
## 7  283.51 284.01 284.51 285.11 285.51 285.71 285.71 285.71 285.71
## 8  283.51 283.91 284.31 284.91 285.51 285.81 285.91 285.81 285.61
## 9  283.31 283.61 284.01 284.51 285.01 285.41 285.71 285.81 285.71
## 10 283.01 283.21 283.61 284.01 284.51 284.81 285.11 285.31 285.31

Process 20th netcdf file of prmsl.

Connect to the NetCDF file using the nc_open() function from the ncdf4 package. We save the connection to the file as a variable called nc.

setwd("C:/Users/Met Eireann/Desktop/2OTH_&_ERA5/May2022/1970S")

nc_data_msl <- nc_open('MSLP_197607.nc')

names(nc_data_msl$var)
## [1] "prmsl"
lon <- ncvar_get(nc_data_msl, "lon")
lon[lon > 180] <- lon[lon > 180] - 360


lat <- ncvar_get(nc_data_msl, "lat", verbose = F)
t <- ncvar_get(nc_data_msl, "time")

head(lon) 
## [1] -22 -21 -20 -19 -18 -17
ndvi.array_msl <- ncvar_get(nc_data_msl, "prmsl") # store the data in a 3-dimensional array
dim(ndvi.array_msl) 
## [1] 30 17
fillvalue_msl <- ncatt_get(nc_data_msl, "prmsl", "_FillValue")
fillvalue_msl
## $hasatt
## [1] FALSE
## 
## $value
## [1] 0
nc_close(nc_data_msl) 

ndvi.array_msl[ndvi.array_msl == fillvalue_msl$value] <- NA


ndvi.slice_msl <- ndvi.array_msl

dim(ndvi.array_msl)
## [1] 30 17
r_msl <- raster(t(ndvi.slice_msl), xmn=min(lon), xmx=max(lon), ymn=min(lat), ymx=max(lat), crs=CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs+ towgs84=0,0,0"))

r_msl <- flip(r_msl, direction='y')

mslrasdf <-as.data.frame(r_msl, xy=TRUE)
head(mslrasdf)
##           x        y    layer
## 1 -21.51667 61.52941 101967.5
## 2 -20.55000 61.52941 102005.5
## 3 -19.58333 61.52941 102036.5
## 4 -18.61667 61.52941 102061.5
## 5 -17.65000 61.52941 102086.5
## 6 -16.68333 61.52941 102108.5
mslrasdf$layer <-c (mslrasdf$layer /100 )

Set up plots

Next we set up the plot by overlayin both parameters on top of each other.

Set up plots

Run world package

world <- ne_coastline(scale = "large", returnclass = "sf")
class(world)
## [1] "sf"         "data.frame"

Save raster layer as data frames

rasdf <-as.data.frame(r, xy=TRUE)
head(rasdf)
##           x        y  layer
## 1 -21.51667 61.52941 277.01
## 2 -20.55000 61.52941 276.91
## 3 -19.58333 61.52941 276.81
## 4 -18.61667 61.52941 276.91
## 5 -17.65000 61.52941 277.31
## 6 -16.68333 61.52941 277.81

Change kelvin to temp below is the formula to convert F = (9/5)c + 32

rasdf$layer <- kelvin.to.celsius(rasdf$layer)
names(nc_data$var)
## [1] "air"

Create ggplot object

p<- ggplot()+
  geom_raster(aes(x=x,y=y, fill= layer),data=rasdf)+
  geom_sf(fill="transparent", data=world)+
  scale_fill_viridis_c(name="TEMP [°C]",direction=+1,option = "turbo")+
  labs( x = "LONGITUDE", y = "LATITUDE") +
  ggtitle("(02-07-1976)") +
  coord_sf(xlim = c(-20.5, 5.5), ylim = c(47,61), expand = T)+
  theme(text=element_text(family="Arial",face="bold", size=9))+
  theme_bw()

Overlay and Plot

B<- p +
  geom_contour(aes(x=x,y=y,z=layer),color="black",size=0.0005,data=mslrasdf) + 
  geom_text_contour(aes(x=x,y=y,z = layer),color="black", data=mslrasdf)   

B

###set up all plots on one pahe or in lattice plot
#NCEP_plot <- grid.arrange(A,B, nrow = 1, top = "NOAA 20CRv3 MSLP & 850-hPa Temperatures")+
  #theme(text=element_text(family="Arial",face="bold", size=9))

#NCEP_plot