El nivel de las tasas aeroportuarias varía perceptiblemente de un Estado miembro a otro y a menudo incluso de un aeropuerto a otro en un mismo Estado miembro.

Question

Estoy trabajando con el procesado de datos LiDAR de áreas boscosas en ArcMap 10, con varios basa en el área de las estadísticas extraídas de cada uno de 30x30m polígono de la cuadrícula de células en un shapefile. Hay aproximadamente 20.000 de estas células. También tengo una capa de puntos relativos a los árboles individuales, lugares para todo el sitio.

Lo que tengo que hacer es calcular el promedio de la distancia entre el árbol y sus desviaciones estándar dentro de cada cuadrícula de células poligonales. Y, a continuación, escribir los dos valores a la correspondiente cuadrícula de células.

He utilizado el GME (caja de herramientas http://www.spatialecology.com/gme/) para obtener este tipo de datos para áreas correspondientes a mi campo de las ubicaciones del sitio; sin embargo, esto simplemente no es factible sobre un área tan grande.

¿Alguien tiene alguna idea de cómo puedo solucionar este problema?

Answer 1

Aquí está una R solución, la intención de funcionar como pseudocódigo para la aplicación en cualquier plataforma apropiada (C++, Python, etc) y para ser un prototipo de trabajo. Se inicia con una función para calcular la media y SD distancias de una matriz de puntos:

#
# Compute distance statistics for points (x,y).
#
dist.stats <- function(x,y) {
    #
    # Compute all distances between distinct points.
    #
    xy <- cbind(x,y)
    f <- function(w) apply(xy, 1, function(v) sum((v - w)^2))
    distances <- apply(xy, 1, f)
    distances <- sqrt(distances[lower.tri(distances)])
    #
    # Return the distance statistics.
    #
    c(mean(distances), sd(distances), length(x))
}

Esta función debe ser aplicado celda por celda. Aquí, supongo que las células están dispuestas en una cuadrícula paralelo a los ejes de coordenadas. Esto permite que los puntos se agrupan por medio de operaciones aritméticas. (Si ya han sido agrupados por polígono (en virtud de una unión espacial), el código sería más simple: dos líneas de división de las coordenadas x e y por el polígono de identificación y una tercera línea sería la de aplicar el bloque de estadísticas para cada grupo).

blockstats <- function(f,x,y,cellsize=1,n.cols=1,n.rows=1,origin.x=0,origin.y=0) {
    #
    # Split points by column.
    #
    cols <- factor(floor((x-origin.x)/cellsize) + 1, 1:n.cols)
    x.cols <- split(x, cols)
    y.cols <- split(y, cols)
    rows.cols <- split(rows, cols)
    #
    # Split columns by rows.
    #
    g <- function(z,g) split(z, factor(g, 1:n.rows))
    x.groups <- mapply(g, x.cols, rows.cols)
    y.groups <- mapply(g, y.cols, rows.cols)
    #
    # Apply summary function `f` to each group.
    #
    mapply(dist.stats, x.groups, y.groups)
}

Para ilustrar su uso, vamos a crear un pequeño conjunto de datos de ejemplo:

cellsize <- 30
n.rows <- 10
n.cols <- 20
n <- n.rows * n.cols * 5
set.seed(17)
x <- rnorm(n, mean=1/2, sd=1/4) * cellsize * n.cols
y <- rnorm(n, mean=1/2, sd=1/6) * cellsize * n.rows

Con estos puntos en la mano, el bloque de estadísticas son calculadas, convertidos a formato raster, y se representa:

system.time(results <- blockstats(dist.stats, x, y, cellsize, n.cols, n.rows))
r.mean <- matrix(results[1,], nrow=n.rows, ncol=n.cols)
r.sd <- matrix(results[2,], nrow=n.rows, ncol=n.cols)
r.n <- matrix(results[3,], nrow=n.rows, ncol=n.cols)

require("raster")
raster.plot <- function(r,s) {
    plot(raster(r, xmn=0, xmx=n.cols*cellsize, ymn=0, ymx=n.rows*cellsize), main=s)
    points(x, y, cex = min(1, 10/sqrt(n)))
}
par(mfrow=c(2,2))
raster.plot(r.mean, "Mean block distance")
raster.plot(r.sd, "SD block distance")
raster.plot(r.n, "# points")

El aumento de n.rows de 10 a 100 y n.cols de 20 a 200 simula la situación en la pregunta: alrededor de 100.000 puntos que cubre 20.000 células en un 30 m cellsize de la cuadrícula. El tiempo es de 20 segundos. Dedicado código compilado debería ser varios órdenes de magnitud más rápido, pero incluso a esta velocidad lenta puede ser adecuada para el problema.