
# title         : getRADAR.R
# purpose       : formating and resampling of radar images (GIFs) from the Bilogora station;
# reference     : [http://spatial-analyst.net/book/Bilogora]
# producer      : Prepared by T. Hengl
# last update   : In Amsterdam, NL, May 2010.
# inputs        : gifs are publish on-line at [http://vrijeme.hr/aktpod.php?id=bradar&param=stat];
# outputs       : geotiff images projected in the Azimuthal Equidistant system;
# remarks 1     : Name of the hdf file "BIL_10_05_04_12_0000" is "BIL-year-month-day-hour-0000";
# remarks 2     : Obtain and install ImageMagick [http://www.imagemagick.org/script/download.php] and FWTools [http://fwtools.maptools.org/];

library(rgdal)
# library(caTools) # reading of GIF images (it had some problems!)
library(colorspace)
library(RSAGA)
AEQD <- "+proj=aeqd +lat_0=45.88356 +lon_0=17.20565"
AEQD2 <- "+proj=aeqd +a=6371000 +lat_0=45.88356 +lon_0=17.20565"  # SAGA GIS
# coordinates of the bounding box for GIFs (estimated by manually georeferencing image in ILWIS; based on country borders!):
bbox.gif <- data.frame(x=c(-200518,198316), y=c(-199993,198478)) # corner coordinates

# ------------------------------------------------------------
# list GIF files & prepare the legend
# ------------------------------------------------------------

gif.list <- dir(path=getwd(), pattern=glob2rx("RN1*.gif"), recursive=TRUE, full.names=TRUE)
# takes cca 1 minutes to list all files (8736);
gif.list[1:5]

# prepare local legend:
legend.y.pnt <- data.frame(y=24+c(7,22,37,52,66,81,97,112,128,143,158,173,188,203,218,233), x=rep(550.5,16), prec=c(.15,.35,.65,.9,1.5,3.5,6.5,9,15,35,65,90,150,350,650,950))
lm.y.legend <- glm(prec~y, legend.y.pnt, family=gaussian(log)) 

# read the legend file (this assumes that all legends are exactly the same):
gif.legend <- readGDAL(gif.list[i], offset=c(215, 550), region.dim=c(241, 1), silent=TRUE)
gif.legend$prec <- predict(lm.y.legend, gif.legend, type="response")
gif.legend$idx <- as.factor(gif.legend$band1)
gif.legend <- aggregate(gif.legend$prec, by=list(gif.legend$idx), FUN="mean");names(gif.legend) <- c("idx","prec")

# ------------------------------------------------------------
# read and convert to geotiff
# ------------------------------------------------------------

# get some idea about the gif files:
system("C:\\PROGRA~1\\IMAGEM~1\\identify -help")
gif.info <- system(paste("C:\\PROGRA~1\\IMAGEM~1\\identify -verbose", gif.list[1]), intern=TRUE) 
str(gif.info)
# color palette:
substr(gif.info[94:157], 31, 37)

options(warn=-1)  # there are some useless warning messages
# convert to geotiffs (takes cca 30 secs per image!)
for(i in length(gif.list):1){  # 2570:2592
 fname <- paste("RN1_", substr(gif.list[i], 36, 45), substr(gif.list[i], 57, 59), sep="")
 if(length(list.files(getwd(), pattern=set.file.extension(fname, ".tif")))==0) {
 # get color pallete (each GIF has another palette :(( ):
 # grid.cols <- data.frame(cols=read.gif(gif.list[i])$col)
 gif.info <- data.frame(cols=system(paste("C:\\PROGRA~1\\IMAGEM~1\\identify -verbose", gif.list[i]), intern=TRUE))
 grid.cols$idx <- as.factor(0:(length(grid.cols$cols)-1))
 # merge colors and prec values:
 grid.cols.prec <- merge(x=data.frame(gif.legend[c("idx", "prec")]), y=grid.cols, by=c("idx"), all.y=TRUE)
 # read image:
 gridgif <- readGDAL(gif.list[i], offset=c(1, 1), region.dim=c(478,478), silent=TRUE)
 # fix coordinates (needs to be corrected because these are the 240 km radius images):
 gridgif@bbox[1,1] <- bbox.gif[1,1]; gridgif@bbox[1,2] <- bbox.gif[2,1]; gridgif@bbox[2,1] <- bbox.gif[1,2]; gridgif@bbox[2,2] <- bbox.gif[2,2]
 gridgif@grid@cellsize <- c(833.35,833.35)
 gridgif@grid@cellcentre.offset[1] <- gridgif@bbox[1,1]+gridgif@grid@cellsize[1]/2; gridgif@grid@cellcentre.offset[2] <- gridgif@bbox[2,1]+gridgif@grid@cellsize[1]/2
 gridgif@coords[,1] <- gridgif@bbox[1,]+gridgif@grid@cellsize/2; gridgif@coords[,2] <- gridgif@bbox[2,]-gridgif@grid@cellsize/2
 # convert to prec values (match values from table):
 gridgif$idx <- as.factor(gridgif$band1)
 gridgif.pnt <- merge(data.frame(gridgif["idx"]), grid.cols.prec[c("idx", "prec")], by=c("idx"), all.x=TRUE)  # takes cca 10 secs!
 gridgif.pnt$prec <- round(ifelse(gridgif.pnt$prec<0, 0, gridgif.pnt$prec), 1)
 coordinates(gridgif.pnt) <- ~x+y
 gridded(gridgif.pnt) <- TRUE
 # write to GIS format:
 proj4string(gridgif.pnt) <- CRS("+proj=aeqd +lat_0=45.88356 +lon_0=17.20565")
 unlink("tmp.tif")
 writeGDAL(gridgif.pnt["prec"], "tmp.tif")
 # resample to UTM (AEQD is not really a popular CRS!)
 unlink(set.file.extension(fname, ".tif"))
 system(paste("C:\\PROGRA~1\\FWTOOL~1.7\\bin\\gdalwarp tmp.tif -s_srs \"+proj=aeqd +lat_0=45.88356 +lon_0=17.20565\" -t_srs \"+proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs\" -srcnodata \"-99999\" -dstnodata \"-99999\" ", set.file.extension(fname, ".tif"), " -r near -te 466500 4878500 874500 5286500 -tr 1000 1000", sep=""), show.output.on.console=FALSE) # takes cca 5 secs
}
}


# ------------------------------------------------------------
# derive daily sums
# ------------------------------------------------------------

# read mask maps:
bilogora <- readGDAL("mask_bilogora.asc")
bilogora$bmask <- readGDAL("bmask.asc")$band1 # black pixels in gif file!
bilogora$bmask <- ifelse(bilogora$bmask<0.15, NA, bilogora$bmask)
proj4string(bilogora) <- CRS("+proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
writeGDAL(bilogora["band1"], "mask_bilogora.sdat", "SAGA")

dates.list <- as.Date("2007-12-31")+0:367
for(j in 1:length(dates.list)){
  if(length(list.files(getwd(), pattern=paste("RN1_", gsub("-", "_", dates.list[j]), ".tif", sep="")))==0) {
  day.list <- dir(path=getwd(), pattern=glob2rx(paste("RN1_", gsub("-", "_", dates.list[j]), "_**.tif", sep="")), recursive=FALSE, full.names=FALSE)
  if(!length(day.list)==0){
  for(i in 1:length(day.list)){ 
  bilogora@data[,paste("h", i, sep="")] <- readGDAL(day.list[i], silent=TRUE)$band1
  # replace NA values with 0 values! 
  bilogora@data[,paste("h", i, sep="")] <- ifelse(is.na(bilogora$band1), NA, ifelse(is.na(bilogora@data[,paste("h", i, sep="")])&is.na(bilogora$bmask), 0,  bilogora@data[,paste("h", i, sep="")]))
 } # takes 10 secs; 
 # derive summary precipitation:
 bilogora$m <- rowSums(bilogora@data[,paste("h", 1:i, sep="")], na.rm=FALSE, dims=1)
 writeGDAL(bilogora["m"], "m.sdat", "SAGA")
 # filter missing pixels (close gap) using a mask map:
 rsaga.geoprocessor(lib="grid_tools", module=7, param=list(INPUT="m.sgrd", RESULT="mf.sgrd", MASK="mask_bilogora.sgrd"), show.output.on.console=FALSE)
 system(paste("C:\\PROGRA~1\\FWTOOL~1.7\\bin\\gdal_translate mf.sdat", paste("RN1_", gsub("-", "_", dates.list[j]), ".tif", sep=""), "-a_srs \"+proj=utm +zone=33 +ellps=WGS84 +datum=WGS84 +units=m +no_defs\""), show.output.on.console=FALSE) 
 }}
}

# create an animation (gif) using a range of gifs:
setwd(paste(getwd(), "2008_05_01", sep="/"))
system("C:\\PROGRA~1\\IMAGEM~1\\convert -delay 20 -loop 0 RN1_Lq_000_*.gif  RN1_2008_05_01.gif")

# end of script;
