¿Cómo producir una cuadrícula espacial a partir de la trama?

9

Necesito una cuadrícula espacial como una cuadrícula maestra para diversos mapas temáticos. ¿Cómo produzco una cuadrícula espacial a partir de un ráster que descarta todos los píxeles de NA?

Janbob Squarebrains
fuente
66
Tíranos algunas sobras aquí. ¿Un pequeño código para hacer una trama y qué esperas de él?
Spacedman

Respuestas:

14

Puede obtener todas las coordenadas de celdas que no son NA en un ráster con:

r = raster(matrix(runif(20),5,4))
r[r>.5]=NA

coordinates(r)[!is.na(values(r)),]
          x   y
 [1,] 0.375 0.7
 [2,] 0.125 0.5
 [3,] 0.375 0.5
 [4,] 0.625 0.5
 [5,] 0.875 0.5
 [6,] 0.125 0.3
 [7,] 0.375 0.3
 [8,] 0.625 0.3
 [9,] 0.375 0.1
[10,] 0.875 0.1

esas son las células que no son NA. Entonces, probablemente pueda alimentarlos a SpatialPixels

SpatialPixels(SpatialPoints(coordinates(r)[!is.na(values(r)),]))
Object of class SpatialPixels
Grid topology:
  cellcentre.offset cellsize cells.dim
x             0.125     0.25         4
y             0.100     0.20         4
SpatialPoints:
          x   y
 [1,] 0.375 0.7
 [2,] 0.125 0.5
 [3,] 0.375 0.5
 [4,] 0.625 0.5
 [5,] 0.875 0.5
 [6,] 0.125 0.3
 [7,] 0.375 0.3
 [8,] 0.625 0.3
 [9,] 0.375 0.1
[10,] 0.875 0.1
Coordinate Reference System (CRS) arguments: NA 

Aunque personalmente cualquier cosa en una cuadrícula lo guardaría como una trama.

Todavía no estoy totalmente seguro de qué es lo que quiere: los SpatialGridobjetos son cuadrículas rectangulares completas, por lo que uno sin los píxeles de NA no tiene sentido.

Hombre espacial
fuente
gracias amigo, tampoco estoy seguro de qué querer. Intento desarrollar un flujo de trabajo para producir mapas ráster idénticos (en términos de resolución, coordenadas, crs, ...) a partir de varios datos de puntos espaciales. rasterizar no es una opción debido a la dispersión de puntos espaciales. Necesita usar iwd o similar. Tiene razón, la cuadrícula es rectengular: lo que necesito como plantilla maestra probablemente sea un SpatialPointsDataFrame. que puede ser cuadriculado S.
Janbob Squarebrains
Probablemente desee un ráster maestro (¡ja!) Y luego crear todas las cuadrículas posteriores concurrentes.
Radar
2
Convenido. Si las cosas están cuadriculadas, la trama es casi siempre la respuesta. O eso, o pasear por unas 20 páginas de Análisis de datos espaciales aplicados en R para descubrir cómo hacer SpatialPixels. Hoy le di esa opción a un estudiante de doctorado, de hecho. Él eligió ... sabiamente!
Spacedman
20

Para transformar un RasterLayer en un objeto espacial * (cuadrícula o píxeles), puede usar la función de coerción "como"

library(raster)
r <- raster(matrix(runif(20),5,4))
r[r>.5] <- NA
g <- as(r, 'SpatialGridDataFrame')
p <- as(r, 'SpatialPixels')

plot(r)
points(p)
Robert Hijmans
fuente
7

Sus dos requisitos parecen ser sobre cosas diferentes:

1) Algún tipo de plantilla de cuadrícula ráster confiable.

2) Una cuadrícula dispersa que no almacena explícitamente las celdas faltantes.

sp :: GridTopology proporciona la primera, es solo una descripción de la cuadrícula basada en la coordenada de la celda inferior izquierda (cellcentre.offset), el espaciado de la celda (cellize) y las dimensiones de la cuadrícula (cells.dim).

La clase sp :: SpatialPixelsDataFrame le permite almacenar cuadrículas dispersas, pero por sí solo almacena mucho más que la "plantilla"; también almacena cada coordenada explícitamente. Esto se debe a que hace dos trabajos, uno le permite preservar las coordenadas originales que provienen de la cuadrícula y posiblemente son ligeramente irregulares, dos le permite almacenar solo aquellas celdas que tienen valores válidos. (Posiblemente * estos dos objetivos deberían haberse separado, pero esa es otra historia).

No creo que el paquete ráster tenga un análogo explícito a GridTopology, pero puede obtener los componentes para "rodar el suyo":

library(raster)
r1 <- raster(nrows=108, ncols=21, xmn=0, xmx=10)

## "cellsize"
res(r1)
## [1] 0.4761905 1.6666667

## extreme cell corners (just a different convention to sp's cellcentre)
extent(r1)
class       : Extent 
xmin        : 0 
xmax        : 10 
ymin        : -90 
ymax        : 90 

## we can also use bbox to get the same thing
bbox(r1)
min max
s1   0  10
s2 -90  90

## grid dimensions (including number of attributes/layers as 3rd "dim")
dim(r1)
## [1] 108  21   1

Al vincular todo esto, podemos pasar de raster a sp:

GridTopology(bbox(r1)[,1] + res(r1)/2, res(r1), dim(r1)[2:1])

(Observe cómo se deben revertir las dimensiones). Otra forma más simple es coaccionar a SpatialGrid y usar getGridTopology de sp, aunque esto es más costoso ya que las coordenadas terminan siendo generadas en el camino:

getGridTopology(as(r1, "SpatialGrid"))

Esas tres partes de la "topología" de ráster no son todas necesarias, ya que algunas son redundantes, pero de otra manera no hay una manera de capturarlas como un solo objeto, excepto que el ráster creado arriba está "vacío" y así puede hacer el trabajo que hace GridTopology para sp. No estoy seguro de los detalles de cuán "vacío" está, pero ciertamente no almacena explícitamente el espacio de "datos" y es más pequeño de lo que sería con valores en él. En general, el paquete ráster hace todo lo posible para mantener el uso de memoria al mínimo, por lo que es posible que no tenga que preocuparse por ser realmente "escaso".

Eso podría ayudar a explicar un poco más, sé que estoy superponiendo la respuesta de Spacedman, pero todavía no está claro exactamente lo que quieres decir en la pregunta todavía.

  • (Podría decirse que, dado que puede almacenar celdas dispersas simplemente almacenando su índice en lugar de coordenadas explícitas, y originalmente las coordenadas de celda "ligeramente irregulares" solo podrían almacenarse como atributos si realmente las quisiera. Ni el sp ni el ráster tratan con rásteres irregulares, ni siquiera el simple caso rectilíneo: esto está de acuerdo con la mayoría de las herramientas SIG, pero es desafortunado ya que es bastante común en formatos como NetCDF y se maneja con la buena función de gráficos :: imagen de R (aunque no si se usa la opción rasterImage reciente) .)
mdsumner
fuente
Gracias por las declaraciones muy instructivas, que necesito digerir. ¡Espero que me permita resolver mi problema o publicar una pregunta que sea más explícita!
Janbob Squarebrains
También @Spacedman: definí con mayor precisión el fondo de mi pregunta aquí en el foro.
Janbob Squarebrains