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.
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]]
#plot a random rasterfile of the stack \\ plot the shapefile on top as preparation to crop and mask plot(r[]) 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 : https://i.stack.imgur.com/PKXfA.jpg