rm(list=ls())


library(ncdf4)
library(raster)


# http://www.unidata.ucar.edu/software/netcdf/docs/faq.html#notthreadsafe


setwd("N:\\Plumes_interpolointi_201603_Monthly\\Snow_FI_20160304\\")
polku_IN <- c("N:\\Plumes_interpolointi_201603_Monthly\\Snow_FI_20160304\\NetCDF\\")
polku_OUT <- c("N:\\Plumes_interpolointi_201603_Monthly\\Test\\")


hila <- read.csv("N:\\Plumes_interpolointi_201603_Monthly\\hila_Fi.csv")
hila_Fi <-  hila[order(hila$point),]

point_count <- dim(hila)[1]  # number of rows in the dataFrame
zero_date <- as.Date(as.character(17991231),format='%Y%m%d') #

x_vals <- hila_Fi$x
y_vals <- hila_Fi$y

vv <- 2010
paiva <- 20
ncfname <- paste(polku_IN,"Snow_",vv,".nc",sep="")
tiffname <- paste(polku_OUT,"Snow_",vv,".nc",sep="")

 
  #  title <- ncatt_put(ncout,0,paste("Daily mean temperature/ year ",vv,sep=""))
  #  institution <- ncatt_put(ncout,0,"Finnish Meteorological Institute FMI")
  #  datasource <- ncatt_put(ncout,0,"Interpolated data")
  #  references <- ncatt_put(ncout,0,"Journal of Geophysical Research: New gridded daily climatology of Finland: Permutation-based uncertainty estimates and temporal trends in climate")
  #  history <- ncatt_put(ncout,0,"Version 1 2016/05")



ncin <- nc_open(ncfname)
print(ncin)

points <- ncvar_get(ncin,"point")
aika <- ncvar_get(ncin,"Time")

head(points)
head(aika)


n_points <- dim(points)
n_aika <- dim(aika)


z1 <- ncvar_get(ncin,"SD",start=c(1,paiva),count=c(-1,1))
pvm <- as.character(as.Date(aika[paiva],origin=zero_date),format="%Y%m%d")

x <- ncvar_get(ncin,"longitude",start=1,count=-1)
y <- ncvar_get(ncin,"latitude",start=1,count=-1)

testi <- as.data.frame(x)
testi$y <- y
testi$Snow <- z1


coordinates(testi) <- c("x", "y")
gridded(testi) <- T
raster_layer <- raster(testi)
projection(raster_layer) <-CRS("+init=epsg:3067")

file_out<- paste(polku_OUT , "Snow_",pvm,".tif", sep="") 
writeRaster(raster_layer, filename=file_out,overwrite=TRUE)



nc_close(ncin)






