GY667 Lab 2
Task 1: Loading the Data
library(maps)
library(readr)
graphics.off()
rm(list=ls())
malinhead = read.csv(file.path(getwd(),'malin_head.csv'),skip = 1,header=TRUE)Task 2: Mapping the Data
names(malinhead)<-c("lon","lat","time","name","h","flag")
malin.lat <- malinhead$lat[1]
malin.lon <- malinhead$lon[1]
map("world",c("ireland","uk"),fill=TRUE,xlim=c(-12,-4),ylim=c(51,56))
map.axes(cex.axis=1)
title(main="Location of Malin Head Tide Gauge",xlab="Longitude",ylab="Latitude")
points(malin.lon,malin.lat,pch=21,col="gray",bg="red")
text(malin.lon-.5,malin.lat,"Malin.",col="red")
Task 3: POSIXlt Format
head(malinhead$time)## [1] 2008-11-13T15:06:00Z 2008-11-13T15:12:00Z 2008-11-13T15:18:00Z
## [4] 2008-11-13T15:24:00Z 2008-11-13T15:30:00Z 2008-11-13T15:36:00Z
## 745330 Levels: 2008-11-13T15:06:00Z ... 2019-02-20T15:50:00Zmalin.time<-as.POSIXlt(malinhead$time, format="%Y-%m-%dT%H:%M:%SZ",tz='UTC')
head(malin.time)## [1] "2008-11-13 15:06:00 UTC" "2008-11-13 15:12:00 UTC"
## [3] "2008-11-13 15:18:00 UTC" "2008-11-13 15:24:00 UTC"
## [5] "2008-11-13 15:30:00 UTC" "2008-11-13 15:36:00 UTC"malin <- data.frame("time"=malin.time,
"year"=malin.time$year+1900,
"month"=malin.time$mon+1,
"day"=malin.time$mday,
"hour"=malin.time$hour+1,
"min"=malin.time$min,
"sec"=malin.time$sec,
"h"=malinhead$h,
"flag"=malinhead$flag)Malin.hour
New time series, starting on the hour
head(malin,30)## time year month day hour min sec h flag
## 1 2008-11-13 15:06:00 2008 11 13 16 6 0 0.348 1
## 2 2008-11-13 15:12:00 2008 11 13 16 12 0 0.434 1
## 3 2008-11-13 15:18:00 2008 11 13 16 18 0 0.532 1
## 4 2008-11-13 15:24:00 2008 11 13 16 24 0 0.638 1
## 5 2008-11-13 15:30:00 2008 11 13 16 30 0 0.701 1
## 6 2008-11-13 15:36:00 2008 11 13 16 36 0 0.758 1
## 7 2008-11-13 15:42:00 2008 11 13 16 42 0 0.847 1
## 8 2008-11-13 15:48:00 2008 11 13 16 48 0 0.927 1
## 9 2008-11-13 15:54:00 2008 11 13 16 54 0 0.955 1
## 10 2008-11-13 16:00:00 2008 11 13 17 0 0 1.066 1
## 11 2008-11-13 16:06:00 2008 11 13 17 6 0 1.131 1
## 12 2008-11-13 16:12:00 2008 11 13 17 12 0 1.223 1
## 13 2008-11-13 16:18:00 2008 11 13 17 18 0 1.239 1
## 14 2008-11-13 16:24:00 2008 11 13 17 24 0 1.334 1
## 15 2008-11-13 16:30:00 2008 11 13 17 30 0 1.356 1
## 16 2008-11-13 16:36:00 2008 11 13 17 36 0 1.448 1
## 17 2008-11-13 16:42:00 2008 11 13 17 42 0 1.496 1
## 18 2008-11-13 16:48:00 2008 11 13 17 48 0 1.547 1
## 19 2008-11-13 16:54:00 2008 11 13 17 54 0 1.563 1
## 20 2008-11-13 17:00:00 2008 11 13 18 0 0 1.585 1
## 21 2008-11-13 17:06:00 2008 11 13 18 6 0 1.645 1
## 22 2008-11-13 17:12:00 2008 11 13 18 12 0 1.643 1
## 23 2008-11-13 17:18:00 2008 11 13 18 18 0 1.706 1
## 24 2008-11-13 17:24:00 2008 11 13 18 24 0 1.736 1
## 25 2008-11-13 17:30:00 2008 11 13 18 30 0 1.743 1
## 26 2008-11-13 17:36:00 2008 11 13 18 36 0 1.749 1
## 27 2008-11-13 17:42:00 2008 11 13 18 42 0 1.782 1
## 28 2008-11-13 17:48:00 2008 11 13 18 48 0 1.790 1
## 29 2008-11-13 17:54:00 2008 11 13 18 54 0 1.786 1
## 30 2008-11-13 18:00:00 2008 11 13 19 0 0 1.755 1tail(malin,n=11)## time year month day hour min sec h flag
## 745320 2019-02-20 15:00:00 2019 2 20 16 0 0 -0.709 0
## 745321 2019-02-20 15:05:00 2019 2 20 16 5 0 -0.662 0
## 745322 2019-02-20 15:10:00 2019 2 20 16 10 0 -0.546 0
## 745323 2019-02-20 15:15:00 2019 2 20 16 15 0 -0.538 0
## 745324 2019-02-20 15:20:00 2019 2 20 16 20 0 -0.402 0
## 745325 2019-02-20 15:25:00 2019 2 20 16 25 0 -0.332 0
## 745326 2019-02-20 15:30:00 2019 2 20 16 30 0 -0.221 0
## 745327 2019-02-20 15:35:00 2019 2 20 16 35 0 -0.164 0
## 745328 2019-02-20 15:40:00 2019 2 20 16 40 0 -0.069 0
## 745329 2019-02-20 15:45:00 2019 2 20 16 45 0 -0.003 0
## 745330 2019-02-20 15:50:00 2019 2 20 16 50 0 0.069 0malin.sub <- malin[10:745320,]
malin.hour <- aggregate(h~hour+day+month+year,malin.sub,mean)
malin.hour <- cbind(malin.hour, NA)
names(malin.hour) <- c("hour","day","month", "year", "h", "time")
malin.hour$time <- as.POSIXlt(sprintf("%s/%s/%s %s",
malin.hour$year, malin.hour$month,
malin.hour$day, malin.hour$hour),
format="%Y/%m/%d %H",tz='UTC')
summary(malin.hour)## hour day month year
## Min. : 1.00 Min. : 1.00 Min. : 1.000 Min. :2008
## 1st Qu.: 7.00 1st Qu.: 8.00 1st Qu.: 3.000 1st Qu.:2011
## Median :13.00 Median :16.00 Median : 6.000 Median :2013
## Mean :12.51 Mean :15.78 Mean : 6.469 Mean :2013
## 3rd Qu.:19.00 3rd Qu.:23.00 3rd Qu.:10.000 3rd Qu.:2016
## Max. :24.00 Max. :31.00 Max. :12.000 Max. :2019
## h time
## Min. :-2.25990 Min. :2008-11-13 17:00:00
## 1st Qu.:-0.67800 1st Qu.:2011-01-01 09:45:00
## Median : 0.04070 Median :2013-04-11 07:30:00
## Mean : 0.01898 Mean :2013-07-08 01:10:54
## 3rd Qu.: 0.71790 3rd Qu.:2016-05-20 00:15:00
## Max. : 2.70610 Max. :2019-02-20 16:00:00head(malin.hour)## hour day month year h time
## 1 17 13 11 2008 1.3403 2008-11-13 17:00:00
## 2 18 13 11 2008 1.7165 2008-11-13 18:00:00
## 3 19 13 11 2008 1.6967 2008-11-13 19:00:00
## 4 20 13 11 2008 1.2760 2008-11-13 20:00:00
## 5 21 13 11 2008 0.5830 2008-11-13 21:00:00
## 6 22 13 11 2008 -0.2266 2008-11-13 22:00:00tail(malin.hour)## hour day month year h time
## 74175 11 20 2 2019 -0.435250 2019-02-20 11:00:00
## 74176 12 20 2 2019 -1.201000 2019-02-20 12:00:00
## 74177 13 20 2 2019 -1.613750 2019-02-20 13:00:00
## 74178 14 20 2 2019 -1.599182 2019-02-20 14:00:00
## 74179 15 20 2 2019 -1.169417 2019-02-20 15:00:00
## 74180 16 20 2 2019 -0.709000 2019-02-20 16:00:00Task 4: Malin.mon
New time series, starting with the first full month and ending with the last.
min(which(malin.hour$year==2008&malin.hour$month==12))## [1] 417max(which(malin.hour$year==2019&malin.hour$month==1))## [1] 73708malin.sub2 <- malin.hour[417:73708,]
malin.mon <- aggregate(h~month+year,malin.sub2,mean)
malin.mon <- cbind(malin.mon, NA)
names(malin.mon) <- c("month", "year", "h", "time")
malin.mon$time <- as.POSIXlt(sprintf("%s/%s/15", malin.mon$year, malin.mon$month),tz="UTC")
malin.mon <- cbind(malin.mon, NA, NA, NA, NA)
names(malin.mon) <- c("month", "year", "h", "time",
"mhw","mlw","mtl","n")Task 5: Plotting the Monthly Timeseries
plot(malin.mon$time, malin.mon$h, type="o", col="red", main="Monthly Tide Gauge", xlab="Years", ylab="Tide")
legend("topleft",
legend=c("Malin Head"),
fill=c("red"), lty=1, lwd=1, bty="n")
Task 5 Data Automation (This Code Works!)
install.packages(“devtools”) devtools::install_github(“troyhill/VulnToolkit”, force = TRUE) library(VulnToolkit)
for(j in 1:length(malin.mon\(month)){ ind <- which(malin.hour\)year == malin.mon\(year[j] & malin.hour\)month==malin.mon$month[j])
t <- malin.hour\(time[ind] h <- malin.hour\)h[ind]
hl<-HL(h,t) HL.plot(h,t)
cat(“you are on index”, j, " of “, length(malin.mon\(month), "\n") cat("in ", malin.mon\)year[j],”month = “, malin.mon\(month[j],"\n") cat("you have ", length(hl\)level[hl$tide == “H”]),” high tides“) cat(”you have “, length(hl\(level[hl\)tide ==”L“]),” low tides“)
HH<-hl\(level[hl\)tide == “H”] LL<-hl\(level[hl\)tide == “L”]
question1 <- readline(“Do the numbers of high tides match low (Y/N)”) if(regexpr(question1, ‘y’, ignore.case = TRUE) == 1){ malin.mon\(mhw[j] <- mean(HH) malin.mon\)mlw[j] <- mean(LL) malin.mon\(mtl[j] <- mean(c(HH,LL)) malin.mon\)n[j] <- length(LL) continue = TRUE next } else if (regexpr(question1, ‘n’, ignore.case = TRUE) == 1){ question2 <- readline(“Would you like to cut out a high or a low tide? (h/l/o)”) if(regexpr(question2, ‘h’, ignore.case = TRUE) == 1){ malin.mon\(mhw[j] <- mean(HH[1:(length(HH)-1)]) malin.mon\)mlw[j] <- mean(LL) malin.mon\(mtl[j] <- mean(c(HH[1:(length(HH)-1)],LL)) malin.mon\)n[j] <- length(LL) cat(“in”, malin.mon\(year[j], "month = ", malin.mon\)month[j],“”) cat(“you now have”, length(HH)-1," high tides“) cat(”you now have “, length(LL),” low tides“) continue = TRUE
} else if (regexpr(question2, ‘l’, ignore.case = TRUE) == 1){ malin.mon\(mhw[j] <- mean(HH) malin.mon\)mlw[j] <- mean(LL[1:(length(LL)-1)]) malin.mon\(mtl[j] <- mean(c(LL[1:(length(LL)-1)],HH)) malin.mon\)n[j] <- length(LL) cat(“in”, malin.mon\(year[j], "month = ", malin.mon\)month[j],“”) cat(“you now have”, length(HH)," high tides“) cat(”you now have “, length(LL)-1,” low tides“) continue = TRUE }else {malin.mon\(mhw[j] <- NA malin.mon\)mlw[j] <- NA malin.mon\(mtl[j] <- NA malin.mon\)n[j] <- length(LL)}
} }
Task 6: Saving the Data
write.csv(malin.mon, file="malin_mon.csv")