# Filename: plot_eco_av.R
# This file serves to plot the variables contained in the files eco_av.YYYYMM.nc. with R
# 
# Script written by Eelke Folmer (Eelke.Folmer@nioz.nl) on the basis of a python script by 
# Matias Duran Matute (matias.duran-matute@nioz.nl; m.duran.matute@tue.nl)
#
#'This script is made available in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.'
library(ncdf4)
library(raster)
library(rgdal)
library(gdalUtils)

# load the data
nc.file <- nc_open('/home/eelke/work/waddenzee/abiotics/matias/eco_av.200901.nc')
var    <- ncvar_get(nc.file, varid='salt_surf')
lat    <- ncvar_get(nc.file, varid='lat')
lon    <- ncvar_get(nc.file, varid='lon')

# make it a SpatialPointsDataFrame
df           <- data.frame(x=c(lon), y=c(lat), z=c(var))
df           <- df[-which(is.na(df$x) | is.na(df$y) | is.na(df$z) ), ]
df.spdf.wgs  <- SpatialPointsDataFrame(cbind(df$x, df$y), df["z"], proj4string =  CRS("+proj=longlat +datum=WGS84"), match.ID = TRUE, bbox = NULL)
spplot(df.spdf.wgs)

# make it a raster in rd (epsg=28992); due to the geographic transformation of a vector gaps emerge if the resolution is taken too small
# an alternative would be to interpolate as an intermediate step
# rasters can be written to disk as geotif so that they can be loaded in a GIS
df.spdf.rd  <- spTransform(df.spdf.wgs, CRS("+init=epsg:28992") )
numcells.x  <- round( (extent(df.spdf.rd)@xmax - extent(df.spdf.rd)@xmin ) / 250 )
numcells.y  <- round( (extent(df.spdf.rd)@ymax - extent(df.spdf.rd)@ymin ) / 250 )
rt          <- raster(ncols=numcells.x, nrows=numcells.y, crs=("+init=epsg:28992"), ext=extent(df.spdf.rd))
r           <- rasterize(df.spdf.rd, rt, df.spdf.rd$z, fun=mean)
plot(r2)