1

Are there any known bugs with the getCover argument in the rasterize function? (R package 'raster'). or has anyone got a fix for the following:

I have a series of landcover polygons, here is a small sample of them in shp format, that i want to rasterize based on coverage at 25m resolution. There are 5 types of landcover.

I would like to perform a rasterization that shows the percentage of each cell that is covered, per Land Cover type (i get 5 resulting rasters that I can stack), using the getCover option in 'rasterize' in R.

Here is some code, im very sorry it is only a download dataset;

require(raster)
require(rgdal)

setwd("C:\\your_wd")

# read in shapefiles, extract .prj
fields <- readOGR(".", "sample")
projGB <- projection(sample)

# make a 25m raster (made specifically for above data)
rOS <- raster(nrows=16, ncols=27, xmn=254110.156238, xmx=254785.156238, ymn=251530.203086, ymx=251930.203086, crs=projGB)

# rasterize the polygons per LandCover code in a loop, using getCover

cats <- unique(fields$LAND_CAT)
st <- stack()

for (i in cats){
 x <- fields[fields$LAND_CAT == i,]
 y <- rasterize(x, rOS,getCover=TRUE)
 writeRaster(y, filename=paste("fields_",i,"_25m"),  format="GTiff", overwrite=TRUE)
 st <- stack(st, y)
}

Unfortunately, the produced GeoTiffs do not show all coverage stats across the polygons. Only some polygons appear to have been used, with no detectable pattern. The biggest category 'GP' is the most problematic.

I have:

  1. plotted all polygons as a group, per Land Cover and individually to check they all plot (they do)
  2. used fortify to check for holes etc, there are none

    require(ggplot2) fields$group <- 1 fields.f <- fortify(fields,region="group")

  3. rasterized the whole dataset disregarding landcover type and not using getCover, just to check - all is fine.

  4. plotted individual polygons and rasterized each one using getCover - no problem
  5. used 2 other small samples from the big dataset but the problem persists
  6. used another shapefile from another project, problem persists

I wonder why getCover in rasterize does not work and if anyone has had this problem?

R version 3.0.3


If and when I get this fixed, I will use the raster stack to then determine which land category has the biggest value in every cell to make a raster that is based on majority cover instead of centroid. Doing this in R seems long winded to me.

  • furthermore, ive tried rasterizing to .grd, .asc and .tif with the same errors. I've also digitized up from scratch a dummy dataset and the same problems occur. The errors are inconsistent as well, some of the polygons are used in the getCover function and some arent, and they seem to change each time i run it! – Sam Nov 6 '15 at 9:52
  • 1
    The link to the example files does not work – Robert Hijmans Nov 9 '15 at 19:00
  • Thanks RobertH, I have no idea why the link was broken – Sam Nov 10 '15 at 11:15
2

This is a bug in the current version or raster. It has been fixed, I think, in version 2.4-25 (will be available on R-Forge in a couple of hours).

  • This is great, thanks a lot RobertH. had a bit of a fight with version of R etc but all up and running now and getCover seems to work :) – Sam Nov 12 '15 at 9:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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