Crear un raster con los valores de las celdas de la cuadrícula basados en el área cubierta por los polígonos

Question

Tengo datos de un polígono con el área cubierta por los bosques (los datos están aquí - https://www.dropbox.com/s/zgckliydalljw6a/sp_data.zip?dl=0 ). Quiero convertir los polígonos en raster. El valor de cada celda de la cuadrícula debe basarse en el área cubierta por el polígono. Por ejemplo, si el tamaño de la celda de la cuadrícula es de 100m x 100m (10000m2) y los polígonos cubren 5000m2 de esta celda - el valor debe ser 50 (o 0,5). A continuación, puede encontrar mi solución propuesta:

Datos en bruto:

library('rgdal')
library('raster')
library('rgeos')
library('maptools')

sp_data <- readOGR(".", "sp_data")
spplot(sp_data)

He creado un nuevo raster con una resolución de 100 metros y una extensión igual a la del cuadro delimitador del polígono:

r <- raster()
bb <- extent(290000, 300000, 500000, 510000)
extent(r) <- bb
res(r) <- 100

También he creado un polígono de fondo - con los valores de las áreas sin bosque iguales a 0. Después, he fusionado estos dos datos poligonales:

sp_back <- as(extent(r), "SpatialPolygons")
proj4string(sp_back) <- proj4string(sp_data)
sp_back <- gDifference(sp_back, sp_data)
plot(sp_back, col='red')
plot(sp_data, col='green')
sp_back <- SpatialPolygonsDataFrame(sp_back, data=data.frame(value=rep(0, length(sp_back)), row.names=row.names(sp_back)))
sp_back <- spChFIDs(sp_back, "new_id")
sp_bind <- spRbind(sp_data, sp_back)
spplot(sp_bind)

Además, he creado una nueva trama con una resolución reducida (500 metros):

r2 <- r
res(r2) <- 500

He rasterizado los datos de los polígonos en el primer raster:

sp_raster <- rasterize(sp_bind, r, field="value", fun=max)
spplot(sp_raster, aspect='iso')

Al final, he remuestreado los valores del primer raster en el segundo raster:

sp_raster2 <- resample(sp_raster, r2, method='bilinear')
spplot(sp_raster2, aspect='iso')

Basándome en los resultados, tengo algunas preguntas:

1. En primer lugar, ¿mi solución y el resultado son correctos?
2. La conversión entre vector y raster simplifica de alguna manera mis datos (por ejemplo, la pérdida de polígonos pequeños) y probablemente añade algunos errores. El resultado del remuestreo se ve afectado por esta simplificación. ¿Existe alguna solución alternativa para calcular la parte del área de los polígonos en la celda de la cuadrícula?

Las respuestas con R, Grass o SAGA serán bienvenidas.

Answer 1

Ya que has puesto GRASS en las etiquetas, aquí tienes una solución basada en GRASS:

En primer lugar, hay que saber exactamente en qué sistema de coordenadas están los datos originales (como siempre con GRASS). Veo que el archivo *.prj contiene "TRANSVERSE MERCATOR" pero no es una de las zonas UTM estándar. Ya que no has mencionado el código EPSG de estos datos, aquí está la cadena proj4 para crear un nuevo GRASS Location que coincida:

+proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=WGS84 +units=m +no_defs

Ahora inicie GRASS, cree un LOCATION (basado en lo anterior) y un MAPSET.

Y ahora importa el shapefile a GRASS ejecutando

v.in.ogr input="sp_data.shp" output=sp_data

Añade una columna para el área de cada característica, y calcula su área:

~~v.db.addcolumn map=sp_data columns="area_sqm INT" v.to.db sp_data option=area column=area_sqm unit=meters~~

Ahora convierta a raster, utilizando la columna area_sqm como valor raster. Pero primero establecer la región de cálculo de GRASS :

g.region -p vector=sp_data res=100

Ahora cree una malla vectorial de polígonos a la resolución de la trama:

v.mkgrid map=grid

Interseca con los polígonos del bosque, y calcula las áreas de estos polígonos de intersección: Interseca con los polígonos del bosque, disuelve los polígonos dentro de la misma celda de la cuadrícula, calcula las áreas de estos polígonos de intersección y traslada el valor del área a las celdas de la cuadrícula:

v.overlay ain=grid bin=sp_data operator=and output=forest_grid_intersect
v.dissolve --overwrite input=forest_grid_intersect@jn column=a_cat output=forest_grid_intersect_dissolve
v.db.addtable map=forest_grid_intersect_dissolve columns="area_sqm INT"
v.to.db forest_grid_intersect_dissolve option=area column=area_sqm unit=meters 
v.db.join map=grid@jn column=cat other_table=forest_grid_intersect_dissolve other_column=cat

v.db.addcolumn map=forest_grid_intersect columns="area_sqm INT" v.to.db forest_grid_intersect option=área column=área_sqm unit=metros

y convertir a raster, utilizando el atributo de área como valor de raster:

v.to.rast input=grid type=area output=sp_data_rast use=attr attribute_column=area_sqm

Y ya está. Espero haberlo hecho bien esta vez :-) .

Answer 2

Estaba tratando de averiguar cómo hacer lo mismo en r. Este es el enfoque que se me ocurrió:

-Crear una trama que tenga las dimensiones y el tamaño de celda que desee

-Utilizar la función 'rasterize' del paquete raster, utilizando la opción "getCover"

library(raster)
wetlands <- readShapePoly("data/wetlands.shp", proj4string = CRS("+init=EPSG:26907))

r <- raster(ncols = 200, nrows = 200, resolution = 5, crs = CRS("+init=EPSG:26907))

wetlandsAREA <- rasterize(wetlands, r, getCover = TRUE)

Creo que esto hace esencialmente el mismo proceso que en el post original. La función rasterizar divide la célula raster en 100 subunidades, y la salida es la proporción de los que fueron cubiertos por el polígono.