La superposición de cuadrícula vs. polígono en R da resultados diferentes que en QGIS

Question

Tengo un shapefile de polígonos de rango. Cada conjunto de polígonos coloreados representa un grupo diferente de especies/géneros/lo que sea.

A continuación superpongo una cuadrícula de 5x5 a este mapa. Y luego uso atributos de unión por ubicación para encontrar el número de cuadrículas que intersecan cada polígono

El resultado final me da una columna en la tabla de atributos COUNT que indica el número de cuadrículas intersecadas por los polígonos.

EDIT: En otras palabras, me da la ocupación del sitio de cada especie.

################## FASE SIGUIENTE

Quiero duplicar este mismo proceso completamente dentro de R. Primero cargo el mismo shapefile de polígonos de rango en R usando el paquete rgdal.

RangeShape<-readOGR(".","overflowExample07162015")

A continuación, genero una cuadrícula análoga de 5x5 en R utilizando la siguiente función y el paquete raster.

# Generate lat/lng grid for use with mammals dataset
generateGrid<-function(XMin=-180,XMax=180,YMin=-90,YMax=90,BinSize=5) {
    LatitudeRows<-length(seq(YMin,YMax,BinSize))
    LongitudeColumns<-length(seq(XMin,XMax,BinSize))
    Grid<-raster(extent(matrix(c(XMin,YMin,XMax,YMax),nrow=2)),nrow=LatitudeRows,ncol=LongitudeColumns,crs="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
    Grid[]<-1:ncell(Grid)
    Grid<-as(Grid,"SpatialPixelsDataFrame")
    return(Grid)
    }

# Create a GridShape object
GridShape<-generateGrid(XMin=-180,XMax=180,YMin=-90,YMax=90,BinSize=5)

A continuación intento superponerlos utilizando la función over() del paquete sp con la siguiente función.

# Find Occupancy of taxa using shapefile
shapeOccupancy<-function(RangeShape,GridShape) {
    Occupancies<-vector("numeric",length=nrow(RangeShape))
    for (i in 1:nrow(RangeShape)) {
        Intersect<-over(GridShape,RangeShape[i,])
        Present<-apply(is.na(Intersect),1,all) # identify grids that do not intersect range polygons
        Occupancies[i]<-nrow(Intersect[Present!=TRUE,])
        }
    return(Occupancies)
    }

 Occupancy<-shapeOccupancy(RangeShape,GridShape)

################## El problema

Las cuentas dadas por mi proceso R no son ni siquiera cerca de las cuentas dadas por el proceso QGIS. Podría entender ligeras diferencias, pero estas son bastante grandes.

# Sample of differences
                     RResults  QGISResult
Abrothrix                7         23     
Alouatta                31         69     
Aplodontia               1         6     
Artibeus                47         88     
Ateles                  24         48     
Avahi                    1         7     
Balantiopteryx           3         17

Mi mejor suposición actualmente es que esto es un resultado del hecho de que algunos de los polígonos de rango son discontinuas y R tiene un tiempo difícil manejar eso, pero no estoy seguro. Basándome en algunas comprobaciones manuales (es decir, contando manualmente las intersecciones de cuadrículas y polígonos en el mapa), los resultados de QGIS parecen ser precisos y son los resultados de R los que están equivocados.

EDIT: Un subconjunto de los polígonos de rango se incluye en el siguiente repositorio github. https://github.com/aazaff/StackOverflowExamples/

Answer 1

Si he entendido bien el resultado que desea, le gustaría disponer de un recuento (riqueza) de especies para cada celda de la cuadrícula en el raster definido. No puedo hablar de las diferencias entre R y QGIS, pero se me ocurrió una forma mucho más optimizada y rápida de llevar a cabo su análisis. Aprovecho el paquete raster y utilizo una pila raster para acumular especies. El flujo de trabajo es el siguiente.

Cargar las bibliotecas necesarias y leer los polígonos de distribución de las especies.

library(raster)
library(sp)
library(rgdal)

setwd("D:/TMP/spp")
spp <- readOGR(getwd(), "overflowSubset07162015")

Crear trama con la resolución o filas/columnas deseadas. Este es el raster de referencia que se utilizará para rasterizar los polígonos de las especies.

xmin=-180; xmax=180; ymin=-90; ymax=90; bin=5
bins <- raster(extent(matrix(c(xmin,ymin,xmax,ymax),nrow=2)),
                 nrow=length(seq(ymin, ymax, bin)), 
                 ncol=length(seq(xmin, xmax, bin)),
                 crs= "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
bins[] <- 1:ncell(bins)

Ahora podemos crear rásters binarios para cada especie correspondientes al ráster de referencia. Primero creamos una pila vacía y luego la rellenamos subdividiendo cada especie y rasterizándola con el raster de referencia. Para mayor claridad, asignamos los nombres de cada especie a cada ráster de la pila.

spp.pa <- stack()
  for(i in unique(spp@data$genus_name) ) {
    s <- as(spp[spp$genus_name == i,], "SpatialPolygons")
    spp.pa <- addLayer(spp.pa, rasterize(s, bins, fun = 'count', background = 0))
   } 
names(spp.pa) <- unique(spp@data$genus_name)

Ahora podemos calcular la suma (recuentos de especies en cada píxel) utilizando la función calc en la pila raster.

spp.sum <- calc(spp.pa, fun = sum) 
plot(spp.sum)  

Si desea obtener la cobertura fraccional de una especie en un píxel determinado, puede utilizar el argumento getCover en la función rasterize. De este modo, puede crear un argumento condicional en el bucle for para que, si una especie tiene una cobertura inferior a un porcentaje definido, el píxel sea 0. A continuación se muestra un ejemplo rápido que establece los píxeles con una cobertura < 30% en 0 y, si no, en 1.

p = 30
spp.pct <- stack()
  for(i in unique(spp@data$genus_name) ) {
    s <- as(spp[spp$genus_name == i,], "SpatialPolygons")
    sfp <- rasterize(s, bins, getCover = TRUE, background = 0)
    sfp[sfp < p] <- 0
    sfp[sfp >= p] <- 1
    spp.pct <- addLayer(spp.pct, sfp)
   } 
names(spp.pct) <- unique(spp@data$genus_name)

spp.pct.sum <- calc(spp.pct, fun = sum)

par(mfrow=c(2,1))
  plot(spp.sum)
  plot(spp.pct.sum)