| 1 |
be7b1081
|
Adam M. Wilson
|
|
| 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)
|
