Data processing and data visualisation of Historic station data
Ciaran Kelly
8/25/2022
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. The following script will produce Time series plot, Station Location Map and Frequency Distribution Plots. This script was produced by Ciaran Kelly & Natascha Seifert. All data used in this paper is attached as .csv files in Station_Investigation.zip folder
Time Series 1
To complie the time series plot tmax data for the month of June 1887 was extracted from a .xsl file named (All_Tmax_data) sourced from knmi explorer and met eireann. The data was then inserted into a new .csv file for ease of use and named (Time_1887.csv”). The folllowing codes produce time series.
# install.packages("gghighlight")
# install.packages("extrafont")
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.1.3library(reshape2)
library(gghighlight)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyverse)## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --## v tibble 3.1.8 v purrr 0.3.4
## v tidyr 1.1.4 v stringr 1.4.0
## v readr 2.1.2 v forcats 0.5.1## Warning: package 'readr' was built under R version 4.1.3## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()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##
## Attaching package: 'raster'## The following object is masked from 'package:dplyr':
##
## selectlibrary(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
#install.packages("weathermetrics")
library(weathermetrics)
library(rnaturalearth)
library(rnaturalearthdata)
#install.packages("rcartocolor")
library(rcartocolor)## Warning: package 'rcartocolor' was built under R version 4.1.3library(metR)##
## Attaching package: 'metR'## The following object is masked from 'package:purrr':
##
## crosslibrary(gridExtra)##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combinelibrary(grid)Important .csv file that contains 1887 tmax data and save as a named data file ie:Time
Time <-read.csv("C:/Eire_Tmax_Investigation/Data_sets/Time_1887.csv")
Time <- Time[,1:8]#change number to suit amount of coloums
names(Time)## [1] "Year" "Armagh" "Markree" "Phoenix.Park" "Sheffield"
## [6] "Kilkenny" "Birr" "Roches.Point"Reshape 1887 tmax data from wide to long
df.melt <- melt(Time, id.vars=c("Year"))
head(df.melt)## Year variable value
## 1 1 Armagh 16.9
## 2 2 Armagh 14.4
## 3 3 Armagh 15.3
## 4 4 Armagh 19.2
## 5 5 Armagh 20.0
## 6 6 Armagh 17.3df.melt <- na.omit(df.melt)
df.melt = rename(df.melt, Stations = variable)Plot Time series for 1887
df.melta <- ggplot(df.melt, aes(x=Year , y=value, col=Stations)) +
geom_line(size = 1)+
xlab('June') +
ylab('Temperature (°C)')+
scale_y_continuous(limits = c(10,36), breaks=seq(10,36,2)) +
scale_x_continuous(limits = c(1,31), breaks=seq(1,31,1))+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+
ggtitle("Max Temperature 1887") +
theme(plot.title = element_text(hjust = 0.5))+
theme(plot.title = element_text(hjust = 0.5))+
theme_bw()+ geom_point(aes(x=26, y=33.3),
shape = 21, colour = "black", fill = "red", size = 2.5, stroke =1.5)
df.melta
Overlay point data frame A ontop of data frame gg1 and save in plots folder with a res=300
png("Time_Series_1887.png", width = 12, height = 9, units = 'in', res = 300)
plot(df.melta)Time Series 2
Important .csv file that contains 1887 tmax data and save as a named data file ie:Time
Time <-read.csv("C:/Eire_Tmax_Investigation/Data_sets/Dublin_1876.csv")
Time <- Time[,1:4]#change number to suit amount of coloums
names(Time)## [1] "Year" "Fitzwilliam_Square_Dublin"
## [3] "Phoenix_Park" "Botanic_Gardens"Reshape 1887 tmax data from wide to long
df.melt <- melt(Time, id.vars=c("Year"))
head(df.melt)## Year variable value
## 1 1 Fitzwilliam_Square_Dublin 22.6
## 2 2 Fitzwilliam_Square_Dublin 21.8
## 3 3 Fitzwilliam_Square_Dublin 23.2
## 4 4 Fitzwilliam_Square_Dublin 17.4
## 5 5 Fitzwilliam_Square_Dublin 19.8
## 6 6 Fitzwilliam_Square_Dublin 20.1df.melt <- na.omit(df.melt)
df.melt = rename(df.melt, Stations = variable)Plot Time series for 1887
df.melta <- ggplot(df.melt, aes(x=Year , y=value, col=Stations)) +
geom_line(size = 1)+
xlab('June') +
ylab('Temperature (°C)')+
scale_y_continuous(limits = c(10,36), breaks=seq(10,36,2)) +
scale_x_continuous(limits = c(1,31), breaks=seq(1,31,1))+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+
ggtitle("Max Temperature Dublin 1876") +
theme(plot.title = element_text(hjust = 0.5))+
theme(plot.title = element_text(hjust = 0.5))+
theme_bw()
df.melta
Overlay point data frame A ontop of data frame gg1 and save in plots folder with a res=300
png("Time_Series_1876.png", width = 12, height = 9, units = 'in', res = 300)
plot(df.melta)Station Location Map 1
The following codes will produce a Station Location Map for all the station used in this reaserch.
library(ggspatial)
library(sf)## Warning: package 'sf' was built under R version 4.1.3## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 7.2.1; sf_use_s2() is TRUElibrary(rnaturalearth)
library(rnaturalearthdata)
library("ggplot2")
library("rgeos")## rgeos version: 0.5-9, (SVN revision 684)
## GEOS runtime version: 3.9.1-CAPI-1.14.2
## Please note that rgeos will be retired by the end of 2023,
## plan transition to sf functions using GEOS at your earliest convenience.
## GEOS using OverlayNG
## Linking to sp version: 1.4-6
## Polygon checking: TRUELoad data frame from the (rnaturalearthdata) library
world <- ne_countries(scale = "medium", returnclass = "sf")
class(world)## [1] "sf" "data.frame"generate world map set the xlim & ylim to your coordinates and save as a named data file ie: gg1
gg1 <-ggplot(data=world)+
geom_sf(color = "black", fill = "lightgreen") +
labs( x = "Longitude", y = "Latitude") +
ggtitle("Station Location 1887") +
coord_sf(xlim = c(-11, 2), ylim = c(50,60), expand = FALSE) +
annotation_scale(location = "bl", width_hint = 0.3) +
annotation_north_arrow(location = "tl", which_north = "true",
height = unit(2.5, "cm"),
width = unit(2.5, "cm"),
pad_x = unit(0.01, "in"), pad_y = unit(0.22, "in"),
style = north_arrow_fancy_orienteering) +
theme(text=element_text(family="Arial",face="bold", size=9))+
theme_bw()Create point data to be overlayed on the map as a named data frame called A
A <- data.frame(
long = c(-6.65, -8.47, -6.3192, -7.88, -8.2208,-1.49, -7.254),
lat = c(54.35, 54.13, 53.3639, 53.08,51.7953,53.381, 52.661 ),
Stations = c("Armagh","Markree","Phoenix Park","Birr", "Roches Point","Sheffield", "Kilkenny"),
stringsAsFactors = FALSE
) Overlay point data frame A ontop of data frame gg1
gg1 +
geom_point(data = A, aes(x = long, y = lat, color=Stations) , size = 2.5)## Scale on map varies by more than 10%, scale bar may be inaccurate
# Station Location Map 2 #
Load data frame from the (rnaturalearthdata) library
world <- ne_countries(scale = "medium", returnclass = "sf")
class(world)## [1] "sf" "data.frame"generate world map set the xlim & ylim to your coordinates and save as a named data file ie: gg1
gg1 <-ggplot(data=world)+
geom_sf(color = "black", fill = "lightgreen") +
labs( x = "Longitude", y = "Latitude") +
ggtitle("Station Location 1876") +
coord_sf(xlim = c(-11, -5), ylim = c(51,56), expand = FALSE) +
annotation_scale(location = "bl", width_hint = 0.3) +
annotation_north_arrow(location = "tl", which_north = "true",
height = unit(2.5, "cm"),
width = unit(2.5, "cm"),
pad_x = unit(0.01, "in"), pad_y = unit(0.22, "in"),
style = north_arrow_fancy_orienteering) +
theme(text=element_text(family="Arial", size=9))+
theme_bw()Create point data to be overlayed on the map as a named data frame called A
A <- data.frame(
long = c(-7.73, -7.36, -7.2, -7.25, -7.13, -7.82, -7.88, -7.93),
lat = c(53.22, 53.47, 53.19, 53.39, 53.11 , 52.65, 53.08, 53.66),
Stations = c("Boora","Mullingar","Clonsast Bord Na Mona","Derrygreenagh", "Ballybrittas", "Thurles", "Birr", "Lanesboro"),
stringsAsFactors = FALSE
) Overlay point data frame A ontop of data frame gg1
gg1 +
geom_point(data = A, aes(x = long, y = lat, color=Stations) , size = 2.5)## Scale on map varies by more than 10%, scale bar may be inaccurate
#** Station Investigation Histogram **# #*
setwd("C:/Eire_Tmax_Investigation/Data_sets/")To following codes complie the Frequency Distribution Plots
#preamble: load package, clear variables, setwd
library(tidyverse)
library(gridExtra)
library(grid)
library(ggplot2)
library(lattice)# Makree
M <-read.csv("C:/Eire_Tmax_Investigation/Data_sets/markree.csv")
Mhist <- ggplot(M, aes(ï..Markree)) +
geom_histogram() +labs(y= "Frequency", x = "Markree (JJA, 2010-2019)") +
geom_point(aes(x=9.9, y=2.5),
shape = 21, colour = "black", fill = "red", size = 2.5, stroke =1.5)+
theme(text=element_text(family="Arial",face="bold", size=9))+
theme_bw()
# Sheffield
S <-read.csv("C:/Eire_Tmax_Investigation/Data_sets/sheffield.csv")
Shist <- ggplot(S, aes(ï..sheffield)) +
geom_histogram() + labs(y= "Frequency", x = "Sheffield (JJA, 2010-2017)") +
geom_point(aes(x=9.4, y=2.5),
shape = 21, colour = "black", fill = "red", size = 2.5, stroke =1.5)+
theme(text=element_text(family="Arial",face="bold", size=9))+
theme_bw()
# PP
P <-read.csv("C:/Eire_Tmax_Investigation/Data_sets/pp.csv")
Phist <- ggplot(P, aes(ï..phoenixpark)) +
geom_histogram() +labs(y= "Frequency", x = "Phoenix Park (JJA, 2013-2019)") +
geom_point(aes(x=9, y=2.5),
shape = 21, colour = "black", fill = "red", size = 2.5, stroke =1.5)+
theme(text=element_text(family="Arial",face="bold", size=9))+
theme_bw()
# Roches Roint
R <-read.csv("C:/Eire_Tmax_Investigation/Data_sets/rpoint.csv")
Rhist <- ggplot(R, aes(ï..rochespoint)) +
geom_histogram() +labs(y= "Frequency", x = "Roches Point (JJA, 2010-2019)") +
geom_point(aes(x=6.6, y=2.5),
shape = 21, colour = "black", fill = "red", size = 2.5, stroke =1.5)+
theme(text=element_text(family="Arial",face="bold", size=9))+
theme_bw()
# Armagh
A <-read.csv("C:/Eire_Tmax_Investigation/Data_sets/armagh.csv")
Ahist <- ggplot(A, aes(ï..Armagh)) +
geom_histogram() +labs(y= "Frequency", x = "Armagh (JJA, 2010-2018)") +
geom_point(aes(x=6.4, y=2.5),
shape = 21, colour = "black", fill = "red", size = 2.5, stroke =1.5)+
theme(text=element_text(family="Arial",face="bold", size=9))+
theme_bw()combine them together
library(ggplot2)
grid.arrange(Mhist, Shist, Phist, Rhist, Ahist, nrow = 2, top = "Tmax (°C) JJA Frequency Distribution 2010-2019")+
theme(text=element_text(family="Arial",face="bold", size=9))+theme_bw()## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## NULL