#Bulk interpolation script in R

#The source data files must be in the same location as the script

#remember to set the working directory in R to the script location



library("automap")

library("sp")

library("gstat")

library("maptools")

library("rgdal")

library("raster")

library("RColorBrewer")



#set working directory

#setwd('E:/2013/Jirka/09')



projekce = "+proj=utm +zone=33"

variable = "snih"

aggregate = "max"

begin_date = as.Date("2013-01-01")

end_date = as.Date("2013-01-10")



#---COLOR RAMP SETTINGS (change as needed..)---

grid_colors = brewer.pal(6, "Blues")

intervaly = c(0,1,5,10,20,50,100)

intervaly2 = seq(1:length(intervaly))



#---READ ELEVATION GRID

srtm = raster('srtm_utm2.asc')

#(note): to get elevation for your area use: 

#

#---END READ ELEVATION GRID



#---CREATE DEFAULT EMPTY GRID

gt = getGridTopology(as(srtm, "SpatialGrid"))

defaultgrid = SpatialGrid(gt)

proj4string(defaultgrid) = CRS(projekce)

#---END CREATE DEFAULT EMPTY GRID



#---READ SHAPEFILE LAYERS

hranice = readShapeLines("hranice.shp")

proj4string(hranice) <- CRS("+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333")

hranice_proj <- spTransform(hranice, CRS(projekce))

vrstvahranice = list("sp.lines", hranice_proj, lwd=2.0,col="red")



toky = readShapeLines("toky.shp")

proj4string(toky) <- CRS("+proj=krovak +lat_0=49.5 +lon_0=24.83333333333333")

toky_proj <- spTransform(toky, CRS(projekce))

vrstvatoky = list("sp.lines", toky_proj, lwd=.6,col="blue")



#******************************

#***ADD YOUR SHAPEFILES HERE***

#

#******************************



#---END READ SHAPEFILE LAYERS



#---START THE FOR LOOP---

all_dates = seq(begin_date, end_date, 1)

for (j in 1:length(all_dates)) {

  selected_date = all_dates[j]

  

  #read data from the file with (id, lat, lon, elev, value) columns

  datapointfile = paste(variable, "_", aggregate, "_", format.Date(selected_date, "%Y-%m-%d"), ".txt", sep="")

  

  datapoints_wgs84 = read.table(datapointfile, header=TRUE)

  coordinates(datapoints_wgs84) = ~lon + lat

  proj4string(datapoints_wgs84) <- CRS("+init=epsg:4326")

  datapoints <- spTransform(datapoints_wgs84, CRS(projekce))



  #INTERPOLATION STARTS HERE!

  

  #calculate linear regression intercept and slope(observ=B+A*elev)

  observlm <- lm(value ~ elev, datapoints)

  datapoints$res = observlm$residuals

  

  #calculate observ value raster using linear model (elevation / observ value)

  intercept = observlm$coefficients[1]

  slope = observlm$coefficients[2]

  regression_grid <- intercept + srtm * slope

  

  #interpolate the residuals using idw

  idw_test = idw(res~1, datapoints, defaultgrid)

  residual_grid = raster(idw_test, "var1.pred")

  

  #add regression grid and residuals

  finalgrid <- regression_grid + residual_grid

  

  #INTERPOLATION ENDS HERE!



  

  #START OF THE PLOT!  

  layout = list(vrstvatoky, vrstvahranice)

  outfile = paste("obr", selected_date, ".png", sep="")

  png(filename = outfile, width = 1500, height = 1000, pointsize = 25, bg = "white", res = 150)

  nadpis = paste("snow", selected_date)

  print(spplot(finalgrid, at=intervaly, col.regions = grid_colors, 

               sp.layout=layout, 

               main=list(nadpis, cex=2, col="black", font=2), 

               colorkey=list(at=intervaly2, labels = list( at=intervaly2, cex = 1.5, labels = intervaly, lab = intervaly2), space="bottom")))

  dev.off()

  #END OF THE PLOT!

  

}

#---END THE FOR LOOP---