Muestreo ráster con vectorial tipo polígono en R

El lenguaje estadístico R es muy potente y versátil y permite la instalación, a través de repositorios, de una serie de librerías (módulos) para manipular formatos vectorial y ráster provenientes de Sistemas de Información Geográfica (SIG).

Se basa en álgebra vectorial y, aunque admite programación en el sentido clásico, algunos autores sugieren que se agoten otras opciones antes de intentar incluir loops que ralentizan y hacen ineficiente el lenguaje.

Como es software libre no existen intenciones de enmascarar su potencia en versiones “más económicas”, como se suele hacer con el software privativo, y los métodos prácticamente admiten y soportan una cantidad ingente de parámetros que permiten realizar cualquier tarea y nuestro deber es leer con detenimiento la documentación para que ello sea posible.

El ejemplo que se va a contemplar a continuación consiste en seleccionar los values de un ráster (aleatorio.tif) que está cubierto por un poligono (single_polygon.shp) y posteriormente con las coordenadas de cada píxel obtener el correspondiente objeto de puntos. Realizar esto con Python requiere cierto esfuerzo de programación con las librerías disponibles (GDAL/OGR principalmente) pero en R está prácticamente todo hecho.

Los comandos para cargar y visualizar los objetos se indican a continuación:

setwd('c:/Users/Usuario/pyqgis_data')
library(raster)
r <- raster('aleatorio.tif')
s_pol <- shapefile('single_polygon.shp')
plot(r)
plot(s_pol, add=T)

y este es el resultado:

R2

Para el muestreo ráster se añaden cellnumber y df (ambos a TRUE) como parámetros al método extract de la librería raster. Los values se encuentran en la columna “cell” del Data Frame.

require(rgeos)
df <- extract(r, s_pol, cellnumber=T, df=T)
names(df)
  [1] "ID"        "cell"      "aleatorio"
values <- df$cell
values
  [1] 128 129 130 131 142 143 146 147 148 149 150 151 152 162 163 164 165 166
  [19] 167 168 169 170 171 172 173 175 176 177 182 183 184 185 186 187 188 189
  [37] 190 191 192 193 194 195 196 197 202 203 204 205 206 207 208 209 210 211
  [55] 212 213 214 215 216 217 218 225 226 227 228 229 230 231 232 233 234 235
  [73] 236 237 238 247 248 249 250 251 252 253 254 255 256 257 258 270 271 272
  [91] 273 274 275 276 277 291 292 293 294 295 296 312 313 314

Ahora sólo resta extraer las coordenadas del ráster y las del polígono serán las correspondientes a los values; tal como se indica a continuación:

coor <- coordinates(r)
new_coor <- coor[values,]

Los objetos de puntos espaciales en R pertenecen a la clase SpatialPointsDataFrame. La secuencia de comandos para construirlo (incluyendo proyección) es la siguiente:

data <- data.frame(id = 1:length(values))  #añade un id
pts <- SpatialPointsDataFrame(new_coor, data, match.ID = TRUE)
projection(pts) <- projection('+init=EPSG:32612')

Finalmente, para añadir el objeto pts como imagen al conjunto anteriormente visualizado se tiene:

plot(pts, add=T)

obteniéndose lo siguiente:

R3

Por otra parte, si lo quieren disponible como shapefile en el directorio de trabajo:

shapefile(pts, 'polygon_points.shp').

Fácil y rápido.

Esta entrada fue publicada en Lenguaje R. Guarda el enlace permanente.

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión /  Cambiar )

Google photo

Estás comentando usando tu cuenta de Google. Cerrar sesión /  Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión /  Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión /  Cambiar )

Conectando a %s