isciences / exactextractr

R package for fast and accurate raster zonal statistics
https://isciences.gitlab.io/exactextractr/
274 stars 26 forks source link

Extraction of weighted mean from raster generating “row <= x@nrows is not TRUE” error #12

Closed NB-Bio closed 4 years ago

NB-Bio commented 4 years ago

Hi, I am using the exactextractr R package to calculate the mean value of raster cells overlapped by a set of multi-polygon geometries (each multi-polygon represents the distribution of a species).

Here, the example files.

Since the raster is in geographic coordinates and the cells don't have all the same area, I actually need to calculate the weighted mean of the raster values, using the coverage fraction and the cell area as weights. So:

I read the raster and the shapefile:

Read the raster

r <- raster("raster.grd") r class : RasterLayer dimensions : 278, 715, 198770 (nrow, ncol, ncell) resolution : 0.5, 0.5 (x, y) extent : -179, 178.5, -55.5, 83.5 (xmin, xmax, ymin, ymax) crs : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 source : C:/Users/MS/Analysis/raster.grd names : w001001 values : 0, 81.45476 (min, max)

Read the shapefile

all <- st_read(dsn = 'all.shp') Reading layer all' from data sourceC:\Users\MS\Analysis\all.shp' using driver `ESRI Shapefile' Simple feature collection with 5 features and 1 field geometry type: MULTIPOLYGON dimension: XY bbox: xmin: -180 ymin: -23.60646 xmax: 180 ymax: -4.725997 epsg (SRID): 4326 proj4string: +proj=longlat +datum=WGS84 +no_defs

Then, following the package instructions, I create a two-layer RasterStack, with the first layer containing my raster and the second layer containing the area of each cell:

Create stack

stk <- stack(list(r= r, a= area(r)))

Finally, I used the exact_extract function from the exactextractr package, with a customized function to calculate the weighted mean:

all$wm <- exact_extract(stk, all, function(values, coverage_frac) weighted.mean(values$r, values$a*coverage_frac, na.rm = T)) | | 0%

but I obtained the following error:

Error in .local(x, ...) : row <= x@nrows is not TRUE

The error is due to a single geometry(field "scntfcN" = Mirimiri acrodonta, here the shapefile).

It's a tiny multipolygon located on the Fiji Island (a theoretically single polygon split into two by the anti-meridian line) not intersecting the original raster.

Strangely, if I try to calculate the (no-weighted) mean through the predefined function of the package, it works well:

all$m <- exact_extract(r, all, 'mean') returning NaN for the problematic geometry (since it doesn't intersect the raster)

What could be the cause of the error?

Thanks, NB

dbaston commented 4 years ago

Thank you for the excellent bug report. I've reproduced it locally and hope to have a fix available shortly.

NB-Bio commented 4 years ago

Hi, I found the problem. The problematic multi-polygon has its coordinates exceeding the extent of the input raster, so it's not a a package issue!

Thank you and congrats for this very useful package! NB

dbaston commented 4 years ago

The problematic multi-polygon has its coordinates exceeding the extent of the input raster

The package is supposed to automatically extend the input raster with NA values in that case, though. As you pointed out, this works in the unweighted calculation.

dbaston commented 4 years ago

I have a fix, if you'd like to test it:

remotes::install_git('https://gitlab.com/isciences/exactextractr', ref='gh-12')

NB-Bio commented 4 years ago

Now it works perfectly, thank you so much!