¿Convertir ESRI Shapefile Raster utilizando R - mayor área de polígono a valor de la celda raster?

Question

Quiero crear un ráster de un shapefile de ESRI. Varios polígonos puede existir en una celda ráster.

Si es así, quiero asignar el atributo de la variable que ocupa la mayor área de la celda ráster. Si, por ejemplo, Polygon1 ocupa el 10%, Polygon14 25% y Polygon5 65% de la trama de la célula, la célula debe contener el atributo variable de Polygon5.

El uso de la R-pacakge rasterizar{trama}, los valores son transferidos en caso de que el polígono cubre el centro de una celda ráster. Por supuesto, las funciones de, como mínimo, máximo, media, o el carácter de valores como la 'primera', 'apellido', 'suma', 'contar' están disponibles. Por defecto es 'último'.

Puede descargar el ejemplo shapefile aquí.

Hasta ahora mi código es este:

setwd("your_wd")

require(sp)
require(raster)
require(rgdal)

sp_df <- readOGR(getwd(), "stok_subs") # reads the shapefile into a large polygon
                                       # data frame

extent(sp_df) # shows the extent of sp_df

summary(sp_df) # simple summary

# Object of class SpatialPolygonsDataFrame
# Coordinates:
#        min      max
# x 32539350 32544695
# y  5732607  5735439
# Is projected: TRUE 
# proj4string :
# [+proj=tmerc +lat_0=0 +lon_0=9 +k=0.9996 +x_0=32500000 +y_0=0 +ellps=GRS80+units=m+no_defs]
# Data attributes:
# X_CENTR           Y_CENTR            ATTRIBUTE 
# Min.   :3539714   Min.   :5734725   9.4-.2.3 :21  
# 1st Qu.:3540436   1st Qu.:5735323   23.4-.2.3:19
# Median :3541226   Median :5735830   9.4.2.3  :19  
# Mean   :3541263   Mean   :5735824   23.4.3.5 :18  
# 3rd Qu.:3542031   3rd Qu.:5736358   23.4.2.3 :16  
# Max.   :3543461   Max.   :5737024   19.4.3.5 :12  
#                                     (Other)  :89  

r <- raster(extent(sp_df)) # creates a raster with the extent of sp_df
projection(r) <- proj4string(sp_df) # uses the projection of the shapefile
                                    # for the raster    
res(r)=10.0 # sets the resolution of the raster to 10 m

r10 <- rasterize(sp_df, field="ATTRIBUTE", r) # converts sp_df to a raster
                                              # assigning the last attribute
                                              # variable in a cell (default)

writeRaster(r10, filename="stok_subs.tif", format="GTiff") # writes the raster
                                                           # to a file

plot(r10) # plots raster

Ninguna de las funciones parecen cubrir mis necesidades. Sin embargo, prefiero operar con R, como las operaciones con cantidades muy grandes de datos raster (por ejemplo, toda la Baja Sajonia @ una resolución de 10 m) pueden ser más fáciles de manejar que con soluciones SIG.

Alguien por ahí que tiene una (o ninguna) de solución para mi problema?

Answer 1

No estoy seguro de donde se obtiene la mencionada 'atributo' variable de, así que supuse rownames(sp_df@data) (que van de 0 a 193) como la salida deseada. De todos modos, aquí está mi enfoque que probablemente no es el más rápido (o el más conveniente en general), pero funciona.

Básicamente, el guión comienza a recorrer todas las celdas de la trama de la plantilla (a 100 m de resolución espacial) y a los cultivos sp_df por la magnitud de la actual píxel (de los que se extrae mediante la configuración de todos, pero la celda actual a NA y, a continuación, convertir a los polígonos). Después de eso, comprueba si no, uno o más de un polígono IDs se encuentran dentro de la corriente de baldosa, calcula las áreas de los Identificadores de referencia y devuelve el ID que cubre el área más grande. Posteriormente, los valores extraídos se insertan en la referencia celdas ráster.

# Required libraries
library(raster)
library(rgdal)
library(rgeos)

# Import shapefile
sp_df <- readOGR("data/", "stok_subs") 

# Raster template (100 m spatial resolution)
r <- raster(extent(sp_df))
projection(r) <- proj4string(sp_df)
res(r) <- 100

# Per pixel, identify ID covering largest area
r_val <- sapply(1:ncell(r), function(i) {

  r_dupl <- r
  r_dupl[i] <- 1
  p <- rasterToPolygons(r_dupl) # Current cell -> polygon
  sp_df_crp <- crop(sp_df, p)   # Crop initial polygons by current cell extent

  # Case 1: no polygon intersecting current cell
  if (is.null(sp_df_crp)) {
    return(NA)
  # Case 2: one polygon intersecting current cell  
  } else if (nrow(sp_df_crp@data) < 2)  {
    return(rownames(sp_df_crp@data)) 
  # Case 3: multiple polygons intersecting current cell
  } else {
    areas <- gArea(sp_df_crp, byid = TRUE)
    index <- which.max(areas)
    return(rownames(sp_df_crp@data)[index])
  }
})

# Write ID values covering the largest area per pixel into raster template
r[] <- as.numeric(r_val)
plot(r)
plot(sp_df, border = "grey45", add = TRUE)