I am trying to aggregate at a 0.1º cell level, using zonal statistics, two variables that come from two different raster files (one in .tif, other in .nc format, links here and here).

Both seem to work fine in the R raster environment as, by loading and plotting them with

r.works <- raster::raster('MOD13Q1_2000-02-18.250m_16_days_NDVI.tif')
r.dwork <- brick("spei01.nc") #this is a multilayer stack

rasterVis::levelplot(r.dwork[[1]]) #plot just the first layer

I obtain clear graphs, respectively attached. However, if I use the ZonalPipe function to do fast zonal statistics within a grid, I obtain a lot of NA values for the latter (not shown/plotted here for neatness), even though both rasters cover 100% of the grid of interest. I obviously set both rasters' crs equal to the grid's one before doing zonal stats. I thus investigated both objects and they seem to differ

> r.works
class       : RasterLayer 
dimensions  : 5000, 10444, 52220000  (nrow, ncol, ncell)
resolution  : 0.01, 0.01  (x, y)
extent      : -26.10815, 78.33185, -10, 40  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
data source : /Users/hannesmueller/Documents/NDVI Data/MOD13Q1_2000-02-18.250m_16_days_NDVI.tif 
names       : MOD13Q1_2000.02.18.250m_16_days_NDVI 
values      : -32768, 32767  (min, max)

> r.dwork[[1]]
class       : RasterLayer 
band        : 1  (of  1380  bands)
dimensions  : 360, 720, 259200  (nrow, ncol, ncell)
resolution  : 0.5, 0.5  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
data source : /Users/hannesmueller/Dropbox/Locus Project Augustin Lavinia Bruno/gis data/spei/spei01.nc 
names       : X1901.01.16 
z-value     : 1901-01-16 
zvar        : spei 

in the band element of the second layer, whose meaning I do not understand. Finally, I exported them with

writeRaster(r.works, "rworks", format = "GTiff")
writeRaster(r.dwork, "rdwork", format = "GTiff")

and the former opens perfectly in a GIS software, whereas the latter does not. That makes me very suspicious about the second raster, which was obtained with the .nc file. Am I importing it into R wrongly and/or missing any step when using zonal statistics/exporting it as geotiffs?

plot of r.work plot of r.dwork

  • When you say you set the crs to match, did you actually reproject one raster to match the other, or just alter the crs directly with e.g. crs(x) <- value?
    – obrl_soil
    Mar 1, 2018 at 0:16
  • 1
    I did set as you mentioned, using crs(x) <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"), which is the CRS of the polygon I am trying to aggregate the data in. But beyond the aggregation, it seems very weird to me that exporting the raster obtained from the .nc stack does not provide a 'readable' raster in a GIS software. Any hints on that? Mar 1, 2018 at 9:02

1 Answer 1


Per your comment above, you have not actually reprojected your rasters. crs() is essentially a metadata editing and reporting function, but its documentation is brief so many people don't realise at first. To alter the actual data, you need to use projectRaster() instead. Your GIS is most likely confused about a mismatch between the claimed CRS and the actual data.

For large rasters, projectRaster() can be slow. You might consider using GDALwarp instead. If you have GDAL installed and on your system path, you can call it in R using system2(), e.g.


# just some random data i have in EPSG:3577
test <- raster('C:/DATA/rmap2.tif') 

# gdalwarp docs - http://www.gdal.org/gdalwarp.html
test_4326 <- system2('gdalwarp',
                     args = c(
                       # source projection (proj4 string, mind the quotes)
                       '-s_srs', paste0('"', crs(test), '"'),
                       # desired projection (EPSG code)
                       '-t_srs', 'EPSG:4326',
                       # resampling method
                        '-r', 'bilinear',
                       # src file - can't just pass in R object, must refer to on disk file
                       # dest file - must also be on disk
                       file.path('C:', 'data', 'test_4326.tif')))

test_4326 <- raster(file.path('C:', 'data', 'test_4326.tif'))
  • Reprojecting the raster as you mentioned (either with projectRaster() or with GDALwarp does not seem to solve my issue of problems with the final exported tif file. I noticed that, if I substitute the many NAvalues in the raster with 0es, I manage to export it well. Nevertheless, if running zonal stats, I get the same problem of 98% of my cells with NA. Mar 1, 2018 at 10:54
  • hmm ok. The only other problem I can see is that the cells in your second grid are 0.5 x 0.5 degree but you're trying to 'aggregate' to a smaller size, 0.1 x 0.1 degree. The zonal stats function is probably confused as a result - things would be going pear-shaped at the use of crop(). That zonal stats script is...not exactly robust, although its certainly well-intended. You should raster::disaggregate() r.dwork.
    – obrl_soil
    Mar 1, 2018 at 11:16
  • That is actually my guess here. I did aggregation at coarser spatial levels and it seems to work. Mar 1, 2018 at 13:34

Your Answer

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

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