library(ncdf)
library(sp)
library(gstat)
files <- dir("C:/Users/Sasha/Desktop/f", recursive=TRUE, full.names=TRUE, pattern="\\.nc$")
# "C:/Users/Sasha/Desktop/f" in previous row. It is path where stored my data in the following way:
# year1966.nc
# year1967.nc
# year1968.nc...till year2014.nc (49 years)
# creating function for processing these files
processFile <- function(f) {
nc <- open.ncdf(f) # opening files
startdate <- basename(f) # get file name ("year1987.nc" for example).
unlist(strsplit(startdate, "[.]"))->x # separating file name ("year1987" "nc")
x[1]->x # take just first part ("year1987")
unlist(strsplit(x, "year"))->x # separating "year1987" (""     "1987")
unlist(strsplit(x, "year"))->x # repeat previous row in order to get only "1987"
paste(x,"01-01", sep="-")-> startday # create startday ("1987-01-01")
paste(x,"12-31", sep="-")-> endday # create endday ("1987-12-31")
dates = seq(as.Date(startday), as.Date(endday), by="1 day") # generate dates for current year
dates -> date # date will need for creating file lonlatdate (29 rows); dates - for interpolation (47 rows)
rotation <- length(date) # set the value for further loop. This is important because there are leap years
print(nc) # check attributes of nc file
lon <- get.var.ncdf(nc,"lon") # get longitude from nc file
lat <- get.var.ncdf(nc,"lat") # get latitude from nc file
var <- get.var.ncdf(nc,"tasmax") # get values of max temperature from nc file
merge(lon,lat)->lonlat # merge coordinates
as.data.frame(lonlat)->lonlat
as.data.frame(date)->date
merge(lonlat,date)-> lonlatdate # merge coordinates with dates
as.vector(var)->tasmax
cbind(lonlatdate,tasmax)->data # add values of max temperature to lonlatdates
data[complete.cases(data), ]->out # get rid of rows which contain NA
# Now we have array(out)with max temperature for each grid per each day of 1987 year:
out[1,]
         x      y       date   tasmax
120 33.625 44.625 1987-01-01 279.2693
station <- read.table("C:\\Users\\sasha\\Desktop\\station.txt",sep='',header=T) # download stations data
station[1,] # how it looks like
    wmo     x     y
1 33998 34.08 44.43
# And now we must to interpolate tasmax values from out to station
# but before creating file for results (res)
res <- c(1:5)
dim(res)<-c(1,5)
res[-1,]->res
colnames(res)<- c("wmo","x","y","tasmax","date") # setting names for columns of res file
#interpolation
for (i in seq(rotation)){
  subset(out, date==dates[i])->test;station->test2;coordinates(test)<-~x+y;coordinates(test2)<-~x+y;
  idw(tasmax~1,test,test2)->z;zdf<- as.data.frame(z);cbind(station$wmo,zdf)->me;me[5]<-dates[i];
  colnames(me)<- c("wmo","x","y","tasmax","date");rbind(res,me)->res} # test,test2,z,zdf,me - these are temporary files
as.Date(res$date,origin="1970-01-01")-> res$date # convert values in res$date to Date format
# This is what we have:
head(res)
    wmo     x     y   tasmax       date
1 33998 34.08 44.43 275.6095 1987-01-01
2 33959 34.43 44.68 276.8527 1987-01-01
3 34622 38.51 47.80 273.2416 1987-01-01
4 33958 34.35 44.75 275.8180 1987-01-01
5 34510 37.98 48.60 273.4068 1987-01-01
6 33915 33.88 46.45 276.8582 1987-01-01
write.table(res,f,row.names = FALSE) # will rewriting all nc files in folder f
file.info(f)$size
}
result <- sapply(files,processFile) # apply processFile function to all nc files in folder f