Usando el Analista Espacial, podemos esperar que cada iteración de IDW, umbral y acumulación de resultados tome desde una fracción de segundo hasta unos pocos segundos. Multiplicar por 24 * 365 = 8760 horas en un año sugiere este cálculo llevará horas como mínimo.
Puede ser mucho más rápido. Más rápido no siempre es mejor, pero aquí -porque hay cuestiones importantes relativas al método de interpolación, ya que es poco probable que el IDW sea preciso para esos datos- será una verdadera ventaja poder realizar los cálculos rápidamente a una escala realista para ver los resultados y experimentar con otras opciones.
Para ilustrar esto y para estimar cuán rápido se puede hacer realmente, implementé el cálculo a escala completa en R usando datos simulados. (Es una cuestión simple leer sus datos en R y exportar la cuadrícula resultante en formato ESRI usando, por ejemplo, el Raster paquete.) Después de los preliminares habituales de conseguir los datos en el formato correcto (los valores se almacenan en una matriz de 60 por 8760 y las ubicaciones de los puntos en una matriz de 60 por 2 de coordenadas (x,y)), el manejo de los valores perdidos, y así sucesivamente, todo el cómputo se reduce a dos líneas de código:
thresh <- threshold * xy %*% one # Compute the denominators
a <- matrix(rowMeans(xy %*% z >= thresh), n.x, n.y)
xy ha precalculado los pesos de los IDW (no normalizados). one indica dónde no faltan datos y z es la matriz de todos los datos. Las dos multiplicaciones de la matriz xy %*% one y xy %*% z computar los denominadores y numeradores de todos los 8760 cálculos de IDW a la vez. La comparación en la segunda línea crea las cuadrículas de 8760 de unos y ceros y rowMeans los promedia para todo el año. El resultado, a es una única cuadrícula de dimensiones n.x y n.y .
Reduciendo este gran número de interpolaciones lineales a unas pocas operaciones de matriz, se puede lograr la mayor velocidad y eficiencia posibles. Además, este enfoque es general: literalmente cualquier El método de interpolación lineal (como el Kriging) puede ser implementado de esta manera.
En mi estación de trabajo (ejecutando esto en un solo núcleo, aunque es posible una paralelización masiva) se tarda poco menos de 60 segundos en procesar una cuarta parte de la red. Al reducir la resolución de la red se pierden pocos detalles y los datos de todo el año pueden ser procesados en segundos.
Para ilustrar, Aquí están los resultados para 60 puntos y 8760 horas en una cuadrícula de 125 por 310 (media resolución; 58 segundos de tiempo total). Para simplificar, utiliza todos los datos disponibles para la interpolación. Si se desea que los vecindarios locales más pequeños participen en un paso de precalculo (muy rápido), donde los pesos xy son computados, pero de otra manera nada más cambiaría.
![Plot]()
Las velocidades del viento durante la primera hora se trazan como círculos cuyas áreas son proporcionales a la velocidad media durante el año; los promedios van de 1 a 5,9 m/seg.
Hay un precio que pagar: este método precomputa algunos valores y los almacena en la RAM. Este cálculo de un cuarto de tamaño requería 6,75 GB de RAM. Por lo tanto, un cálculo de tamaño completo requeriría 27 GB de RAM. Sin embargo, si no se dispone de tanto, es sencillo tallar la cuadrícula en rectángulos más pequeños, hacer el cálculo para cada rectángulo por separado y hacer un mosaico de las piezas después. El tiempo total de cálculo no aumentaría de forma apreciable.
Sigue el código completo.
set.seed(23)
n.points <- 60
n.hours <- 24 * 365
n.x <- 620/2 # Number of columns in grid
n.y <- 250/2 # Number of rows in grid
threshold <- 5 # Wind speed threshold
#
# The data: point locations.
#
points <- matrix(rnorm(n.points*2, 10^5, 10^4), ncol=2)
points <- points[order(points[,1]), ] # Sort by x-coordinate for testing
#
# Expand the extent of the points by about 10% to define the extent of the grid.
#
ranges <- apply(points, 2, range)
d.r <- diff(ranges)
x <- seq(ranges[1,1] - .05*d.r[1], ranges[2,1] + 0.05*d.r[1], length.out=n.x)
y <- seq(ranges[1,2] - .05*d.r[2], ranges[2,2] + 0.05*d.r[2], length.out=n.y)
#
# Obtain the inverse squared distances to the points.
#
distance.to <- function(u, x, y, eps=10^(-8)) {
e <- eps * (max(x) - min(x) + max(y) - min(y))^2
return(1 / (e + (x-u[1])^2 + (y-u[2])^2))
}
system.time(
xy <- apply(points, 1,
function(u) outer(x, y, function(x0, y0) distance.to(u, x0, y0)))
)
#
# The data: an array of values, one per column.
#
p.x <- points[,1]
z.scale <- rgamma(n.points, 5, scale=((p.x - min(x))/10^4 + 2)/50)
z <- matrix(rgamma(n.points * n.hours, 10, scale=z.scale), nrow=n.points)
z[z > 10] <- NA # Introduce some missing values
#
# Perform the interpolation, thresholding, and averaging.
#
system.time({
one <- !is.na(z) # Indicates the valid values
z[is.na(z)] <- 0 # Replace missing with 0
thresh <- threshold * xy %*% one # Compute the denominators
a <- matrix(rowMeans(xy %*% z >= thresh), n.x, n.y)
})
#
# Plot the results.
#
z.sum <- rowSums(z)
z.max <- max(z.sum)
z.root <- sqrt(z.sum / z.max) * 3
filled.contour(x=x, y=y, a, color.palette=terrain.colors,
plot.axes={axis(1); axis(2);
points(points, col="Black", cex=z.root)},
main="Proportion of Year at or above 5 m/sec",
xlab="Easting")
