How to extract pixel values of a polygon shape, using R for annual (365 files per year), daily NetCDF raster files?

Question

I am trying to extract a daily value from a study area (given polygon). The files are provided as annual NetCDF files (365 of them per year) for a daily value. However, with GIS programs it would be far too tidious, hence I am looking if R for example can perform an extraction of the needed values as a so to say, time series. I only need the values or the value within a pre-defined polygon (for example with long and lat of it e.g.). Thus, I need a time series, I am not looking into processing one file at a time. Preferably, I want all the Metadata can provide (long, lat, etc) as a Excel file or Txt file.

Is there a way to script it in R, that the files are processed and the respective values are extracted? With my research I already found the package"raster" for R, however I have limited skills considering R. Probably someone can help or give hints how to proceed with my problem.

add update:

install.packages("raster")
install.packages("sf")
install.packages("rgdal")
install.packages("exactextractr")
install.packages("ncdf4")
install.packages("ggplot2")

# load necessary libraries
library(raster)
library(sf)
library(rgdal)
library(exactextractr)
library(ncdf4)
library(ggplot2)

setwd('D:\\Studium\\Masterarbeit\\05_aris_sBarley_2018\\sBarley')

# read NetCDF raster files in given directory
files <- list.files('D:\\Studium\\Masterarbeit\\05_aris_sBarley_2018\\sBarley',pattern='*.nc',full.names=TRUE)
files

#generate raster stack
vname <-"RSS_TOP"
r <- stack(files)
r
newr <- brick(r, CRS=TRUE, values=TRUE)
newr
plot(r)
crs(r)<-("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

#generate raster brick 
crs(newr)<-("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
newr
plot(newr)

# load shapefile for area of interest
field <- read_sf(dsn = 'D:\\Studium\\Masterarbeit\\Ebod und Parrot', layer = "BOKU_Feld")
field
#change projection
fieldnew <- sf::st_transform(field, crs = "+proj=longlat +datum=WGS84 +no_defs") #same projection as raster stack
fieldnew
summary(fieldnew)
plot(fieldnew)

Unfortunately, my projections of the crs of either, the shapefile as well as the rasters dont match. The shapefile is centered in my raster instead of its actual location --> check picture. Is there a way to change the crs of my shapefile to project correctly, or can I read the extent of it in order to project it as raster instead?[![Shapefile in wrong location][1]][1]

#plot a random rasterfile of the stack \\ plot the shapefile on top as preparation to crop and mask
plot(r[[50]])
plot(fieldnew, add=TRUE)

r.mask =mask(r, fieldnew)
r.mask

# mask areas other than area of interest
r.crop <- crop(r, fieldnew)
r.crop

# calculate mean over area of interest
r.mean <- calc(r.crop, mean)
r.mean


  [1]: https://i.sstatic.net/PKXfA.jpg

Answer 1

You would be much better served by thinking through your analysis and keeping it all in R. It is well know that Excel is very limited in statistical analysis and performing task such as aggregated statistical summaries is very simple, and repeatable, in R. The analytical tools that you have on hand in R, for spatial and statistical analysis, are far superior to attempting to perform an analysis in a spreadsheet type software.

First, I would not recommend "extracting" your time-series data to a polygon per se but, rather keeping it in a raster format masked to the polygon. By extracting to a polygon you introduce unnecessary difficulties into your analysis, and potentially lose the "pixel-level" information. You also make your analysis much less tractable.

Here is a quick run-through example of what you are asking. First, add the required libraries

library(raster)
library(exactextractr)
library(sf)

Now we create some example data that is a stack of 50, 100-row by 100 -column rasters with random values assigned. We also create a vector of dates to work with, that theoretically correspond to our stack of rasters. You would use the raster stack or brick function to read in your raster data into a single object representing the time-series.

i=100;j=100;n=50
r <- do.call(stack, replicate(n,raster(matrix(runif(i*j), i, j))))
  ( dates <- seq(as.Date("2000/1/1"), by = "month", length.out = n) )

Here we create a polygon that would represent an area that we may want to subset. I create an sp class object but you can use the sf st_read function to read in something like a shapefile.

poly <- as(extent(0.1045409,0.9528445,0.1265942,0.8459615), "SpatialPolygons") 
  plot(r[[1]])
    plot(poly,add=TRUE)

Here is where you can simply mask your raster time-series to the polygon AOI. This allows one to perform operations on the time-series in very simple ways.

r.mask <- mask(r, poly) plot(r.mask[[1:12]])

For instance, we can take the pixel means across the time-series as such

r.mean <- calc(r.mask, mean)    
  plot(r.mean)

Or, the mean of a subset of dates (ie., only 2002).

r.mean2002 <- calc(r.mask[[grep("2002", dates)]], mean)

These are very simple examples and for each resulting summary raster you could then easily assign the results to your polygon. If you are resolutely on the path of "extracting" the raster values to a polygon and exporting the data, I would still suggest masking the data and then leveraging one of the sp classes "SpatialPixelsDataFrame" because you can get pixel center coordinates as well as easily write out results from the @data slot.

Here we coerce our raster stack to an sp class, assign the date names, and add columns for the [x,y] coordinates. You can then use write.csv to export the results (a data.frame in the @data slot) to a comma-separated flat file. The columns represent each raster (date) and the rows are the associated pixel values.

r.sp <- as(r, "SpatialPixelsDataFrame")
  names(r.sp) <- as.character(dates)
    r.sp@data <- data.frame(coordinates(r.sp), r.sp@data) 
write.csv(r.sp@data, "my_results.csv") 

I do feel compelled to demonstrate a method for "extracting" the data to a polygon but, in your case, it adds an additional processing step. This would be much more relevant if you had several polygons. You can use the exact_extract function, from the exactextractr package, to actually extract the raster values for the polygon. Please note that there is no efficient way to include [x,y] pixel coordinates so, you just get the pixel values.

r.values <- exact_extract(r, as(poly, "sf"))[[1]]
  names(r.values)[-ncol(r.values)] <- as.character(dates)
    str(r.values)