Extraer ráster de ráster utilizando el archivo de forma Polígono en R

13

Soy nuevo en R y estoy usando el paquete raster. Tengo un problema para extraer polígonos de un archivo ráster existente. Si yo uso

extract(raster, poly_shape)

La función en el ráster siempre crea una lista con los datos. Lo que realmente quiero es extraer otro archivo ráster que pueda cargar con ArcGIS nuevamente. Después de leer un poco más, creo que la función de recorte es lo que realmente necesito. Pero cuando trato de usar esta función

crop(raster, poly_shape)

Me sale este error:

Error in .local(x, y, ...) : extents do not overlap
In addition: Warning message:
In intersect(extent(x), extent(y)) : Objects do not overlap

Los archivos ráster y poly_shape son los mismos para ambas funciones. ¿Puedes decirme qué podría estar mal aquí? ¿Es correcto que la función de recorte cree otro ráster y no una lista?

EDITAR : La función de extensión () no funciona para mí. Todavía recibo el mismo error. ¡Pero estoy seguro de que los 2 conjuntos de datos se superponen! Con el

extract(raster, poly_shape)

Obtengo los datos correctos. Solo como una lista y no como una trama como quiero tenerla. Acabo de cargar los conjuntos de datos en ArcGIS antes y se ajustan muy bien, así que no verifiqué la proyección. Ahora lo intenté

projection(raster) # "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m +no_defs"
projection(poly_shape) # "+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs"

y puedes ver que las proyecciones no encajan. La función de extracción parece ser capaz de transformar automáticamente los archivos de la manera correcta. Lo sé porque hice lo siguiente:

  • Corté la parte exacta del polígono que extraje en R también en ArcGIS
  • Calculé la suma de todos los valores del polígono R extraído (lista)
  • Calculé la suma de todas las celdas ráster que corté en ArcGIS

Los 2 tienen exactamente el mismo resultado, así que supongo que la conclusión debería ser que la función de extracción funcionó correctamente. Ahora tengo 2 opciones, supongo:

  1. Necesito una forma de volver a sacar un Ráster de la lista extraída o
  2. Los 2 conjuntos de datos (raster + poly_shape) deben usar la misma proyección y la función de recorte debería funcionar

¿Qué sugerirías hacer aquí?

Lars
fuente
¿Qué pasa si es un rgbi ráster de 4 bandas? Las bandas están perdidas hasta ahora ...
Doris
¡Bienvenido a GIS SE! Como nuevo usuario, asegúrese de hacer un breve recorrido . Luego considere editar su respuesta para proporcionar información adicional y referencias. Vea ¿Cómo escribo una buena respuesta? para más información.
Andy
Si tiene una nueva pregunta, hágala haciendo clic en el botón Hacer pregunta . Incluya un enlace a esta pregunta si ayuda a proporcionar contexto. - De la opinión
Erik

Respuestas:

19

La función de extracción se comporta exactamente como debería. Puede forzar la función de recorte para usar la extensión del polígono y luego enmascarar el objeto para devolver el ráster exacto que representa el área del polígono. Si continúa recibiendo el error, significa que sus datos, de hecho, no se superponen.

Tenga en cuenta que R no realiza la proyección "sobre la marcha", por lo tanto, verifique sus proyecciones. Puede verificar si sus extensiones se superponen usando la función "extensión ()".

Aquí hay un ejemplo de recorte usando la extensión del polígono y luego enmascarando el ráster resultante usando el polígono "rasterizado".

# Add required packages
require(raster)
require(rgeos)
require(sp)

# Create some data using meuse 
data(meuse)
  coordinates(meuse) <- ~x+y
    proj4string(meuse) <- CRS("+init=epsg:28992")
data(meuse.grid)
  coordinates(meuse.grid) = ~x+y
    proj4string(meuse.grid) <- CRS("+init=epsg:28992")
      gridded(meuse.grid) = TRUE    
        r <- raster(meuse.grid) 
          r[] <- runif(ncell(r))

# Create a polygon
f <- gBuffer(meuse[10,], byid=FALSE, id=NULL, width=250, 
                         joinStyle="ROUND", quadsegs=10)   

# Plot full raster and polygon                       
plot(r)
  plot(f,add=T)

# Crop using extent, rasterize polygon and finally, create poly-raster
#          **** This is the code that you are after ****  
cr <- crop(r, extent(f), snap="out")                    
  fr <- rasterize(f, cr)   
    lr <- mask(x=cr, mask=fr)

# Plot results
plot(lr)
  plot(f,add=T)
Jeffrey Evans
fuente
44
extracto () no realizar en la reproyección mosca, pero los cultivos () no lo hace. Eso podría explicar cierta confusión
mdsumner
@Jefferey crop () y mask () solo recortan el ráster de acuerdo con las extensiones rectangulares del polígono, no lo recorta dentro del límite del polígono. ¿Alguna idea de qué comandos podrían recortar el ráster dentro del límite del polígono dado?
csheth
1
@Chintan Sheth, para que la máscara se subconjunte dentro del polígono, debe tener un ráster que represente los valores dentro del polígono. Es por eso que rasteriza el polígono del subconjunto y luego lo enmascara. El paso de recorte es reducir la extensión del ráster al mismo que el polígono del subconjunto, que en el pasado, si se omite, arrojaría un error de desajuste de extensión.
Jeffrey Evans
spTransformdesde el sppaquete (que a veces se carga automáticamente con otros paquetes R espaciales) permite la reproyección para que ambos archivos estén en la misma proyección, por ejemplo. good_poly=spTransform(spolygon, CRSobj=crs(raster_file))
user3386170
@ user3386170, ¿eh? No estoy seguro de a qué te refieres. Esta pregunta se produjo en un momento en que el paquete ráster acaba de agregar "sobre la marcha la proyección" dentro de algunas de sus funciones. Estas funciones previamente arrojaron un error cuando las proyecciones no coincidían, pero esta publicación fue de 2014. También debe tener en cuenta que siempre carga rgdal cuando usa sp :: spTransform () ya que agrega funcionalidad adicional importante detrás de escena.
Jeffrey Evans
8

Lo que realmente busqué fue la mask()función.

mask(raster, poly_shape)

funciona sin errores y devuelve lo que busqué.

Lars
fuente
2
Vuelva a proyectar sus datos para que estén en el mismo espacio de proyección. Incluso en ArcGIS, donde la proyección sobre la marcha es automática, es una práctica muy mala realizar análisis en diferentes proyecciones. Con datos en diferentes proyecciones, sus datos no compartirán extensiones comunes, que es el error que está recibiendo.
Jeffrey Evans
Use projectExtent () para obtener solo la extensión para recortar el ráster.
mdsumner
Para adaptarse al formato de preguntas y respuestas del sitio, esto debe colocarse en el cuerpo principal de la pregunta como edición / actualización (y luego agregar un comentario a la respuesta a la que está en "respuesta" para que sepan que hay más para ver).
Matt wilkie
@mattwilkie Lo siento por no ajustarse al formato, pero mi texto era demasiado largo para publicarlo como comentario aquí. @JeffreyEvans He intentado lo siguiente: projection(raster) = projection(poly_shape)y al revés projection(poly_shape) = projection(raster), pero en ambos sentidos producir el mismo error: Error in .local(x, y, ...) : extents do not overlap In addition: Warning message: In intersect(extent(x), extent(y)) : Objects do not overlap. ¿Hay alguna manera de ver qué proyección se usa sobre la marcha usando la función extract () (porque obviamente funciona)?
Lars
1
Lo que realmente busqué fue la función mask (). mask(raster, poly_shape)funciona sin errores y devuelve lo que busqué.
Lars
-1

La extensión funciona bien ... Creo que Xmin, Xmax, Ymin e Ymax de su extensión son diferentes de las X e Y de su ráster, es decir, están opuestas

ToNoY
fuente