library("gstat")
library("maps")

# read stationlist
stnlist = read.table("rlr_monthly/rlr_monthly/filelist.txt", sep=";", header=FALSE)

# remove last line (blank)
stnlist = stnlist[1:(nrow(stnlist)-1),]

# get stnids for filenames
stnids = stnlist[,1]

# dim variables to hold trend data
trend = stnids
year1 = stnids
year2 = stnids
datacnt = stnids

# create array of filenames
filenames = paste("rlr_monthly/rlr_monthly/data/", stnids,".rlrdata", sep="")

# ----------------------
# fourth naive cut
# ----------------------
years = seq(1940,2010)
yrtr30 = 31:(length(years)) # dim
yrtr30 = NA # NA

for (y in 31:(length(years)) ) {

  yearbeg = years[y-30]
  yearend = years[y]

  n=1
  for (fn in filenames) {
	arlr = as.matrix(read.table(as.character(fn), sep=";", header=F, na.strings="-99999"))

		trend[n] = NA
		year1[n] = NA
		year2[n] = NA
		datacnt[n] = NA

	# filter on year1, year2
	#arlr = arlr[! is.na(arlr[,2]),]
	arlr = arlr[arlr[,1]>yearbeg & arlr[,1]<yearend & !is.na(arlr[,2]),]

	if (is.integer(nrow(arlr))) {
	  if (nrow(arlr) > ((yearend - yearbeg) * 12 * .80)) {
		# calculate trend and range of years
		trend[n] = lm(arlr[,2] ~ arlr[,1])$coeff[2]
		year1[n] = arlr[1,1]
		year2[n] = arlr[nrow(arlr),1]
		datacnt[n] = nrow(arlr)
	} }
        
	n = n+1
  }
  #hist(trend,breaks=20)
  yrtr30[y-30] = mean(trend,na.rm=TRUE)
  print (paste(yearbeg, "-", yearend, ":", mean(trend,na.rm=TRUE)))
}

trend30 = trend
mean(trend30,na.rm=TRUE)
sum(!is.na(trend30))

# bind trend, start_date, end_date, and data_counts to stationlist
stntrend = cbind(stnlist, trend30,year1,year2,datacnt)

# find largest risers
stntrend[stntrend$trend30>10 & !is.na(stntrend$trend30),]
stntrend[stntrend$trend30< -10 & !is.na(stntrend$trend30),]

# extract stnid, lat, lon, trend
tdata= stntrend[,c(1,2,3,8)]

# remove stations without trends (NA)
tdata = tdata[!is.na(tdata[,4]),]

d = tdata
e = na.omit(d)
names(e)= c("id", "x", "y", "trend")
coordinates(e) <- ~ x+y

x11()
png("psmsl-30-ave.png", bg="white", width=400)
plot(yrtr30 ~ years[(31-15):(length(years)-15)], xlab="",ylab="mm / yr", 
main="PSMSL Average of 30 Year Trends for Data Set",cex.sub="0.8",
sub="not a mean sea level rise / no geographic distribution\ndata: http://www.psmsl.org/data/\nimage: http://rhinohide.wordpress.com")
lines(yrtr30 ~ years[(31-15):(length(years)-15)])
dev.off()

x11()
png("psmsl-30-map.png", bg="white", width=400)
map("world", fill=FALSE)
title(main="PSMSL Stations with Continuous 30 Year Trend", cex.main=1)
rising = coordinates(e[e$trend>0,])
falling = coordinates(e[e$trend<0,])
points(rising[,2], rising[,1], col="red")
points(falling[,2], falling[,1], col="blue")
legend.txt <- c("rising", "falling")
legend.col <- c("red", "blue")
legend(-175,-30, legend.txt, fill=legend.col, cex=0.7)
title(sub="data: http://www.psmsl.org/data/\nimage: http://rhinohide.wordpress.com", cex.sub=0.9)
dev.off()

x11()
png("psmsl-30-hist.png", bg="white",width=400)
hist(trend30, breaks=100, xlab="mm / yr",
main="PSMSL Stations 30 Year Trends", cex.main=1,
sub="data: http://www.psmsl.org/data/\nimage: http://rhinohide.wordpress.com", cex.sub=0.7)
dev.off()

# ----------------------
# fourth naive cut - 17 yr trend
# ----------------------

years = seq(1940,2010)
yrtr17 = 17:(length(years)) # dim
yrtr17 = NA # NA

for (y in 18:(length(years)) ) {

  yearbeg = years[y-17]
  yearend = years[y]

  n=1
  for (fn in filenames) {
	arlr = as.matrix(read.table(as.character(fn), sep=";", header=F, na.strings="-99999"))

		trend[n] = NA
		year1[n] = NA
		year2[n] = NA
		datacnt[n] = NA
	# filter on year1, year2
	#arlr = arlr[! is.na(arlr[,2]),]
	arlr = arlr[arlr[,1]>yearbeg & arlr[,1]<yearend & !is.na(arlr[,2]),]

	if (is.integer(nrow(arlr))) {
	  if (nrow(arlr) > ((yearend - yearbeg) * 12 * .80)) {
		# calculate trend and range of years
		trend[n] = lm(arlr[,2] ~ arlr[,1])$coeff[2]
		year1[n] = arlr[1,1]
		year2[n] = arlr[nrow(arlr),1]
		datacnt[n] = nrow(arlr)
	} }
        
	n = n+1
  }
  #hist(trend,breaks=20)
  yrtr17[y-17] = mean(trend,na.rm=TRUE)
  print (paste(yearbeg, "-", yearend, ":", mean(trend,na.rm=TRUE)))
}

trend17 = trend
mean(trend17,na.rm=TRUE)
sum(!is.na(trend17))
 
# bind trend, start_date, end_date, and data_counts to stationlist
stntrend = cbind(stnlist, trend17,year1,year2,datacnt)

# find largest riser, fastest sinker
stntrend[stntrend$trend17>20 & !is.na(stntrend$trend17),]
stntrend[stntrend$trend17< -20 & !is.na(stntrend$trend17),]

# extract stnid, lat, lon, trend
tdata= stntrend[,c(1,2,3,8)]

# remove stations without trends (NA)
tdata = tdata[!is.na(tdata[,4]),]

d = tdata
e = na.omit(d)
names(e)= c("id", "x", "y", "trend")
coordinates(e) <- ~ x+y

x11()
png("psmsl-17-ave.png", bg="white", width=400)
plot(yrtr17 ~ years[(18-9):(length(years)-9)], xlab="",ylab="mm / yr", 
main="PSMSL Average of 17 Year Trends for Data Set",cex.sub="0.8",
sub="not a mean sea level rise / no geographic distribution\ndata: http://www.psmsl.org/data/\nimage: http://rhinohide.wordpress.com")
lines(yrtr17 ~ years[(18-9):(length(years)-9)])
dev.off()

x11()
png("psmsl-17-map.png", bg="white",width=400)
map("world", fill=FALSE)
title(main="PSMSL Stations with Continuous 17 Year Trend", cex.main=1)
rising = coordinates(e[e$trend>0,])
falling = coordinates(e[e$trend<0,])
points(rising[,2], rising[,1], col="red")
points(falling[,2], falling[,1], col="blue")
legend.txt <- c("rising", "falling")
legend.col <- c("red", "blue")
legend(-175,-30, legend.txt, fill=legend.col, cex=0.7)
title(sub="data: http://www.psmsl.org/data/\nimage: http://rhinohide.wordpress.com", cex.sub=0.9)
dev.off()

x11()
png("psmsl-17-hist.png", bg="white",width=400)
hist(trend17, breaks=100, xlab="mm / yr",
main="PSMSL Stations 17 Year Trends", cex.main=1,
sub="data: http://www.psmsl.org/data/\nimage: http://rhinohide.wordpress.com", cex.sub=0.7)
dev.off()


#==============
# hong kong
#==============

hkids=stnlist[stnlist$V5==611,1]
hklist=stnlist[stnlist$V5==611,]
trend = NULL; trend = hkids; trend = NA
year1 = NULL; year1 = hkids; year1 = NA
year2 = NULL; year2 = hkids; year2 = NA
files = paste("rlr_monthly/rlr_monthly/data/", hkids,".rlrdata", sep="")
n=1
for (f in files) {
 arlr = as.matrix(read.table(as.character(f), sep=";", header=F, na.strings="-99999"))
 trend[n] = lm(arlr[,2] ~ arlr[,1])$coeff[2]
 year1[n] = arlr[1,1]
 year2[n] = arlr[nrow(arlr),1]
 n = n+1
 }
trend
mean(trend, na.rm=T)
hk = cbind(hklist, trend, year1, year2)