library(xtable)

library(ellipse) # for correlation matrix

library(maptools) # for rgshhs

## load coastline data

roill=bbox(projectExtent(topo,CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")))

coast=getRgshhsMap("/home/adamw/acrobates/Global/GSHHS_shp/gshhs/gshhs_f.b",xlim=c(roill[1,1],roill[1,2]),ylim=c(roill[2,1],roill[2,2]))

coast=as(as(coast,"SpatialLines"),"SpatialLinesDataFrame")

coast@data=data.frame(id=1:length(coast@lines))  #convert to dataframe

coast=spTransform(coast,projs)

rcoast=rasterize(coast,topo)

dist=distance(rcoast)/1000 #get distance to coast and convert to km

### transform distance to coast using b-log(dist), where b \approx log(DTC of the farthest point on earth from the sea)

### Garcia-Castellanos, Daniel, and Umberto Lombardo. 2007. “Poles of Inaccessibility: A Calculation Algorithm for the Remotest Places on Earth.” Scottish Geographical Journal 123 (3): 227–233. doi:10.1080/14702540801897809.

### farthest point is ~2510km from ocean

log(2510) #round up to 8 to be sure we won't get negative numbers

dist2=brick(dist,calc(dist,function(x) (8-log(x+1))))

layerNames(dist2)=c("DTCkm","logDTCkm")

colnames(dist2@data@values)=c("DTCkm","logDTCkm")

topo2=stack(subset(topo,"dem"),EW,NS,slope,dist2)

##topo=brick(list(topo,demb))

st2=SpatialPointsDataFrame(st2,data=cbind.data.frame(st2@data,demb=extract(demb,st2),extract(topo,st2),extract(topo2,st2),extract(lulc,st2),extract(lst,st2)))

s=SpatialPointsDataFrame(s,data=cbind.data.frame(x=coordinates(s)[,1],y=coordinates(s)[,2],demb=extract(demb,s),extract(topo,s),extract(topo2,s),extract(lulc,s),extract(lst,s)))

save(s,file=paste("output/mod_sample.Rdata",sep=""))

Sys.setenv(MKL_NUM_THREADS=24)

system("export MKL_NUM_THREADS=24")

niter=12500

lapply(months,function(m){

  f1=formula(paste(m,"~logDTCkm+y+Shrub+Grass+Crop+Mosaic+Urban+Barren+Snow+Wetland+dem+eastwest+northsouth",sep=""))

    n.samples=niter,sub.samples=c(2500,niter,10),verbose=TRUE, n.report=100)

### Read in results

load(paste("output/mod_sample.Rdata",sep=""))

ms=lapply(months,function(m) {print(m); load(paste("output/mod_",m,".Rdata",sep="")) ; get(paste("mod_",m,sep="")) })

### generate summaries

ms1=do.call(rbind.data.frame,lapply(1:length(ms),function(i){

  mi=ms[[i]]

  m1.s=mcmc(t(mi$p.samples),start=round(niter/4),thin=1)

  m1.s2=as.data.frame(t(apply(m1.s,1,quantile,c(0.025,0.5,0.975))))

  colnames(m1.s2)=c("Q2.5","Q50","Q97.5")

  m1.s2$parm=rownames(m1.s)

  m1.s2$month=factor(months[i],levels=months,ordered=T)

  m1.s2$type=ifelse(m1.s2$parm%in%c("tau.sq","sigma.sq","phi"),"Spatial",ifelse(m1.s2$parm%in%c("(Intercept)","eastwest","northsouth","dem","logDTCkm","y"),"Topography","LULC"))

  return(m1.s2)

}))

## drop spatial parameters

#ms1=ms1[!ms1$type%in%"Spatial",]

## Convert dem to degrees/km to make it comparable to other topographical parameters

#ms1[ms1$parm%in%c("x","y","dem"),c("Q2.5","Q50","Q97.5")]=ms1[ms1$parm%in%c("x","y","dem"),c("Q2.5","Q50","Q97.5")]

## DEM

at=unique(quantile(as.matrix(subset(topo2,c("dem"))),seq(0,1,length=100),na.rm=T))

levelplot(subset(topo2,c("dem")),at=at,col.regions=bgyr(length(at)),

          main="Elevation",sub="",margin=F)+

  layer(sp.lines(roi, lwd=1.2, col='black'))

## DTC

at=unique(quantile(as.matrix(subset(topo2,c("logDTCkm"))),seq(0,1,length=100),na.rm=T))

levelplot(subset(topo2,c("logDTCkm")),at=at,col.regions=bgyr(length(at)),

          main="Elevation",sub="",margin=F)+

  layer(sp.lines(roi, lwd=1.2, col='black'))

       sub="One point per pixel, per month",ylab="LST",

## just forest, grouped by distance to coast

xyplot(value~dem|variable,groups=cut(DTCkm,breaks=c(0,300,1000)),melt(s@data[s$lulc=="Forest",c("DTCkm","dem","lulc",months)],id.vars=c("DTCkm","dem","lulc")),pch=16,cex=.5,

       main="Month-by-month scatterplots of Elevation and LST, grouped by Distance to Coast",

       sub="Showing only forest class",ylab="LST",

       par.settings = list(superpose.symbol = list(pch=16,cex=.5)),auto.key=list(space="right",title="Distance to \n Coast (km)"))+

  layer(panel.abline(lm(y~x),col="red"))+

  layer(panel.text(500,50,bquote(beta==.(round(coefficients(lm(y~x))[2],4)))))

useOuterStrips(xyplot(value~dem|variable+lulc,groups=cut(DTCkm,breaks=c(0,300,1000)),

                      melt(s@data[s$lulc!="Snow",c("DTCkm","dem","lulc",months)],id.vars=c("DTCkm","dem","lulc")),pch=16,cex=.5,

       main="Month-by-month scatterplots of Elevation and LST, grouped by LULC and colored by distance to coast",

       sub="One point per pixel, per month",ylab="LST",

       par.settings = list(superpose.symbol = list(pch=16,cex=.5)),auto.key=list(space="right",title="Distance to \n Coast (km)")))+

  layer(panel.abline(lm(y~x),col="red"))

##########################################

### Model output

### show fitted values

spplot(s,zcol="demb",col.regions=terrain.colors(5),cex=.5)+layer(sp.lines(roi,lwd=1.2,col="black"))+

  layer(sp.points(s2,col="blue",pch=13,cex=2,lwd=2))+

  layer(sp.points(st2,col="black",lwd=1.5))

## look at correlation of LULC variables

lulct=s@data[,colnames(s@data)[colnames(s@data)%in%unique(s$lulc)]]

lulcc=cor(lulct)[order(cor(lulct)[1,]),order(cor(lulct)[1,])]

round(cor(lulct),2)

#print(xtable(round(cor(lulct),2),include.row.names=F))

plotcorr(lulcc,col=colors <- c("#A50F15","#DE2D26","#FB6A4A","#FCAE91","#FEE5D9","white","#EFF3FF","#BDD7E7","#6BAED6","#3182BD","#08519C")[5*lulcc+6],

         main="Correlation Matrix")

t1=as.data.frame(table(s$lulc));colnames(t1)=c("class","samplefreq")

t1$sampleprop=t1$samplefreq/sum(t1$samplefreq)

t2=as.data.frame(table(st2$lulc))

t1$stationfreq=0

t1$stationfreq[t1$class%in%t2$Var1]=t2$Freq[match(t2$Var1,t1$class[t1$class%in%t2$Var1])]

t1$stationprop=t1$stationfreq/sum(t1$stationfreq)

print(xtable(t1[,c(1,3,5)]),include.rownames=F)

### two plots of model parameters by month to show the effects of LULC and topography on LST

## Spatial

xyplot(Q50~month|parm,data=ms1[ms1$type=="Spatial",],panel=function(x,y,subscripts){

  tms1=ms1[ms1$type=="Spatial",][subscripts,]

  panel.segments(tms1$month,tms1$Q2.5,tms1$month,tms1$Q97.5,groups=tms1$parm,subscripts=1:nrow(tms1))

  panel.xyplot(tms1$month,tms1$Q50,pch=16,type="l")

  panel.abline(h=0,lty="dashed",col="grey")

},auto.key=list(space="right"),par.settings = list(superpose.symbol = list(pch=16,cex=.5,col=1:15),superpose.line=list(col=1:15)),

       subscripts=T,scales=list(y=list(relation="free",rot=90)),ylab="Parameter Coefficient (slope) with 95% Credible Intervals",xlab="Month",

       main="Spatial Parameters",

       sub="")

## LULC

xyplot(Q50~month,groups=parm,data=ms1[ms1$type=="LULC",],panel=function(x,y,subscripts){

  tms1=ms1[ms1$type=="LULC",][subscripts,]

  sig=ifelse(tms1$Q2.5<0&tms1$Q97.5<0|tms1$Q2.5>0&tms1$Q97.5>0,"red","black")

  panel.segments(tms1$month,tms1$Q2.5,tms1$month,tms1$Q97.5,groups=tms1$parm,subscripts=1:nrow(tms1))

  panel.xyplot(tms1$month,tms1$Q50,groups=tms1$parm,pch=16,type="l",subscripts=1:nrow(tms1))

  panel.abline(h=0,lty="dashed",col="grey")

},auto.key=list(space="right"),par.settings = list(superpose.symbol = list(pch=16,cex=.5,col=1:15),superpose.line=list(col=1:15)),

       subscripts=T,scales=list(y=list(relation="free",rot=90)),ylab="Parameter Coefficient (slope) with 95% Credible Intervals",xlab="Month",

       main="Effects of LULC on LST",

       sub="Coefficients are unstandardized and represent the change in LST expected with a 1% increase in that class from 100% Forest")

## Topography

xyplot(Q50~month|parm,data=ms1[ms1$type=="Topography",],panel=function(x,y,subscripts){

  tms1=ms1[ms1$type=="Topography",][subscripts,]

  panel.segments(tms1$month,tms1$Q2.5,tms1$month,tms1$Q97.5,subscripts=1:nrow(tms1))

  panel.xyplot(tms1$month,tms1$Q50,groups=tms1$parm,pch=16,type="l",subscripts=1:nrow(tms1))

  panel.abline(h=0,lty="dashed",col="grey")

},auto.key=list(space="right"),par.settings = list(superpose.symbol = list(pch=16,cex=.5,col=1:15),superpose.line=list(col=1:15)),

       subscripts=T,scales=list(y=list(relation="free",rot=90)),ylab="Parameter Coefficient (slope) with 95% Credible Intervals",xlab="Month",

       main="Effects of Topography on LST",

       sub="Coefficients are unstandardized. Intercept is degrees C, eastwest and northsouth range \n from -1 (90 degree slope) to 0 (flat) to 1 (90 degree slope), and dem is m \n logDTCkm is in log(km), and y (lat) is m")

### Capture same pattern using only station data?

dev.off()

 shrinkpdf<-function(pdf,maxsize=1,suffix="_small",verbose=T){

   require(multicore)

   wd=getwd()

   td=paste(tempdir(),"/pdf",sep="")

   if(!file.exists(td)) dir.create(td)

   if(verbose) print("Performing initial compression")

   setwd(td)

   if(verbose) print("Splitting the PDF to parallelize the processing")

   system(paste("pdftk ",wd,"/",pdf," burst",sep=""))

   mclapply(list.files(td,pattern="pdf$"),function(f){

     ## loop through all pages, perform compression with with ps2pdf

     if(verbose) print("Performing initial compression")

       system(paste("ps2pdf -dUseFlateCompression=true ",td,"/",f," ",td,"/compressed_",f,sep=""))

        file.rename(paste(td,"/compressed_",f,sep=""),paste(td,"/",f,sep=""))

     ## get sysmte size

        size=file.info(paste(td,"/",f,sep=""))$size*0.000001 #get sizes of individual pages

        toobig=size>=maxsize

        if(verbose&toobig)  print(paste("Resizing ",basename(paste(td,"/",f,sep="")),sep=""))

        system(paste("gs -dBATCH -dTextAlphaBits=4 -dNOPAUSE -r300 -q -sDEVICE=png16m -sOutputFile=- -q ",paste(td,"/",f,sep="")," | convert -background transparent -quality 100 -density 300 - ",f,sep=""))

      })

        if(verbose) print("Compiling the final pdf")

        setwd(wd)

        system(paste("gs -dBATCH -dNOPAUSE -q -sDEVICE=pdfwrite -sOutputFile=",strsplit(pdf,".",fixed=T)[[1]][1],suffix,".pdf ",td,"/*.pdf",sep=""))

        file.remove(list.files(td,full=T))

       if(verbose) print("Finished!!")

}

shrinkpdf("output/lst_lulc.pdf")