Project

General

Profile

Download (2.8 KB) Statistics
| Branch: | Revision:
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)