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

# 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")
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")

#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)

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

  • 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.
    – Chr
    Nov 20, 2020 at 13:28
  • @Chr thank you for your answer, I tried your suggestion. However, I think I used it once before. My raster files lack an input projection which is NA. Error message: Error in projectRaster projectRaster(r, crs = "+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs ") : input projection is NA I read in another thread before, that if the data lacks a projection, it is useless. Do you know of anything alike?
    – daves_ma
    Nov 20, 2020 at 16:21
  • The data is not useless. Sometimes the projection is mentioned somewhere in the data set's documentation. Otherwise you can try to guess it. The raster package's isLonLat and couldBeLonLat provide 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.
    – Chr
    Nov 20, 2020 at 17:29
  • @Chr well, I did open one NC-file on its own and according to the coordinate reference system it is this: ``` +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 ``` to be honest, I did not generate the NC files, hence I am not awar of its coordinate system. The problem is rather, that I cant relate this projection to any common I would know of. I also tried to project my shapefile in this weird coordinate system. This nevertheless didnt work well.
    – daves_ma
    Nov 20, 2020 at 17:32
  • Then check what CRS the other netCDF files you load have. If they all have the one you posted, you can assign that. If the CRS differs, you should load them separately and project them to a common CRS.
    – Chr
    Nov 20, 2020 at 17:38

2 Answers 2


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")

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

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


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:\\Studium\\Masterarbeit\\05_aris_sBarley_2018\\sBarley',pattern='*.nc',full.names=TRUE)

#generate raster stack

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

# load shapefile for area of interest
field <- read_sf(dsn = 'D:\\Studium\\Masterarbeit\\Studyarea_Polygon\\Polygon for R', layer = "testfield")

#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(fieldnew, add=TRUE)

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

#visual check if raster and shapefile line up 
  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, )

## 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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.