¿Utilizar sampleRandom() desde un raster grande sin valores NA en R?

Question

Tengo un número de rásteres de tamaño variable que necesitan ser muestreados aleatoriamente con el valor de retorno siendo una matriz de x, y, y valor. El paquete raster sampleRandom(raster,n, na.rm=TRUE, xy=TRUE) lo hará bien la mayoría de las veces. Cuando funciona correctamente, esta función devuelve una matriz de valores no-NA para n pares de coordenadas. Cuando aparecen valores NA en la muestra, se descartan y se sustituyen por un valor no NA.

Sin embargo, para mis rásteres (el más pequeño tiene 4e^7 celdas y algunos tienen un alto porcentaje de valores NA), sampleRandom() devuelve una matriz sustancialmente menor que n pares de coordenadas. Presumiblemente, esto se debe a que los valores NA muestreados, no son reemplazados cuando se muestrean.

¿Por qué el sampleRandom ¿la función devuelve resultados incompletos en el ejemplo de datos del mundo real?

Como ha señalado correctamente @Radar, la documentación de los paquetes rasterizados indica: With argument na.rm=TRUE, the returned sample may be smaller than requested

Con esto, mi pregunta se convierte en; ¿cómo puedo dibujar un trabajo alrededor de esto y eficientemente dibujar muestra aleatoria de n ¿pares de coordenadas?

Ejemplo 1: esto funciona correctamente en la recuperación de una muestra aleatoria de n a partir de un raster más grande recortado y enmascarado por polígonos espaciales. devuelve una matriz de 2000 puntos cordiantes.

region1 <- rbind(c(0,0), c(50,0), c(50,50), c(20,20), c(0,0))
region2 <- rbind(c(50,0), c(80,0), c(100,50), c(60,40), c(80,20), c(50,0))
polys <- SpatialPolygons(list(Polygons(list(Polygon(region1)), "region1"),
                        Polygons(list(Polygon(region2)), "region2")))

r <- raster(ncol=1000, nrow=1000)
r[] <- runif(ncell(r),0,1)
extent(r) <- matrix(c(0, 0, 1000, 1000), nrow=2)

r_crop <- crop(r, extent(polys), snap="out", progress='text')
r_mask <- mask(r_crop, polys) 

plot(r_mask)
plot(polys, add=TRUE)

x <- sampleRandom(r_mask,2000, na.rm=TRUE, xy=TRUE)
nrow(x)

>[1] 2000

Ejemplo 2: El siguiente ejemplo es con datos reales que consisten en un raster universal (geo.r) de 2e^8 celdas y un subconjunto de polígonos espaciales (geo.poly) que contiene 1200 polígonos y es de menor extensión que geo.r. Este código resulta incorrectamente en una matriz de mucho menos de n Dependiendo de la muestra aleatoria, algunas ejecuciones producen una matriz de entre 3 y 117 pares de coordenadas que no son AN.

require(maptools)
Prj <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
modeling_areas_SHP <- "C:/.../modeling_areas_dissolve.shp"
geo.polys <- readShapePoly(modeling_areas_SHP, IDvar="area_ID", proj4string=CRS(Prj))
geo.poly <- modeling_areas[modeling_areas$area_ID == i,] #subset the shapefile

geo.r <- raster("C:/.../cost_raster")

geo.r_crop <- crop(geo.r, extent(geo.poly), snap="out", progress='text')
geo.r_mask <- mask(geo.r_crop, geo.poly, progress='text')

plot(geo.r_mask)
plot(geo.poly, add=TRUE)

x <- sampleRandom(geo.r_mask,2000, na.rm=TRUE, xy=TRUE)
nrow(x)

>[1] 117

Al menos para mí, los ejemplos anteriores son iguales salvo por el tamaño total de los rásteres y la complejidad de los polígonos; dos factores muy importantes. Obviamente, no puedo proporcionar los datos del mundo real debido al tamaño de los archivos, pero espero que el código sea suficiente.

¿Cómo puedo solucionarlo?

---

Utilicé este truco de trabajo, pero no es súper eficiente. Sin embargo, fue más eficiente que usar spsample() del paquete "sp".

micro_sample = 50000
tmp_rand_smple <- data.frame(x = numeric(0), y = numeric(0), layer = numeric(0))
            while(nrow(tmp_rand_smple) < micro_sample){
                tmp_smple <- data.frame(sampleRandom(geo.r_mask,10000, na.rm=TRUE, xy=TRUE)) ### 10k is an arbitrary chunk, loops until > micro_sample
                tmp_rand_smple <- rbind(tmp_rand_smple, tmp_smple)
                tmp_rand_smple <- unique(tmp_rand_smple[c("x", "y", "layer")]) # remove any duplicate coordinate pairs
            }
            tmp_rand_smple <- tmp_rand_smple[1:micro_sample,] # trim to length of micro_sample

Ejemplo 3: Este es un ejemplo del código anterior que puede reproducirse con el shapefile vinculado. En mi ordenador, este código no devuelve el número requerido de muestras aleatorias https://www.dropbox.com/s/7poaqcxju808arw/riverine_region_1.zip

Prj <- "+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
geo.poly <- readShapePoly("FILE LOCATION:/riverine_region_1", IDvar="area_ID", proj4string=CRS(Prj))  ## set file location

r <- raster(ncol=5202, nrow=8182)
r[] <- runif(ncell(r),0,1)
extent(r) <- matrix(c( 1533500, 592219.7, 1447689, 537662.6), nrow=2)

r_crop <- crop(r, extent(geo.poly), snap="out", progress='text')
r_mask <- mask(r_crop, geo.poly, progress='text') 

plot(r_mask)
plot(geo.poly, add=TRUE)

x <- sampleRandom(r_mask,2000, na.rm=TRUE, xy=TRUE)
nrow(x)

Answer 1

Parece que esto es un artefacto del sampleRandom paquete que está utilizando.

Si comprueba el documentación , afirma que:

Con el argumento na.rm=TRUE, la muestra devuelta puede ser menor que la solicitada

¿Muestreo aleatorio de raster usando R? puede proporcionarle una forma alternativa de realizar este análisis.

Answer 2

Simplemente rellene el número deseado de muestras aleatorias y luego vuelva a muestrear hasta el n correcto. Esto debería tener en cuenta los NA ocasionales que se producen y que posteriormente se eliminan con el argumento na.rm=TRUE.

    require(raster)
    # Create example data
    r1 <- raster(ncols=500, nrows=500, xmn=0)
      r1[] <- runif(ncell(r1))
    r2 <- raster(ncols=500, nrows=500, xmn=0)
      r2[] <- runif(ncell(r2))  
    r <- stack(r1,r2)

    # Sample size
    n=50

    # Random sample of raster  
    r.samp <- sampleRandom(r, size=(n+20), na.rm=TRUE, sp=FALSE, asRaster=FALSE) 
      dim( r.samp )[1]

   # Create a random sample of n size to subset r.samp
   #   (works with dataframe, matrix and sp objects)
   r.samp <- r.samp[sample( 1:dim(r.samp)[1], n),]
    dim ( r.samp )[1]

Si puedes leer el ráster en la memoria, un enfoque en sp sería utilizar rgdal para crear un SpatialGridDataFrame y luego coaccionarlo en un SpatialPointsDataFrame para que puedas eliminar fácilmente los NA y terminar con un objeto de punto de tu submuestra. A continuación, puede muestrear los siguientes rásteres utilizando este objeto de punto sp. El dataframe @data puede ser extraído y convertido en una matriz para sus propósitos.

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

n=50 # Number of random samples

# Read raster data using rgdal, results in SpatialGridDataFrame 
r <- readGDAL(system.file("external/test.ag", package="sp")[1])
  class(r)
    spplot(r, "band1")

# Coerce into SpatialPointsDataFrame    
r <- as(r, "SpatialPointsDataFrame")      

# remove NA's   
r@data <- na.omit(r.pts@data)
  plot(r, pch=20)

# Create random sample. Object is a SpatialPointsDataFrame     
r.samp <- r[sample(1:dim(r)[1], n),]
  plot(r.samp, pch=20, col="red", add=TRUE)   
    class(r.samp)

#  Use r.samp sp object for additional sampling 
#    Add extra column and coerce to raster stack
r2 <- readGDAL(system.file("external/test.ag", package="sp")[1])
  r2@data <- data.frame(r2@data, band2=runif(dim(r2)[1]) ) 
    r2 <- stack(r2)

# Extract raster values using r.samp object
r.samp@data <- data.frame(r.samp@data, band2=extract(r2[[2]], r.samp))
  str(r.samp@data)

Answer 3

He podido reproducir el problema en el tercer ejemplo.

Una solución para ello es utilizar los procedimientos incorporados. Hay varias opciones, pero un método conveniente es seleccionar cada celda de la cuadrícula de forma uniforme e independiente con una probabilidad lo suficientemente grande como para asegurar que se seleccionen al menos n=2000 (o lo que sea) celdas no nulas, pero no mucho más que eso. Esto se puede conseguir calculando la desviación estándar de la proporción de todas las celdas que se seleccionarán (que tiene una distribución binomial) y añadiendo un pequeño múltiplo de esa desviación estándar a la proporción deseada. Un múltiplo en torno a 6 prácticamente garantiza que se seleccionarán al menos n celdas. En el código de ejemplo que sigue, se seleccionaron 2020 celdas donde se necesitaban 2000.

Este método es un poco ineficiente en comparación con la repetición del sampleRandom procedimiento. Sin embargo, a diferencia de este último, este método muestrea sin de reemplazo.

---

Este código continúa en el contexto del Ejemplo 3 de la pregunta: utiliza el r_mask para la entrada y requiere la raster que debe cargarse para poder utilizar getValues .

set.seed(17)
n.sample <- 2000 # Number of non-null cells to sample
system.time({
  m <- dim(r_mask)[1]
  n <- dim(r_mask)[2]
  k <- sum(!is.na(getValues(r_mask))) # Number of non-null cells
  p <- n.sample / k                   # Proportion of them to be sampled
  pp <- p + 6*sqrt(p*(1-p)/(m*n))     # Proportion to request

  z <- matrix(runif(m*n) < pp, nrow=m)# Indicator of cells to select
  x <- unlist(apply(z, 2, function(col) (1:m)[col]))    # X-coordinates
  y <- unlist(apply(t(z), 2, function(row) (1:n)[row])) # Y-coordinates
  z <- getValues(r_mask)[z]           # Values
  i <- !is.na(z)                      # Indicator of non-null values
  a <- cbind(x[i], y[i], z[i])        # Result (with too many rows)
  print(dim(a))
  a <- a[1:n.sample, ]                # Remove any unneeded rows
})