#This file is for generating data sets and runing WinBugs models for hurricane counts and hurricane intensities from the best track data and SST data #generated in chapter 6 by lat lon. #Start in "C:/Users/ThomasJagger/Dropbox/Book" #First we incldue the hourly interpolated best track data set, the sst data set and the covariates #You can run from any chapter (here chapter 6). Change this (chapter 12) when the data sets are collected load("../chap06/best.use.Rdata") load("../chap06/ncdataframe.Rdata") load("../chap07/annual.Rdata") years<-1866:2009 #for setting up data sets #This completes the needed data sets #Define the hexagon tiling #Make copy of data set and add WmaxSM attribute Wind.df=best.use Wind.df$WmaxSM=Wind.df$WmaxS*.5144 require(sp) #Create spatial points dataframe and coordinates(Wind.df)=c("lon","lat") coordinates(ncdataframe)=c("lon","lat") #Add projection information ll = "+proj=longlat +ellps=WGS84" proj4string(Wind.df)=CRS(ll) proj4string(ncdataframe)=CRS(ll) #New projection, though I would use another to reflect the lower lattitudes of hurricanes maybee lat_1=15 lat_2=45 lcc="+proj=lcc +lat_1=60 +lat_2=30 +lon_0=-60" require(rgdal) #Note the transformation does not rename the old coordinates, as.data.frame(Wind.sdf) shows old columns lat lon but with the new coordinates Wind.sdf = spTransform(Wind.df, CRS(lcc)) SST.sdf=spTransform(ncdataframe,CRS(lcc)) # compare make sure that Wind.sdf contained in SST.sdf bbox(Wind.sdf) bbox(SST.sdf) #Generate our spsample hpt=spsample(Wind.sdf,type="hexagonal", n=250,bb=bbox(Wind.sdf)*1.2,offset=c(1,-1)) #Create polygons special function for hexagons hpg=HexPoints2SpatialPolygons(hpt) #Note this is not a dataframe, but it does have slot ids. #Also note that we have missing slots for the SST. Wind.hexid<-overlay(Wind.sdf,hpg) SST.hexid<-overlay(SST.sdf,hpg) #Split on dataframe only (not on hexid) Wind.split<-split(Wind.sdf@data,Wind.hexid) SST.split<-split(SST.sdf@data,SST.hexid) #reassign names to match those of polygons names(SST.split)<-sapply(hpg@polygons, function(x) x@ID)[as.numeric(names(SST.split))] names(Wind.split)<-sapply(hpg@polygons,function(x) x@ID)[as.numeric(names(Wind.split))] #From now we need to compile mean SST after subsetting the SST data to match our track data set. #Some hurricanes are over land only we remove these hexagons Wind.subset<-Wind.split[names(Wind.split)%in%names(SST.split)] SST.subset<-SST.split[names(SST.split)%in%names(Wind.split)] #Now the datasets are in synch. Alternatively we can subset best.use on M==F to remove storms over land. SST.mean<-data.frame(t(sapply(SST.subset,function(x) colMeans(x)))) #For generating the August-October values by year. SSTYearMonth<-strsplit(substring(colnames(SST.mean),2),"M") SSTYear<-as.numeric(sapply(SSTYearMonth,function(x) x[1])) SSTMonth<-as.numeric(sapply(SSTYearMonth,function(x) x[2])) SSTKeep<-which(SSTMonth%in%c(8,9,10)) #get values SSTYear<-SSTYear[SSTKeep] SST.mean.keep<-SST.mean[,SSTKeep] SST.mean.year<-sapply(unique(SSTYear),function(x) as.vector(rowMeans(SST.mean.keep[,which(x==SSTYear)]))) dimnames(SST.mean.year)<-list(id=rownames(SST.mean.keep),Year=paste("Y",unique(SSTYear),sep="")) SST.mean.year.subset<-SST.mean.year[,paste("Y",years,sep="")] SST.pdf<-SpatialPolygonsDataFrame(hpg[rownames(SST.mean.year.subset)],data.frame(SST.mean.year.subset)) #To plot spplot(SST.pdf["Y2005"]) #Now we need to generate a dataset of storms. #First we generate a list of the maximum storms by region, then count the number of storms per year. get.max <- function(x, in.order=FALSE, idfield="Sid", maxfield="Wmax") { use.max <- x[sapply(split( (1:nrow(x)), x[,idfield]), function(i, x) { xi <- x[i]; ii <- i[max(xi)==xi]; ii[1] }, x=x[,maxfield]), ] if(in.order) use.max <- use.max[rev(order(use.max[, maxfield])), ] use.max } Wind.max<-lapply(Wind.subset,function(x) get.max(x,maxfield="WmaxS")) sapply(Wind.max,nrow) Wind.Count<-t(sapply(Wind.max,function(x) table(factor(subset(x,WmaxS>=34 & Yr%in%years)$Yr,level=years)))) colnames(Wind.Count)<-paste("Y",colnames(Wind.Count),sep="") Wind.Count<-data.frame(Wind.Count) Wind.pdf<-SpatialPolygonsDataFrame(hpg[rownames(Wind.Count)],Wind.Count) #To generate data for WinBugs we need to determine the Years of interest and create hexagons require(maptools) #Write the WinBugs adjacency map out require(spdep) sp2WB(Wind.pdf,"hexagons.txt") #gets winbugs map out hexagon.nb<-poly2nb(Wind.pdf) hexagon.splus<-nb2WB(hexagon.nb) #save for latter source("../chap12/writedatafileR.R") writeDatafileR(hexagon.splus,"hexagonadj.txt") #Will write data (so far) for testing WinBugsFields<-subset(alldata,Year%in%years)[,c("soi","sst")] #Removing the local SST mean from the data #GLM Poisson fits to local SST and global SOI hexfit<-lapply(1:nrow(Wind.Count),function(i) glm(unlist(Wind.Count[i,])~unlist(SST.mean.year.subset[i,])+WinBugsFields$soi,family="poisson")) hexfitz<-t(sapply(hexfit,function(x) summary(x)$coef[,3])) #get Z value rownames(hexfitz)<-rownames(Wind.Count) colnames(hexfitz)<-c("Intercept","Local SST","SOI") Wind.t<-SpatialPolygonsDataFrame(hpg[rownames(hexfitz)],data.frame(hexfitz)) spplot(Wind.t) #For example #This plot shows 108 independent hypothesis tests for 108 regions. #Due to the variability from region to region, and the fact that some regions contain small counts, #we can use a spatial ICAR model for the local coefficients and SOI coefficient. #y[year,hexagon]~dpois(lambda[year,hexagon]) #log(lambda[year,hexagon])<-intercept[hexagon]+SST[year,hexagon]*sst[hexagon]+SOI[year]*soi[hexagon] #intercept[hexagon]<-int.u[hexagon]+int.s[hexagon] #soi[hexagon]<-soi.u[hexagon]+soi.s[hexagon] #sst[hexagon]<-sst.u[hexagon]+sst.s[hexagon] #soi.s[hexagon]~ICAR(tau_s_soi,neighborhood=hexagonnb,weights=1) #psuedo code car.normal(adj[],weights[],num[]) #sst.s[hexagon]~ICAR(tau_s_sst,...) #int.s[hexagon]~ICAR(tau,...) #soi.u~dnorm(int_soi,tau_soi) #sst.u~dnorm(int_sst,tau_sst) #... #The initial model with only the local SST and global SOI glm for initial values: #Initial data: WinBugsData<-c(list(H=nrow(Wind.Count),T=ncol(Wind.Count),N=as.matrix(Wind.Count),SST=as.matrix(SST.mean.year.subset),SOI=WinBugsFields$soi),hexagon.splus) writeDatafileR(WinBugsData,"../chap12/WinBugsData.txt") #Run WinBugs collect output codaout<-read.coda("../chap12/CODAchain1.txt","../chap12/CODAindex.txt") #Example for sst values #Collect sst terms sstsamples<-codaout[,substring(colnames(codaout),1,4)=="sst["] #For example .05,.10,.25,.5,.75,.90,.95 quantiles' sstquantiles<-t(apply(sstsamples,2,function(x) quantile(x,c(.05,.10,.25,.5,.75,.90,.95)))) rownames(sstquantiles)<-rownames(Wind.Count) #again sst.polygons<-SpatialPolygonsDataFrame(hpg[rownames(sstquantiles)],data.frame(sstquantiles))