|
1
|
|
|
2
|
setwd("~/acrobates/adamw/projects/interp/data")
|
|
3
|
library(reshape)
|
|
4
|
library(sp);library(rgdal);
|
|
5
|
library(reshape)
|
|
6
|
library(ncdf4)
|
|
7
|
library(geosphere)
|
|
8
|
library(rgeos)
|
|
9
|
library(multicore);library(RSQLite)
|
|
10
|
|
|
11
|
stInfo=function(file,inventory,vars=vars){
|
|
12
|
### read in station location data
|
|
13
|
cs=c(11,-1,8,-1,9,-1,6,-4,30,-1,3,-1,3,-1,5) #column widths, negative means skip
|
|
14
|
print("Importing station data")
|
|
15
|
st=read.fwf(file,widths=cs,colClasses="character",comment.char = "",header=F,skip=0)
|
|
16
|
colnames(st)=c("id","lat","lon","elev","name","gsnflag","hcnflag","wmoid")
|
|
17
|
print("Updating factors")
|
|
18
|
for(i in c(5,6,7,8)) st[,i]=as.factor(st[,i])
|
|
19
|
for(i in 2:4) st[,i]=as.numeric(st[,i])
|
|
20
|
### inventory
|
|
21
|
print("Getting station start-stop years from station inventory file")
|
|
22
|
cs=c(11,-1,8,-1,9,-1,4,-1,4,-1,4) #column widths, negative means skip
|
|
23
|
stinv=read.fwf(inventory,widths=cs,colClasses="character",comment.char = "",header=F,skip=0)
|
|
24
|
colnames(stinv)=c("id","lat","lon","variable","yearstart","yearstop")
|
|
25
|
stinv=stinv[stinv$variable%in%vars,]
|
|
26
|
### reshape yearstart
|
|
27
|
stinvw1=cast(stinv,id+lat+lon~variable,value="yearstart")
|
|
28
|
colnames(stinvw1)[!grepl("id|lat|lon",colnames(stinvw1))]=
|
|
29
|
paste(colnames(stinvw1)[!grepl("id|lat|lon",colnames(stinvw1))],".start",sep="")
|
|
30
|
### reshape yearstop
|
|
31
|
stinvw2=cast(stinv,id+lat+lon~variable,value="yearstop")
|
|
32
|
colnames(stinvw2)[!grepl("id|lat|lon",colnames(stinvw2))]=
|
|
33
|
paste(colnames(stinvw2)[!grepl("id|lat|lon",colnames(stinvw2))],".stop",sep="")
|
|
34
|
stinvw=merge(stinvw1,stinvw2)
|
|
35
|
for(i in colnames(stinvw[,!grepl("id|lat|lon",colnames(stinvw))])) stinvw[,i]=as.numeric(as.character( stinvw[,i]))
|
|
36
|
colnames(stinvw)=sub("TMAX","tmax",colnames(stinvw))
|
|
37
|
colnames(stinvw)=sub("TMIN","tmin",colnames(stinvw))
|
|
38
|
colnames(stinvw)=sub("PRCP","ppt",colnames(stinvw))
|
|
39
|
### merge the data
|
|
40
|
print("Merging start-stop dates with station locations")
|
|
41
|
st2=merge(st,stinvw[,!grepl("lat|lon",colnames(stinvw))],all.x=T,by=c("id"))
|
|
42
|
## add station type field
|
|
43
|
st2$type=ifelse((!is.na(st2$tmax.start)|!is.na(st2$tmin.start))&!is.na(st2$ppt.start),"T&P",
|
|
44
|
ifelse((!is.na(st2$tmax.start)|!is.na(st2$tmin.start))&is.na(st2$ppt.start),"T",
|
|
45
|
ifelse((is.na(st2$tmax.start)|is.na(st2$tmin.start))&!is.na(st2$ppt.start),"P","None")))
|
|
46
|
### Convert to spatial points and return
|
|
47
|
print("Convert to spatial points")
|
|
48
|
coordinates(st2)=c("lon","lat")
|
|
49
|
st2@data[,c("lon","lat")]=coordinates(st2)
|
|
50
|
proj4string(st2)=CRS("+proj=longlat +datum=WGS84")
|
|
51
|
return(st2)
|
|
52
|
}
|
|
53
|
|
|
54
|
|
|
55
|
### Process the station data
|
|
56
|
st=stInfo(file="ghcn/ghcnd-stations.txt",inventory="ghcn/ghcnd-inventory.txt",
|
|
57
|
vars=c("PRCP","TMAX","TMIN"))
|
|
58
|
|
|
59
|
writeOGR(st,".","stationlocations",driver="ESRI Shapefile",overwrite=T)
|
