GY667 Assignment 1
Jordan Brannigan
22 March 2019
Assignment 1: Calculating Sea Level Rise Around Ireland
Introduction
This assignment sets out to accomplish the task of analysing sea level rise around Ireland using data from both the NTGN and PSMSL services. Throughout, a discussion of the R functions executed in order to clean, visualise, analyse, and present data will accompany each part of the document. It is hoped that by doing so, the methods employed here will be sufficiently transparent to allow for reproducible research and results, an important aspect of data analysis and scientific research.
By first investigating tide gauge levels processed by the GY667 class in our two labs, PSMSL data can then be read and compared against overlapping stations and overlapping time periods to assess the difference in the records potentially be used to extend some station records with a degree of error and certainty, similar to the ethos of the ‘virtual’ station approach taken by Jevrejeva et al., (2008; 2014). Lastly, a discussion of the viability of such approaches in the case of irish tide gauges will conclude this document, as well as approaches to creating a mean sea level rise curve for Ireland.
Task 1: Describe the National Tide Gauge Network records
Using the tide gauge data processed by the GY667 class, the monthly mean level measurements could be extracted and read into R in order to map both the locations and measured sea levels.
The working directory was first set to the folder including both the NTGN and PSMSL files.
setwd("C:/Users/Jordan/Desktop/GY667/Assignment 1/GY667_Assignment1")A collection of libraries were loaded in order to provide access to necessary functions for mapping and graphing the gauge data:
library(maps)## Warning: package 'maps' was built under R version 3.5.3Describing the NTGN files - Part 1 (Mapping)
Having inspected the NTGN files on Excel, the majority of them were deemed similar enough to import without column or row modification, bar Howth. As such, a line of code was used to rename the columns in the Howth dataset to match the data from the rest of the gauges. Also, the ‘time’ column contained dates in text strings which were in varying formats e.g. dd-mm-yyyy and dd/mm/yyyy.
Arklow = read.csv(file.path(getwd(),'arklow_monthly.csv'), header=TRUE)
Ballycotton = read.csv(file.path(getwd(),'ballycotton_monthly.csv'), header=TRUE)
Castletownbere = read.csv(file.path(getwd(),'castletownbere_monthly.csv'), header=TRUE)
Dublin = read.csv(file.path(getwd(),'dublin_monthly.csv'), header=TRUE)
Fenit = read.csv(file.path(getwd(),'fenit_monthly.csv'), header=TRUE)
Galway = read.csv(file.path(getwd(),'galway_monthly.csv'), header=TRUE)
Howth = read.csv(file.path(getwd(),'howth_monthly.csv'), header=TRUE)
colnames(Howth) <- c("X", "month", "year", "h", "time")
Malin = read.csv(file.path(getwd(),'malin_monthly.csv'), header=TRUE)
Port_Oriel = read.csv(file.path(getwd(),'port_oriel_monthly.csv'), header=TRUE)
Ringaskiddy = read.csv(file.path(getwd(),'ringaskiddy_monthly.csv'), header=TRUE)
Rossaveel = read.csv(file.path(getwd(),'rossaveel_monthly.csv'), header=TRUE)
Ballyglass = read.csv(file.path(getwd(),'ballyglass_monthly.csv'), header=TRUE)
Dunmore_East = read.csv(file.path(getwd(),'dunmore_east_monthly.csv'), header=TRUE)The majority of datasets had a column named ‘X’ which seemed to be a count of observations or another numeric variable. It was deemed unnecessary for inclusion, and such they were removed:
Arklow$X <- NULL
Ballycotton$X <- NULL
Castletownbere$X <- NULL
Dublin$X <- NULL
Dunmore_East$X <- NULL
Fenit$X <- NULL
Galway$X <- NULL
Howth$X <- NULL
Malin$X <- NULL
Ringaskiddy$X <- NULLWith the ‘X’ column removed, the last thing to do was to format the ‘time’ column in each csv, as various date formats were present and this would be necessary prior to any time series analysis. Interestingly, Ballyglass had no time present but a month and year column with values. The ‘X’ column also proved useful in this case to replace the missing “15” representing the day in for the rest of the gauges. As such, Ballyglass was formatted a little differently below:
library(lubridate)## Warning: package 'lubridate' was built under R version 3.5.3##
## Attaching package: 'lubridate'## The following object is masked from 'package:base':
##
## dateArkTime <-as.POSIXct(Arklow$time,format="%Y-%m-%d",tz='UTC', origin="1970-01-01")
Arklow$time <- ArkTime
BallyCTime <-as.POSIXct(Ballycotton$time,format="%Y-%m-%d",tz='UTC', origin="1970-01-01")
Ballycotton$time <- BallyCTime
Ballyglass$X <- 15
Ballyglass$time <- with(Ballyglass, ymd(sprintf('%04d%02d%02d', year, month, X)))
Ballyglass$X <- NULL
Ballyglass$time <- as.POSIXct(Ballyglass$time, "%Y-%m-%d", tz = 'UTC', origin="1970-01-01")
CastleTime <- as.POSIXct(Castletownbere$time,format="%d/%m/%Y",tz='UTC', origin="1970-01-01")
Castletownbere$time <- CastleTime
DubTime <- as.POSIXct(Dublin$time,format="%Y-%m-%d",tz='UTC', origin="1970-01-01")
Dublin$time <- DubTime
DunTime <- as.POSIXct(Dunmore_East$time,format="%d/%m/%Y",tz='UTC', origin="1970-01-01")
Dunmore_East$time <- DunTime
FenTime <- as.POSIXct(Fenit$time,format="%Y-%m-%d",tz='UTC', origin="1970-01-01")
Fenit$time <- FenTime
GalTime <- as.POSIXct(Galway$time,format="%d/%m/%Y",tz='UTC', origin="1970-01-01")
Galway$time <- GalTime
HowTime <- as.POSIXct(Howth$time,format="%Y-%m-%d",tz='UTC', origin="1970-01-01")
Howth$time <- HowTime
MalTime <- as.POSIXct(Malin$time,format="%d/%m/%Y",tz='UTC', origin="1970-01-01")
Malin$time <- MalTime
PortTime <- as.POSIXct(Port_Oriel$time,format="%d/%m/%Y",tz='UTC', origin="1970-01-01")
Port_Oriel$time <- PortTime
RingTime <- as.POSIXct(Ringaskiddy$time,format="%Y-%m-%d",tz='UTC', origin="1970-01-01")
Ringaskiddy$time <- RingTime
RossTime <- as.POSIXct(Rossaveel$time,format="%d/%m/%Y",tz='UTC', origin="1970-01-01")
Rossaveel$time <- RossTimeWith the individual gauge datasets formatted, the plotting of their locations constituted the next task in describing them. In Excel, a ‘lat’ and ‘lon’ column was added to a csv file containing names of the gauges. The latitude and longitude coordinates were sourced from the Marine Institute’s (2019) Tidal Observations repository. Upon loading this csv file into R, the maps() library was used to display the locations of the gauges around Ireland.
All_Gauges = read.csv(file.path(getwd(),'Gauge_labs.csv'), header=TRUE)library(maptools)## Warning: package 'maptools' was built under R version 3.5.3## Loading required package: sp## Warning: package 'sp' was built under R version 3.5.3## Checking rgeos availability: TRUEmap("world",c("ireland","uk"),fill=TRUE,xlim=c(-12,-4),ylim=c(51,56))
map.axes(cex.axis=1)
title(main="Location of Irish Tide Gauges",xlab="Longitude",ylab="Latitude")
points(All_Gauges$lon,All_Gauges$lat, pch=21,col="gray",bg="red")
pointLabel(All_Gauges$lon, All_Gauges$lat, All_Gauges$gauge, method = "SANN", allowSmallOverlap = FALSE, trace = FALSE, doPlot = TRUE, cex = 0.55, offset = 3, col = "darkgrey")
library(leaflet)
library(sp)
library(ggplot2)## Warning: package 'ggplot2' was built under R version 3.5.3library(RColorBrewer)
library(tidyverse)## Warning: package 'tidyverse' was built under R version 3.5.3## -- Attaching packages --------------------------------------- tidyverse 1.2.1 --## v tibble 1.4.2 v purrr 0.2.5
## v tidyr 0.8.1 v dplyr 0.7.8
## v readr 1.1.1 v stringr 1.3.1
## v tibble 1.4.2 v forcats 0.3.0## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x lubridate::as.difftime() masks base::as.difftime()
## x lubridate::date() masks base::date()
## x dplyr::filter() masks stats::filter()
## x lubridate::intersect() masks base::intersect()
## x dplyr::lag() masks stats::lag()
## x purrr::map() masks maps::map()
## x lubridate::setdiff() masks base::setdiff()
## x lubridate::union() masks base::union()library(dplyr)
library(reshape2)##
## Attaching package: 'reshape2'## The following object is masked from 'package:tidyr':
##
## smithslibrary(classInt)## Loading required package: spData## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source'))`An interactive version of the map was created below using the ‘leaflet’ library. This map also had advantages in that it could allow for zooming in to see the exact location of the gauges and also allowed for identification of the gauge name by hovering over the point as well as zooming in to see the local coastline in which the gauge sat.
leaflet(data=All_Gauges,height=430,width=390) %>% addProviderTiles('CartoDB.Positron') %>%
setView(-8,53.5,6) %>% addCircleMarkers(radius = 5, col = "blueviolet", fillColor = col, fillOpacity = 0.2, label = All_Gauges$gauge) ## Assuming "lon" and "lat" are longitude and latitude, respectivelyDescribing the NTGN files - Part 2 (Plotting)
GaugePlot <- ggplot(Arklow,aes(time,h))+geom_line(aes(color="Arklow"))+ geom_line(data=Ballycotton,aes(color="Ballycotton"))+ geom_line(data=Ballyglass,aes(color="Ballyglass"))+ geom_line(data=Castletownbere,aes(color="Castletownbere"))+ geom_line(data=Dublin,aes(color="Dublin"))+ geom_line(data=Dunmore_East,aes(color="Dunmore_East"))+ geom_line(data=Fenit,aes(color="Fenit"))+ geom_line(data=Galway,aes(color="Galway"))+ geom_line(data=Howth,aes(color="Howth"))+ geom_line(data=Malin,aes(color="Malin"))+ geom_line(data=Port_Oriel,aes(color="Port_Oriel"))+ geom_line(data=Ringaskiddy,aes(color="Ringaskiddy"))+ geom_line(data=Rossaveel,aes(color="Rossaveel"))+
ggtitle("National Tide Gauge Network Sea Level Readings")
labs(color="Legend")## $colour
## [1] "Legend"
##
## attr(,"class")
## [1] "labels"GaugePlot+labs(x = "Year", y = "Sea Level (m)")
With the NTGN sea levels formatted and visualised, the PSMSL gauges could then be loaded and analysed.
Downloading and analysing the PSMSL records
At the www.psmsl.org website, the records for Malin Head, Portrush, Belfast, Bangor, and Dublin were downloaded in rlr format.
path <- file.path(getwd(),"psmsl")
dub.psmsl <- read.table(file.path(path,"Dublin.rlrdata.txt"), header = TRUE)
bel.psmsl <- read.table(file.path(path,"Belfast.rlrdata.txt"), header = TRUE)
mal.psmsl <- read.table(file.path(path,"malin_head.rlrdata.txt"), header = TRUE)
bang.psmsl <- read.table(file.path(path,"Bangor.rlrdata.txt"), header = TRUE)
portrush.psmsl <- read.table(file.path(path,"portrush.rlrdata.txt"), header = TRUE)Plotting PSML
Using the plot () function, exploratory data analysis was undertaken.
plot(dub.psmsl$year, dub.psmsl$rlr - 6500, type = 'l', col = 'purple', xlab = 'Year', ylab = ' ', ylim = c(0, 3500), xlim = c(1850, 2020), yaxt = 'n', main = "PSMSL record lengths")
lines(portrush.psmsl$Year, portrush.psmsl$Rlr - 6000, col = 'red')
lines(bel.psmsl$Year, bel.psmsl$Rlr - 5500, col = 'darkgreen')
lines(bang.psmsl$Year, bang.psmsl$Rlr - 5000, col = 'cyan')
lines(mal.psmsl$Year, mal.psmsl$Rlr - 4500, col = 'brown')
legend("topleft",
legend=c("Malin", "Bangor", "Belfast", "Portrush","Dublin"),
col= c("brown", "cyan", "darkgreen", "red","purple"), cex = 1, lty=1,lwd=1.2, bty="n")
Taking a closer look at Dublin
As can be observed in the plot above (where the gauges were arbitrarily offset), Dublin is the longest continuous record. Malin is the second longest. Both of these stations are also in the NTGN dataset, so it was deemed worthwhile looking at the possibility of using the information contained in both to inform any trends, given that Jevrejeva et al. (2008) reference that 50/60 years is the absolute minimum to plausibly consider for assessing sea level trends outside of variability.
In order to potentially combine these datasets from different sources, the contents of both the NTGN and PSMSL datasets for both Dublin and Malin were investigated in Excel. Points of note included:
PSMSL data being relative to specific ‘rlr’ values based on the local datum corrections for the tide gauge of interest,
NTGN data is given in metres, whereas PSMSL is in millimetres, and
There were overlapping periods of observations for Dublin, but not for Malin.
Following the ethos of the virtual station approach and the use of overlapping residuals as well as removing the mean of common overlapping periods (Dangendorf et al., 2017; Hay et al., 2015; Jevrejeva et al., 2008 ), the following approaches were taken to make estimates of trends in an effort to produce robust trends:
Standardising the Dublin Data
First, the ‘-99999’ values were removed in Excel. Using the reference material on the PSMSL website, the rlr correction for the Dublin data was identified as ‘7007’ mms. With this information, a new column was added to the PSMSL file for Dublin, and the mean of psmsl observations were subtracted from the PSMSL tide gauge observations for Dublin:
dub.psmsl <- cbind(NA, dub.psmsl)
colnames(dub.psmsl) <- c("ala", "year", "month", "rlr")
dub.psmsl$ala <- dub.psmsl$rlr
dub.psmsl$ala <- dub.psmsl$ala - 7007
dub.psmsl$ala <- dub.psmsl$ala - mean(dub.psmsl$ala, na.rm = TRUE)Also, the NTGN files from the GY667 lab had sea level measurements in mm format. As such, multiplying by 1000 was necessary to have the data be in the same unit of measurement.
Dublin$h <- Dublin$h*1000With both files standardised, they could then be merged:
Dublin.merge <- merge(Dublin, dub.psmsl, all = TRUE)
Dublin.merge$gauge = NULL
head(Dublin.merge)## month year h time ala rlr
## 1 1 1938 NA <NA> 3.206776 6973
## 2 1 1939 NA <NA> 183.206776 7153
## 3 1 1940 NA <NA> 15.206776 6985
## 4 1 1941 NA <NA> 27.206776 6997
## 5 1 1942 NA <NA> -194.793224 6775
## 6 1 1943 NA <NA> 116.206776 7086plot(Dublin.merge$year, Dublin.merge$h, type = "l", ylim = c(-350, 500), main = "Dublin Tide Gauge (PSMSL and NTGN)")
lines(Dublin.merge$year, Dublin.merge$ala, col = "red")
legend("topleft",
legend=c("PSMSL", "NTGN"),
col= c("Red", "Black"), cex = 1, lty=1,lwd=1.2, bty="n")
Above, it is clear that there is overlapping years in the mid 2000s but there seems to be an offset between the two datasets still.
Finding overlapping data
By inspecting the dataset in RStudio, 2007 to 2009 was identified as the overlapping period and extracted from the merged dataframe.
Dub_overlap <- Dublin.merge[c(315, 396, 476, 556, 636, 716, 796, 876, 957, 71, 153, 235, 316, 397, 477, 557, 637, 717, 797, 877, 958, 72, 154, 236, 317, 398, 478, 558, 638, 718, 798, 878, 959), c(1:6)]A new column named ‘corrected’ was included and the ‘ala’ column from the PSMSL data (i.e. the column from which the rlr reference and mean of the observations was subtracted) was copied into this column. From that, the NTGN observation values for the overlapping period were then subtracted.
Dub_overlap$corrected <- NA
Dub_overlap$corrected <- Dub_overlap$ala
Dub_overlap$corrected <- Dub_overlap$corrected - Dub_overlap$h
head(Dub_overlap)## month year h time ala rlr corrected
## 315 4 2007 -207.68895 2007-04-15 10.20678 6980 217.8957
## 396 5 2007 -94.45065 2007-05-15 115.20678 7085 209.6574
## 476 6 2007 -36.45109 2007-06-15 168.20678 7138 204.6579
## 556 7 2007 -42.19667 2007-07-15 173.20678 7143 215.4034
## 636 8 2007 -103.93555 2007-08-15 111.20678 7081 215.1423
## 716 9 2007 -113.80456 2007-09-15 97.20678 7067 211.0113The mean of the differences between the NTGN and the PSMSL records was added to the NTGN records to ‘stitch’ the observations together, as such, an investigation of the gauge readings for the overlapping period was performed by effectively creating a timeseries of residuals, as was a part of Hay et al.’s (2015) Kalman smoothing approach.
mean(Dub_overlap$corrected)Dublin.merge$h <- Dublin.merge$h + 206.54
head(Dublin.merge)## month year h time ala rlr
## 1 1 1938 NA <NA> 3.206776 6973
## 2 1 1939 NA <NA> 183.206776 7153
## 3 1 1940 NA <NA> 15.206776 6985
## 4 1 1941 NA <NA> 27.206776 6997
## 5 1 1942 NA <NA> -194.793224 6775
## 6 1 1943 NA <NA> 116.206776 7086plot(Dublin.merge$year, Dublin.merge$ala, type = "l", xlab = "Year", ylab = "Sea Level (mm)", main = "Dublin Gauge Sea Level Trends")
lines(Dublin.merge$year, Dublin.merge$h, col = "red")
legend("topleft",
legend=c("PSMSL", "NTGN"),
col= c("Red", "Black"), cex = 1, lty=1,lwd=1.2, bty="n")
Dublin.merge$slrise <- rowMeans(Dublin.merge[c('ala', 'h')], na.rm = TRUE)plot(Dublin.merge$year, Dublin.merge$slrise, type = "p", col = "red", main = "Merged Dublin Tide Gauge Timeseries (PSMSL and NTGN)", xlab = "Sea Level", ylab = "Year")
Above, the ‘Dublin Gauge Sea Level Trends’ plot shows the shift in the NTGN observations based on a continuation of the long term record rather than a baseline of 0 from the mid-2000s.
Taking a closer look at Malin
In assuming the same approach could be taken with Malin, an overlapping period between the PSMSL and NTGN observations would need to be identified also. Using the contextual rlr reference for Malin, as well as the mean for the corrected observations, the same approach was undertaken:
mal.psmsl <- read.table(file.path(path,"malin_head.rlrdata.txt"), header = TRUE)
mal.psmsl$ala <- NA
mal.psmsl$ala <- mal.psmsl$Rlr
mal.psmsl$ala <- mal.psmsl$ala - 7037
mal.psmsl$ala <- mal.psmsl$ala - mean(mal.psmsl$ala, na.rm = TRUE)Malin$h <- Malin$h*1000colnames(mal.psmsl) <- c("year", "month", "Rlr", "ala")
Malin.merge <- merge(Malin, mal.psmsl, all = TRUE)plot(Malin.merge$year, Malin.merge$ala, type = "l", xlab = "Year", ylab = "Sea Level (mm)", main = "Malin Head Gauge Sea Level Trends")
lines(Malin.merge$year, Malin.merge$h, col = "red")
legend("bottomleft",
legend=c("Malin PSMSL", "Malin NTGN"),
col= c("Black", "Red"), cex = 0.75, lty=1,lwd=1.2, bty="n")
Above, a very different result to the Dublin gauge emerges: there is no overlap between the stations. To address this, an approach to bridge this gap was considered by looking at other gauges in the PSMSL datasets which may have data for that missing gap in the mid-2000s and be close enough to warrant an approach similar to the ‘virtual station’ approach employed by Jevrejeva et al (2014) (albeit in a more complex manner for their study). As such, to answer the question, it was possible to use Portrush as a proxy or virtual station to infill Malin’s missing records and this approach with Malin and Dublin constituted an effort to account for the pitfalls of shorter observations by stitching together timeseries to produce a plot longer than the cited 50 years.
Referring back to the plot of PSMSL record lengths, Portrush was identified as having overlapping records and close geographical proximity to Malin Head.
portrush.psmsl <- read.table(file.path(path,"portrush.rlrdata.txt"), header = TRUE)portrush.psmsl$ala <- NA
portrush.psmsl$ala <- portrush.psmsl$Rlr
portrush.psmsl$ala <- portrush.psmsl$ala - 7027
portrush.psmsl$ala <- portrush.psmsl$ala - mean(portrush.psmsl$ala, na.rm = TRUE)colnames(portrush.psmsl) <- c("year", "month", "rlr", "Ala")plot(mal.psmsl$year, mal.psmsl$ala, type = "l", xlab = "Year", ylab = "Sea Level (mm)", main = "Malin Head Gauge Sea Level Trends")
lines(portrush.psmsl$year, portrush.psmsl$Ala, col = "red")
legend("bottomleft",
legend=c("Malin PSMSL", "Portrush PSMSL"),
col= c("Black", "Red"), cex = 1, lty=1,lwd=1.2, bty="n")
PortRMal_merge <- merge(portrush.psmsl, mal.psmsl, all = TRUE)PortRMal_overlap <- PortRMal_merge[c(450:533), c(1:6)]
PortRMal_overlap$minus <- PortRMal_overlap$Ala - PortRMal_overlap$alamean(PortRMal_overlap$minus, na.rm = TRUE)## [1] 31.09884PortRMal_merge$Ala <- PortRMal_merge$Ala + 31.098plot(PortRMal_merge$year, PortRMal_merge$ala, type = "l", xlab = "Year", ylab = "Sea Level (mm)", main = "Malin and Portrush Gauge Sea Level Trends")
lines(PortRMal_merge$year, PortRMal_merge$Ala, col = "red")
lines(Malin$year, Malin$h, col = "green")
legend("bottomleft",
legend=c("Malin PSMSL", "Portrush PSMSL" ,"Malin NTGN"),
col= c("Black", "Red", "Green"), cex = 0.75, lty=1,lwd=1.2, bty="n")
PortRMal_merge$combined <- rowMeans(PortRMal_merge[c('ala', 'Ala')], na.rm = TRUE)plot(PortRMal_merge$year, PortRMal_merge$combined, type = "l", xlab = "Year", ylab = "Sea Level", main = "Malin and Portrush Gauge Timeseries")
lines(Malin$year, Malin$h, col = "red")
legend("bottomleft",
legend=c("Malin & Portrush PSMSL", "Malin NTGN"),
col= c("Black", "Red"), cex = 0.75, lty=1,lwd=1.2, bty="n")
MalinFinalMerge <- merge(PortRMal_merge, Malin, all = TRUE)MalinFinalOverlap <- MalinFinalMerge[c(611:718), c(1:9)]MalinFinalOverlap$difference <- MalinFinalOverlap$h - MalinFinalOverlap$combinedmean(MalinFinalOverlap$difference, na.rm = TRUE)## [1] -20.44027MalinFinalMerge$h <- MalinFinalMerge$h - 20.44MalinFinalMerge$complete <- rowMeans(MalinFinalMerge[c('h', 'combined')], na.rm = TRUE)plot(MalinFinalMerge$year, MalinFinalMerge$complete, type ="l", xlab = "Year", ylab = "Sea Level (mm)", main = "Merged Malin Gauge Record (Virtual Portrush)")
legend("bottomleft",
legend= "Malin Merged Timeseries",
col= "black", cex = 0.75, lty=1,lwd=1.2, bty="n")
Above, the results of taking the mean of the difference between the Portrush and Malin psmsl datasets is apparent. Similarly with applying the approach again to the newer NTGN data for Malin. In considering this approach, subtraction of the mean for the common period took influence from the approach of Dangendorf et al. (2017) who elaborate on the virtual station method and employ a different approach to that of Jevrejeva et al. (2014).
Linear Models
Dublin
Castletownbere$h <- Castletownbere$h*1000
Dunmore_East$h <- Dunmore_East$h*1000
Galway$h <- Galway$h*1000
Ballyglass$h <- Ballyglass$h*1000lm(Castletownbere$year~Castletownbere$h)##
## Call:
## lm(formula = Castletownbere$year ~ Castletownbere$h)
##
## Coefficients:
## (Intercept) Castletownbere$h
## 2.014e+03 1.075e-02lm(Dunmore_East$year~Dunmore_East$h)##
## Call:
## lm(formula = Dunmore_East$year ~ Dunmore_East$h)
##
## Coefficients:
## (Intercept) Dunmore_East$h
## 2.016e+03 3.274e-03lm(Ballyglass$year~Ballyglass$h)##
## Call:
## lm(formula = Ballyglass$year ~ Ballyglass$h)
##
## Coefficients:
## (Intercept) Ballyglass$h
## 2.013e+03 1.711e-03lm(Galway$year~Galway$h)##
## Call:
## lm(formula = Galway$year ~ Galway$h)
##
## Coefficients:
## (Intercept) Galway$h
## 2.012e+03 6.012e-03lm(Dublin.merge$year~Dublin.merge$slrise)##
## Call:
## lm(formula = Dublin.merge$year ~ Dublin.merge$slrise)
##
## Coefficients:
## (Intercept) Dublin.merge$slrise
## 1.976e+03 8.888e-02lm(MalinFinalMerge$year~MalinFinalMerge$complete)##
## Call:
## lm(formula = MalinFinalMerge$year ~ MalinFinalMerge$complete)
##
## Coefficients:
## (Intercept) MalinFinalMerge$complete
## 1.988e+03 4.176e-03SL_rise = read.csv(file.path(getwd(),'Gauge_SLchange.csv'), header=TRUE)Mapping sl change per year
col_fun <- colorNumeric('YlGnBu',SL_rise$lm)leaflet(data=SL_rise,height=450,width=620,) %>% addProviderTiles('CartoDB.Positron') %>%
setView(-8,53.5,6) %>% addCircleMarkers(fillColor=~col_fun(SL_rise$lm),weight=0,fillOpacity = 0.85, popup=SL_rise$gauge) %>%
addLegend(pal=col_fun,values=SL_rise$lm,title="Sea Level Rise (mm/yr)",position='bottomright')## Assuming "lon" and "lat" are longitude and latitude, respectivelyAbove, an interesting but trend emerges. Dublin appears to have experienced a rate of sea level rise magnitudes above the rest of the gauges. However, it is worth noting that Dublin’s timeseries is much longer than the rest. To have a look at this further, a ‘plotJenks’ function was created using code provided by Alberti (2018):
plotJenks <- function(data, n=3, brks.cex=0.70, top.margin=10, dist=5){
df <- data.frame(sorted.values=sort(data, decreasing=TRUE))
Jclassif <- classIntervals(df$sorted.values, n, style = "jenks") #requires the 'classInt' package
test <- jenks.tests(Jclassif) #requires the 'classInt' package
df$class <- cut(df$sorted.values, unique(Jclassif$brks), labels=FALSE, include.lowest=TRUE) #the function unique() is used to remove non-unique breaks, should the latter be produced. This is done because the cut() function cannot break the values into classes if non-unique breaks are provided
if(length(Jclassif$brks)!=length(unique(Jclassif$brks))){
info <- ("The method has produced non-unique breaks, which have been removed. Please, check '...$classif$brks'")
} else {info <- ("The method did not produce non-unique breaks.")}
loop.res <- numeric(nrow(df))
i <- 1
repeat{
i <- i+1
loop.class <- classIntervals(df$sorted.values, i, style = "jenks")
loop.test <- jenks.tests(loop.class)
loop.res[i] <- loop.test[[2]]
if(loop.res[i]>0.9999){
break
}
}
max.GoF.brks <- which.max(loop.res)
plot(x=df$sorted.values, y=c(1:nrow(df)), type="b", main=paste0("Jenks natural breaks optimization; number of classes: ", n), sub=paste0("Goodness of Fit: ", round(test[[2]],4), ". Max GoF (", round(max(loop.res),4), ") with classes:", max.GoF.brks), ylim =c(0, nrow(df)+top.margin), cex=0.75, cex.main=0.95, cex.sub=0.7, ylab="observation index", xlab="value (increasing order)")
abline(v=Jclassif$brks, lty=3, col="red")
text(x=Jclassif$brks, y= max(nrow(df)) + dist, labels=sort(round(Jclassif$brks, 2)), cex=brks.cex, srt=90)
results <- list("info"=info, "classif" = Jclassif, "breaks.max.GoF"=max.GoF.brks, "class.data" = df)
return(results)
}plotJenks(SL_rise$lm)
## $info
## [1] "The method did not produce non-unique breaks."
##
## $classif
## style: jenks
## one of 10 possible partitions of this variable into 3 classes
## [0.001711,0.006012] (0.006012,0.01075] (0.01075,0.08888]
## 4 1 1
##
## $breaks.max.GoF
## [1] 5
##
## $class.data
## sorted.values class
## 1 0.088880 3
## 2 0.010750 2
## 3 0.006012 1
## 4 0.004176 1
## 5 0.003274 1
## 6 0.001711 1The plotJenks graph above shows that Dublin’s value lies far from the other lm coefficients for each gauge. In fact, the boxplot function was used after and shows this stark difference too:
boxplot(SL_rise$lm)
For purposes of illustration and comparison, a lm function was applied to just the NTGN data for Dublin to show how this affects the relative view of sea level rise per year:
lm(Dublin$year~Dublin$h)##
## Call:
## lm(formula = Dublin$year ~ Dublin$h)
##
## Coefficients:
## (Intercept) Dublin$h
## 2.012e+03 9.718e-03Re-plotting this value for Dublin:
SL_riseDN <- read.csv(file.path(getwd(),'Gauge_SLchange2.csv'), header=TRUE)col_funD <- colorNumeric('YlGnBu',SL_riseDN$lm)leaflet(data=SL_riseDN,height=450,width=620,) %>% addProviderTiles('CartoDB.Positron') %>%
setView(-8,53.5,6) %>% addCircleMarkers(fillColor=~col_funD(SL_riseDN$lm),weight=0,fillOpacity = 0.85, popup=SL_rise$gauge) %>%
addLegend(pal=col_funD,values=SL_riseDN$lm,title="Sea Level Rise (mm/yr)",position='bottomright')## Assuming "lon" and "lat" are longitude and latitude, respectivelyThe difference in the map above shows the effect of removing the longer term lm trend at Dublin as an outlier. However, the first map is important in showing the a noticeably higher level of sea level rise per year at Dublin in recent times.
It is important to consider many contextual factors, especially when the potential outlier is found at the Dublin tide gauge. If one was to consider that methodological or statistical approaches are not responsible for any anomalous values at Dublin, there are many physical For example, Dangendorf et al (2017) reference the role of changes in terrestrial water storage such as ground water depletion and ice melt. Issues such as subsidence of urban land, may also be important. All showing the anthropogenic considerations when accounting for sea level rise.
However, other factors such as seasonal and decadal variability (Hay et al., 2015) as well as Glacial Isostatic Adjustment and geoid changes are paramount in affecting sea level observations on tide gauge records (Dangendorf et al., 2017). Jevrejeva et al. (2008) make explicit reference to the fact that for records shorter than 50 years, decadal variations dominate estimates of acceleration of any sea level rise, and may mask change signals with this noise. Seasonal trends may affect sea level rise also, particularly if temporal coverage is not consistent for one tide gauge or for the sample of the gauge. As Hay et al. (2015) refer to, spatial sparsity and missing data hamper the robustness of approaches such as that used by Jevrejeva et al. It may also be important to consider how tidal cycles and equinoxes may contribute to peaks and troughs in sea level data at tide gauges, as well as seasonal climates such as winter storms in the Irish context. As such, the trends produced above are highly vulnerable to influences from these factors.
Lastly, the impact of glacial isostatic adjustment from the British and Irish Ice Sheet of the last ice age may have produced spatially varying patterns of depression on the crust below. Reflecting on the map above, it is interesting to considering this in a north-east to south-west gradient, particularly along the western half of Ireland, and the increasing trends visible on the map. As such, the ‘fingerprints’ are particularly relevant.
Mean Sea Level Curve for Ireland
Overlaps of all NTGN and average altogether and this mean sea level for 10 years - tricky to merge onto longer records.
The following code and discussion addresses the proposed methods for creating a mean sea level curve for the island of Ireland. Note that this method was only applied to the subset of gauges analysed in this assignment. Ideally, more of the stations would inform this sea level curve.
The first part involved sourcing all of the NTGN sea level observations into one data frame. As detailed above, for Dublin and Malin, the NTGN gauges have been stitched together and comprise much longer datasets now than the rest of the NTGN gauges. As such, a decision to use only NTGN data and specify a time period where overlapping monthly data was available for all of the NTGN gauges was made. The reason being: to minimise the impact of missing values and varying lengths of observations.
As the ‘time’ column had different formats, it was removed from the NTGN datasets.
Malin_N <- Malin
Dublin_N <- Dublin
Dunmore_East_N <- Dunmore_East
Castletownbere_N <- Castletownbere
Galway_N <- Galway
Ballyglass_N <- Ballyglasscolnames(Malin_N) <- c("month", "year", "Malin", "time")
colnames(Dublin_N) <- c("month", "year", "Dublin", "time")
colnames(Dunmore_East_N) <- c("month", "year", "Dunmore East", "time")
colnames(Castletownbere_N) <- c("month", "year", "Castletownbere", "time")
colnames(Galway_N) <- c("month", "year", "Galway", "time")
colnames(Ballyglass_N) <- c("month", "year", "Ballyglass", "time")NTGN <- merge(Malin_N, Dublin_N, na.rm = FALSE, all.x = TRUE, all.y = TRUE)
NTGN <- merge(NTGN, Dunmore_East_N, na.rm = FALSE, all.x = TRUE, all.y = TRUE)
NTGN <- merge(NTGN, Castletownbere_N, na.rm = FALSE, all.x = TRUE, all.y = TRUE)
NTGN <- merge(NTGN, Galway_N, na.rm = FALSE, all.x = TRUE, all.y = TRUE)
NTGN <- merge(NTGN, Ballyglass_N, na.rm = FALSE, all.x = TRUE, all.y = TRUE)NTGN$MSLC <- rowMeans(NTGN[c("Malin", "Dublin", "Dunmore East", "Ballyglass", "Castletownbere", "Galway")], na.rm = TRUE)MSLC <- ggplot(data = NTGN, aes(time, MSLC)) + geom_line(aes(colour = "Irish Mean Sea Level")) +
ggtitle("Mean National Tide Gauge Network Sea Level Readings 2007 - 2019")
MSLC+labs(x = "Year", y = "Sea Level (mm)")
MSL_Trend <- lm(NTGN$year~NTGN$MSLC)
MSL_Trend##
## Call:
## lm(formula = NTGN$year ~ NTGN$MSLC)
##
## Coefficients:
## (Intercept) NTGN$MSLC
## 2.013e+03 8.199e-03A linear model fit to this Irish mean sea level curve gives an annual rise of 0.008199 per year, notably lower than the GSML rise estimates of 3.0 +- 0.7 mm found by Hay et al. (2015).
Global sea level rise
As Hay et al. (2015), Jevrejeva et al. (2008), and Dangendorfet al. (2017) show, regardless of their approach and its associated merits or pitfalls, global sea level rise is a complex process influenced by a myriad of physical processes and also anthropogenic influneces beyond simple greenhouse gas emission. The stuggle to form consensus on 20th century sea level rise often refered to bias and coverage issues with regard to tide gauges. The satellite altimetry era though has been of use in cross validating more recent estimates towards the end of the 20th century. However, going forward, there is no doubt that sea level change will not be uniform. Circulation patterns such as ocean currents, wind, temperature changes, are all expected to change in a warming world and as tipping points and elements of the climate system such as the greenland ice sheet collapse, future acceleration of sea level rise may be uncertain. Local contextual factors will mean that for each tide gauge, the resulting effect may not be uniform.
Bibliography
Alberti(2018) R/plotJenks.R [online]. Available at: https://rdrr.io/github/gianmarcoalberti/GmAMisc/src/R/plotJenks.R (accessed 31 March 2019).
Dangendorf, S., Marcos, M., Wöppelmann, G., Conrad, C. P., Frederikse, T., & Riva, R. (2017). Reassessment of 20th century global mean sea level rise. Proceedings of the National Academy of Sciences, 201616007.
Hay, C. C., Morrow, E., Kopp, R. E., & Mitrovica, J. X. (2015). Probabilistic reanalysis of twentieth-century sea-level rise. Nature, 517(7535), 481.
Jevrejeva, S., Moore, J. C., Grinsted, A., & Woodworth, P. L. (2008). Recent global sea level acceleration started over 200 years ago?. Geophysical Research Letters, 35(8).