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
crs(r)<-("+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs ")does not project the raster object to another CRS. It simply overwrites the raster object's CRS attribute without changing the grid.
r <- projectRaster(r, crs = "+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs")projects the raster object to the mentioned CRS.
projectRaster(r, crs = "+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs ") : input projection is NAI read in another thread before, that if the data lacks a projection, it is useless. Do you know of anything alike?
couldBeLonLatprovide an idea of whether the data might be in lonlat format. According to my experience, data sets with a missing CRS attribute tend to be in lonlat format, in particular
+proj=longlat +datum=WGS84 +no_defs. Simply assigning any CRS you like is not a correct solution.