library(raster)
library(rasterVis)
library(maptools)
library(maps)
library(ncdf4)
library(ncdf.tools)
library(chron)
library(lattice)
library(RColorBrewer)

setwd("D:/...")


ncdf <- nc_open("2013FullYear.nc")

lon <- ncvar_get(ncdf,"Longitude")
nlon <- dim(lon)
head(lon)

lat <- ncvar_get(ncdf,"Latitude")
nlat <- dim(lat)
head(lat)

# Create a vector for Variable week
t <- ncvar_get(ncdf,"Step")
t

# Variable of interest in the NetCDF file
dname <- "Weekly Growth Index"

# Array with
tmp_array <- ncvar_get(ncdf,dname)

m <- 1
tmp_slice <- tmp_array[,,m]

# Create list with matrices for each time slot
tmp_stack <- vector("list",length(t))

for (i in 1:length(t)) {
        tmp_stack[[i]] <- tmp_array[,,i]
}

# Create list for 4 seasons:
tmp_stack_1 <- vector("list",length(t)/4)
for (i in 1:13) {
        tmp_stack_1[[i]] <- tmp_array[,,i]
}

tmp_stack_2 <- vector("list",length(t)/4)
for (i in 14:26) {
        for (j in 1:13){
                tmp_stack_2[[j]] <- tmp_array[,,i]
        }

}

tmp_stack_3 <- vector("list",length(t)/4)
for (i in 27:39) {
        for (j in 1:13){
                tmp_stack_3[[j]] <- tmp_array[,,i]
        }

}

tmp_stack_4 <- vector("list",length(t)/4)
for (i in 40:52) {
        for (j in 1:13){
                tmp_stack_4[[j]] <- tmp_array[,,i]
        }

}

# Summarize the seasons:
tmps_stack_avg1 <- Reduce("+",tmp_stack_1)/length(tmp_stack_1)
tmps_stack_avg2 <- Reduce("+",tmp_stack_2)/length(tmp_stack_2)
tmps_stack_avg3 <- Reduce("+",tmp_stack_3)/length(tmp_stack_3)
tmps_stack_avg4 <- Reduce("+",tmp_stack_4)/length(tmp_stack_4)

tmps_stack_max1 <- apply(simplify2array(tmp_stack_1), 1:2, max)
tmps_stack_max2 <- apply(simplify2array(tmp_stack_2), 1:2, max)
tmps_stack_max3 <- apply(simplify2array(tmp_stack_3), 1:2, max)
tmps_stack_max4 <- apply(simplify2array(tmp_stack_4), 1:2, max)



# Piece of code that gives nice map that looks right:
grid <- expand.grid(lon=lon, lat=lat)
cutpts <- seq(0,1,0.1)
levelplot(tmps_stack_max1 ~ lon * lat, data=grid, at=cutpts, cuts=11, pretty=T,
          col.regions=(rev(brewer.pal(10,"RdBu"))))

# Convert to raster attempt 1 - gives strange result!
raster1 <- raster(
        tmps_stack_max1,
        xmn=min(lon),xmx=max(lon),
        ymn=min(lat),ymx=max(lat),
        #crs=CRS(+proj=longlat +datum=WGS84)
        crs=CRS("+proj=cea  +lat_0=0 +lon_0=0")
)

plot(raster1)

#### Convert to raster... again wrong output!
r <- raster(tmps_stack_max1, xmn=min(lon), xmx=max(lon), ymn=min(lat), ymx=max(lat))

s <- as(r, 'SpatialGridDataFrame')
spplot(s)
writeGDAL(s, "SGDF.tif")