Extract pixel values of dataset (NetCDF) based on shapefile in R

Question

I want to extract certain pixel values from a raster stack/brick from one layer (there are two Layers stored in the data - RSS_TOP and RSS_SUB --> if possible I would take both, or separately as long as it works).
Nevertheless, currently I found a solution to have the same projection for my shapefile as well as the raster stack/brick. My problem is, that if I check a plotted rasterfile and shapefile, they dont line up where they should. Thats why there is an error if I want to mask and crop and extract finally my mean pixel values.
I am not sure what causes the problem, since all my data should have the same projection.
# load necessary libraries
library(raster)
library(sf)
library(exactextractr)
library(ncdf4)

# read NetCDF raster files in given directory
files <- list.files(pattern='*.nc', full.names=TRUE)

#generate raster stack using the RSS_TOP layer
r <- stack(files, varname="RSS_TOP")
r
crs(r) <- "+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs"

# load shapefile for area of interest
field <- read_sf(dsn = 'Polygon for R', layer = "testfield")
field

#change projection
fieldnew <- sf::st_transform(field, crs = "+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs ") #same projection as raster stack

####
#betterr <- projectRaster(r, crs=field)

### here starts my question
r.mask = mask(newr, fieldnew)
r.mask

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

Robert Hijmans · Answer

You are clearly making a mistake by setting the crs of the raster data. Apparently it is +proj=lcc +lat_1 ... (and nothing weird about it) but you set it to +proj=utm +zone=33. Your problem might go away if you just leave that out. But you also make a mistake by projecting the raster data (that is not a good approach if your goal is to extract values).
I edited your question to remove a lot of clutter (like plot statements that are useless for us if you do not include the output), but can you please add to your question the information we need to understand what is going on? That is can you show what R shows you when you do:
r <- stack(files, varname="RSS_TOP")
r

#field <- read_sf(dsn = 'Polygon for R', layer = "testfield")
#field
# read_sf works, but for easier comparison with the RasterStack above use
field <- shapefile("Polygon for R.shp")
field

You do not actually show code that does an extract, whereas you say that is where you have a problem.

daves_ma · Answer

Thank you for your answer @Rober Hijmans.
I got a solution for my problem, I guess. I projected my shapefile in the raster's crs. Which resulted in an overlay of my raster and shapefile as it should be. Your shortening of the script helped as well. In the end I face a little problem, that I cant double check if the data output lines up with my NC-files.
# read NetCDF raster files in given directory
files <- list.files('D:StudiumMasterarbeit5_aris_sBarley_2018sBarley',pattern='*.nc',full.names=TRUE)
files

#generate raster stack

r <- stack(files, varname="RSS_TOP")
r

# load shapefile for area of interest
field <- read_sf(dsn = 'D:StudiumMasterarbeitStudyarea_PolygonPolygon for R', layer = "testfield")
field

#change projection
fieldnew <- sf::st_transform(field, crs = "+proj=lcc +lat_1=46 +lat_2=49 +lat_0=47.5 +lon_0=13.333333 +x_0=0.0 +y_0=0.0 +a=6377397.16 +b=6356078.96 +units=m -te -380501.324 -180498.983 320498.676 220501.017") #same projection as raster stack

#plot raster and shapefile for visual assessment
plot(r[[7]])
plot(fieldnew, add=TRUE)

plot(r)
  plot(fieldnew, add=TRUE)
Warning:
In plot.sf(fieldnew, add = TRUE) : ignoring all but the first attribute

#visual check if raster and shapefile line up 
plot(r[[105]])
  plot(fieldnew, add=TRUE)

## extract by weights of polygon Raster,SpatialPoylgon and relevant observations
output <- extract(r [[101:196]], fieldnew, weights=TRUE, normalizeWeights=FALSE, layer=2, method="bilinear", df=TRUE, )
output

## write csv file with extracted data
write.csv(output, "Barley_Vegetationperiod.csv")

My last step would be to extract my relevant data [101:196]. I want to extract the mean of my pixels within the shapefile. There should be 2 layers with relevant data "RSS_TOP and RSS_SUB", both I would like to extract.
The code "write.csv" could be little more sophisticated with 3 separate columns featuring the data output.

Extract pixel values of dataset (NetCDF) based on shapefile in R

2 Answers

Add your own answers!

Ask a Question