library(mapdata)

library(maptools)

library(raster)

library(sp)

library(dplyr)

library(rgeos)

library(foreach)

## Get all available tiles

tiles=list.files("/nobackupp6/aguzman4/climateLayers/out/series_reg5",recursive=T,pattern="gam_CAI_dailyTmax",full=T)

## reconstruct tile/date metadata

tdata=data.frame(

  path=tiles)

tdata[,c("form","type","var","pred","model","x1","x2","date","x3","lat","lon")]=

  do.call(rbind.data.frame,strsplit(sub("[.]tif","",basename(tiles)),split="_"))

tdata$tile=paste0(tdata$lat,"_",tdata$lon)

tdata$year=as.numeric(substr(tdata$date,1,4))

## Define domain

#data(worldMapEnv)

#wld=map('world' , fill = T,plot=F,xlim=c(-20,55),ylim=c(-35,40))

#wlds=map2SpatialPolygons(wld,IDs= sapply(strsplit(wld$names, ":"), function(x) x[1]))

#projection(wlds)="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

# Get domain of tiles

utile=unique(tdata$tile)

stile=do.call(rbind,lapply(1:length(utile),function(u) tdata[tdata$tile==utile[u],][1,]))

t1=as(extent(raster(as.character(stile$path[1]))), 'SpatialPolygons')

t2=as(extent(raster(as.character(stile$path[2]))), 'SpatialPolygons')

domain=gUnion(t1,t2)

## Select locations

points=spsample(domain,100,type="random")

write.csv(as.data.frame(points),"~/interp/samplePoints.csv")

## Extract timeseries across tifs

#library(doParallel)

#registerDoParallel(2)

ytiles=unique(data.frame(year=tdata$year,tile=tdata$tile))

foreach(i=1:nrow(ytiles)) %do% {

  tyear=ytiles$year[i]

  ttile=ytiles$tile[i]

   outfile=paste0("~/interp/pdata_",ttile,"_",tyear,".csv")

    if(!file.exists(outfile)) {

    rs=stack(as.character(tdata$path[tdata$year==tyear&tdata$tile==ttile]))

    names(rs)=tdata$date[tdata$year==tyear&tdata$tile==ttile]

    pdata=data.frame(tile=ttile,year=tyear,as.data.frame(points),extract(rs,points))

    write.csv(pdata,outfile)

    print(paste(ttile," ",tyear))

}

}

## Make animation

foreach(i=1:nrow(ytiles)) %do% {

  tyear=ytiles$year[i]

  ttile=ytiles$tile[i]

   outfile=paste0("~/interp/pdata_",ttile,"_",tyear,".png")

    if(!file.exists(outfile)) {

    rs=stack(as.character(tdata$path[tdata$year==tyear&tdata$tile==ttile]))

    names(rs)=tdata$date[tdata$year==tyear&tdata$tile==ttile]

    rs@z=list(as.Date(substr(names(rs),2,9),"%Y%m%d"))

    pdata=data.frame(tile=ttile,year=tyear,as.data.frame(points),extract(rs,points))

    write.csv(pdata,outfile)

    print(paste(ttile," ",tyear))

}

}

#animate(rs)