Translate a net-cdf file with different years to polygon regions.

*DATA*

I have a net-cdf file which has raster of 0.5 on 0.5 degrees. You can
easily download it here
<https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_3.23/cruts.1506241137.v3.23/pet/>
and
search for *cru_ts3.23.2001.2010.pet.dat.nc.gz *(it is also downloadable as
a dat.file if this is more handy to work with)
 (Or simply download the net-cdf file directly through:
cru_ts3.23.2001.2010.pet.dat.nc.gz
<https://crudata.uea.ac.uk/cru/data/hrg/cru_ts_3.23/cruts.1506241137.v3.23/pet/cru_ts3.23.2001.2010.pet.dat.nc.gz>
).

I opened the file in ArcMap as well and found that the coordinate system
used is: GCS_WGS_1984. The net-cdf file contains monthly data from 2001-2010

Download from this website
<http://ec.europa.eu/eurostat/web/gisco/geodata/reference-data/administrative-units-statistical-units/nuts>
the
following ZIP file: *NUTS_2006_60M_SH.zip* and save all the *NUTS_RG_60M_2006
files *in a folder where you can easily find the back. In what follows I
refer to this as "NUTS3".

*WHAT I WANT*

- I want to add the information from the raster files to the NUTS3
shapefile (which are polygons and no rasters) in order to obtain a table
per nuts3 region for each monthtly variable.

*WHERE I AM STUCK*

The file appears to be very difficult to work with in ArcGis. Also, I will
have to repeat this process a number of times for different variables. So I
would like to have one code in R that I can run.

I have a number of separate codes to do the following:
- translate the net-cdf file to a cvs file with longitude and latitude as
identifiers (see below under 1)
- translate a cvs file with accompanying empty raster file to nuts3
regions. (this is a code that I have used before when I had a cvs file and
a raster). (see below under 2).

However, I don't have a raster file now. Well, technically I do (since the
net-ncf file is a raster) but I don't know how to use it in this format.

Can somebody help me to link the codes below or suggest a different code to
obtain what I want?

Thanks a lot!

Janka

*1) With the following code, I can make a cvs file and extract all the data
in table format.*

library(fields)
library(chron)

library(ncdf4)
ncname<-"cru_ts3.22.2001.2010.pet.dat"
ncname<-"cru_ts3.23.1991.2000.pet.dat"
ncfname <- paste(ncname, ".nc", sep = "")
dname <- "pet"
ncin <- nc_open(ncfname)
print(ncin)

lon <- ncvar_get(ncin, "lon")
nlon <- dim(lon)
head(lon)

lat <- ncvar_get(ncin, "lat", verbose = F)
nlat <- dim(lat)
head(lat)

print(c(nlon, nlat))

t <- ncvar_get(ncin, "time")
tunits <- ncatt_get(ncin, "time", "units")
nt <- dim(t)

tmp.array <- ncvar_get(ncin, dname)
dlname <- ncatt_get(ncin, dname, "long_name")
dunits <- ncatt_get(ncin, dname, "units")
fillvalue <- ncatt_get(ncin, dname, "_FillValue")
dim(tmp.array)

title <- ncatt_get(ncin, 0, "title")
institution <- ncatt_get(ncin, 0, "institution")
datasource <- ncatt_get(ncin, 0, "source")
references <- ncatt_get(ncin, 0, "references")
history <- ncatt_get(ncin, 0, "history")
Conventions <- ncatt_get(ncin, 0, "Conventions")

nc_close(ncin)

# split the time units string into fields
tustr <- strsplit(tunits$value, " ")
tdstr <- strsplit(unlist(tustr)[3], "-")
tmonth = as.integer(unlist(tdstr)[2])
tday = as.integer(unlist(tdstr)[3])
tyear = as.integer(unlist(tdstr)[1])
chron(t, origin = c(tmonth, tday, tyear))

tmp.array[tmp.array == fillvalue$value] <- NA

length(na.omit(as.vector(tmp.array[, , 1])))

m <- 1
tmp.slice <- tmp.array[, , m]
library(RColorBrewer)
image(lon, lat, tmp.slice, col = rev(brewer.pal(10, "RdBu")))

grid <- expand.grid(lon = lon, lat = lat)
cutpts <- c(-50, -40, -30, -20, -10, 0, 10, 20, 30, 40, 50)
levelplot(tmp.slice ~ lon * lat, data = grid, at = cutpts, cuts = 11,
pretty = T,
          col.regions = (rev(brewer.pal(10, "RdBu"))))

lonlat <- expand.grid(lon, lat)
tmp.vec <- as.vector(tmp.slice)
length(tmp.vec)

tmp.df01 <- data.frame(cbind(lonlat, tmp.vec))
names(tmp.df01) <- c("lon", "lat", paste(dname, as.character(m), sep = "_"))
head(na.omit(tmp.df01), 20)

csvfile <- "cru_tmp_1.csv"
write.table(na.omit(tmp.df01), csvfile, row.names = FALSE, sep = ",")

tmp.vec.long <- as.vector(tmp.array)
length(tmp.vec.long)

tmp.mat <- matrix(tmp.vec.long, nrow = nlon * nlat, ncol = nt)
dim(tmp.mat)

head(na.omit(tmp.mat))

lonlat <- expand.grid(lon, lat)
tmp.df02 <- data.frame(cbind(lonlat, tmp.mat))

names(tmp.df02) <- c("lon","lat","pet_jan_2001",
                     "pet_feb_2001",
                     "pet_mar_2001",
                     "pet_apr_2001",
                     "pet_may_2001",
                     "pet_jun_2001",
                     "pet_jul_2001",
                     "pet_aug_2001",
                     "pet_sep_2001",
                     "pet_oct_2001",
                     "pet_nov_2001",
                     "pet_dec_2001",
                     "pet_jan_2002",
                     "pet_feb_2002",
                     "pet_mar_2002",
                     "pet_apr_2002",
                     "pet_may_2002",
                     "pet_jun_2002",
                     "pet_jul_2002",
                     "pet_aug_2002",
                     "pet_sep_2002",
                     "pet_oct_2002",
                     "pet_nov_2002",
                     "pet_dec_2002",
                     "pet_jan_2003",
                     "pet_feb_2003",
                     "pet_mar_2003",
                     "pet_apr_2003",
                     "pet_may_2003",
                     "pet_jun_2003",
                     "pet_jul_2003",
                     "pet_aug_2003",
                     "pet_sep_2003",
                     "pet_oct_2003",
                     "pet_nov_2003",
                     "pet_dec_2003",
                     "pet_jan_2004",
                     "pet_feb_2004",
                     "pet_mar_2004",
                     "pet_apr_2004",
                     "pet_may_2004",
                     "pet_jun_2004",
                     "pet_jul_2004",
                     "pet_aug_2004",
                     "pet_sep_2004",
                     "pet_oct_2004",
                     "pet_nov_2004",
                     "pet_dec_2004",
                     "pet_jan_2005",
                     "pet_feb_2005",
                     "pet_mar_2005",
                     "pet_apr_2005",
                     "pet_may_2005",
                     "pet_jun_2005",
                     "pet_jul_2005",
                     "pet_aug_2005",
                     "pet_sep_2005",
                     "pet_oct_2005",
                     "pet_nov_2005",
                     "pet_dec_2005",
                     "pet_jan_2006",
                     "pet_feb_2006",
                     "pet_mar_2006",
                     "pet_apr_2006",
                     "pet_may_2006",
                     "pet_jun_2006",
                     "pet_jul_2006",
                     "pet_aug_2006",
                     "pet_sep_2006",
                     "pet_oct_2006",
                     "pet_nov_2006",
                     "pet_dec_2006",
                     "pet_jan_2007",
                     "pet_feb_2007",
                     "pet_mar_2007",
                     "pet_apr_2007",
                     "pet_may_2007",
                     "pet_jun_2007",
                     "pet_jul_2007",
                     "pet_aug_2007",
                     "pet_sep_2007",
                     "pet_oct_2007",
                     "pet_nov_2007",
                     "pet_dec_2007",
                     "pet_jan_2008",
                     "pet_feb_2008",
                     "pet_mar_2008",
                     "pet_apr_2008",
                     "pet_may_2008",
                     "pet_jun_2008",
                     "pet_jul_2008",
                     "pet_aug_2008",
                     "pet_sep_2008",
                     "pet_oct_2008",
                     "pet_nov_2008",
                     "pet_dec_2008",
                     "pet_jan_2009",
                     "pet_feb_2009",
                     "pet_mar_2009",
                     "pet_apr_2009",
                     "pet_may_2009",
                     "pet_jun_2009",
                     "pet_jul_2009",
                     "pet_aug_2009",
                     "pet_sep_2009",
                     "pet_oct_2009",
                     "pet_nov_2009",
                     "pet_dec_2009",
                     "pet_jan_2010",
                     "pet_feb_2010",
                     "pet_mar_2010",
                     "pet_apr_2010",
                     "pet_may_2010",
                     "pet_jun_2010",
                     "pet_jul_2010",
                     "pet_aug_2010",
                     "pet_sep_2010",
                     "pet_oct_2010",
                     "pet_nov_2010",
                     "pet_dec_2010")

options(width = 110)
head(na.omit(tmp.df02, 20))

dim(na.omit(tmp.df02))

csvfile <- "cru_tmp_2.csv"
write.table(na.omit(tmp.df02), csvfile, row.names = FALSE, sep = ",")

*2) translate a cvs-file with accompanying raster file to polygon regions.*

The  "filename.txt" file should contain the variables: lon, latitude, and
all the monthly_yearly variables extracted from point 1 above.

The grid shapefile (*grid_025dd.shp*) can be found through the following
link but it is only an example and not the correct grid for the problem
above :
https://drive.google.com/folderview?id=0By9u5m3kxn9yfjZtdFZLcW82SWpzT1VwZXE1a3FtRGtSdEl1c1NvY205TGpack9xSFc2T2s&usp=sharing

# upload data
mydata<-read.table("filename.txt", header=TRUE,sep=",",dec=".")

# upload empty raster
library(rgdal)
# 40 seconds
grid <- readOGR(".", layer = "grid_025dd")

# concatenate data in R
# 2 seconds
mydata$lonlat<-do.call(paste, c(mydata[c("lon", "lat")], sep=""))
grid at data$lonlat<-do.call(paste, c(grid at data[c("LONGITUDE", "LATITUDE")],
sep=""))

# use common variable lonlat to merge data in raster

###### prepare shapefile #####
library(rgdal)        ## Load geographic info
library(maps)         ## Projections
library(maptools)     ## Data management
#library(sm)           ## Data management
library(spdep)        ## Spatial autocorrelation
library(gstat)        ## Geostatistics
library(splancs)      ## Kernel Density
library(spatstat)     ## Geostatistics
library(pgirmess)     ## Spatial autocorrelation
library(RColorBrewer) ## Visualization
library(classInt)     ## Class intervals
library(spgwr)        ## GWR

# Match polygons with data
idx <- match(grid$lonlat, mydata$lonlat)
# Places without information
idxNA <- which(is.na(idx))
# Information to be added to the SpatialPolygons object
dat2add <- mydata[idx, ]
# spCbind uses row names to match polygons with data
# First, extract polygon IDs
IDs <- sapply(grid at polygons, function(x)x at ID)
# and join with the SpatialPolygons
row.names(dat2add) <- IDs
datPols <- spCbind(grid, dat2add)
# Drop those places without information
datPols <- datPols[-idxNA, ]
# write new shapefile
# 7 seconds
writeOGR(datPols, dsn = ".", layer ='sm2000eu28', driver = 'ESRI Shapefile')
# read new shapefile
# 51 seconds
data <- readOGR(".", layer="sm2000eu28")

############################
# intersect nuts with grid #
############################

library(rgdal)
nuts <- readOGR(".", layer = "NUTS_RG_60M_2006")

library(rgeos)
proj4string(data) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")

#
grid <- data
grid at data$lonlat <- NULL
grid at data$lonlat_1 <- NULL
grid at data$ID <- NULL
grid at data$lat <- NULL
grid at data$lon <- NULL
grid at data$ELEVATION <- NULL
grid at data$DAYS_RAIN_ <- NULL

# First find out which grid cells intersect your NUTS polygons
grid_nuts <- gIntersects(grid,nuts,byid = TRUE)

# use the apply() function to calculate the mean, min, and max of your
value.
# The loop makes

for(i in names(grid at data)){
  nuts at data[[paste(i, 'average_value', sep="_")]] <-
apply(grid_nuts,1,function(x) mean(grid at data[[i]][x]))
  nuts at data[[paste(i, 'min_value', sep="_")]] <-
apply(grid_nuts,1,function(x) min(grid at data[[i]][x]))
  nuts at data[[paste(i, 'max_value', sep="_")]] <-
apply(grid_nuts,1,function(x) max(grid at data[[i]][x]))
}

write.table(nuts at data, "nuts_sm2000eu28_unweighted.txt", sep="\t")