¿Interpolación espacio-temporal en R o ArcGIS?

12

Estoy tratando de calcular el valor promedio de lluvia a partir de varios puntos usando la herramienta Distancia inversa ponderada en ArcGIS 9.3.

Mi problema es que: cada punto tiene su propia serie de tiempo, por lo tanto, el proceso de interpolación debería poder llevarse a cabo durante todos los años (tipo de iteración, por así decirlo).

A continuación se muestra una tabla de atributos de muestra:

ID X Y Name Rain1990 Rain1991 Rain1992 Rain1993 .... Rain2010

1 xx1 yy1 AA 1210 1189 1863 1269 ......  
2 xx2 yy2 BB 1492 1502 2187 1923 ......
......

¿Alguien podría mostrarme cómo hacer eso?

Edición 1: Finalmente hice esto usando el código C ++ que requería la cuadrícula de máscara ArcGIS, archivos de datos y ubicaciones de todos los puntos.

Edición 2: recientemente usé R para hacer esta tarea de interpolación. Puede utilizar cualquiera de hydroTSM, gstato spacetimepaquetes. Pocos enlaces de ejemplo a continuación:

http://spatial-analyst.net/wiki/index.php?title=Spatial_interpolation_exercises_%28NL%29

http://www.geostat-course.org/Topic_Bivand_2012

Edición 3: se agregó un ejemplo de trabajo a continuación para futuros lectores

arcgis-desktop r interpolation spatial-statistics time-series Tung
fuente
¿Esto ayudará? Series temporales
Brad Nesom
Se podría hacer en R, pero imagino que hay una manera simple de hacerlo directamente en ArcMap. Todo lo que el OP quiere es iterar a través de las variables separadas (años) y calcular el ráster interpolado para cada variable separada. El hecho de que los valores en este ejemplo sean años consecutivos no hace ninguna diferencia.
Andy W
Gracias por su respuesta. En realidad, hay una opción por lotes cuando se hace clic con el botón derecho en la herramienta IDW, pero aún así es un trabajo bastante tedioso si tiene datos por hora o por día. KR
Tung
@thecatalyst: si la herramienta IDW por lotes hace el trabajo, debe publicarla como respuesta. Aunque puede ser tedioso, si es poco frecuente (como las estimaciones anuales de lluvia son poco frecuentes), entonces hay pocas razones para buscar otras soluciones.
Andy W
@Andy: la herramienta por lotes ayudaría si tienes un número limitado pero tengo cientos de datos que hacen que la idea de usarla sea un poco poco realista. Todavía estoy buscando la solución de este problema. KR
Tung

Respuestas:

3

Resolví esto insertando un iterador de "Selección de características" en un modelo. (En la ventana ModelBuilder, en el menú Insertar-> Iteradores).

Use su campo de tiempo como su variable "agrupar por". Al hacer esto, el modelo iterará una vez por cada vez en su clase de entidad.

Luego, adjunte su herramienta de interpolación preferida (spline, IDW, lo que sea) a la salida de características del iterador. Ejecute el modelo, vaya de vacaciones durante unas semanas y, cuando regrese, tendrá tantas cuadrículas como puntos de tiempo en la clase de entidad.

Tenga en cuenta que esta solución supone que tiene puntos de muestreo de tiempo discreto con una fecha o un campo numérico que indica un punto de tiempo único para cada registro en su conjunto de características. Si está utilizando el formato "hora de inicio" y "hora de finalización", puede que no sea tan sencillo.


fuente
1
Además, no olvide utilizar la variable "% n%" en el nombre del archivo de salida (o alguna otra forma de generar un nombre de archivo único), o puede sobrescribir su ráster en cada iteración. Para obtener más información, consulte help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//... o simplemente google "Ejemplos de sustitución de variables en línea con las variables del sistema ModelBuilder"
TY. Es bueno saber que hay una forma diferente de hacerlo. ¡Salud!
Tung
2

Parece que este hilo es respondido por la herramienta IDW, pero si tuviera que solicitar e ingresar el año de inicio y luego recorrer los campos del año utilizando una variable en línea en el generador de modelos, esta sería una forma más elegante de manejar el modelado .

PD: Estoy de acuerdo con @AndyW en que si lo resolvió usando IDW, publíquelo como respuesta y luego "marque con la marca"

CDBrown
fuente
1

Agregar mi propia solución usando Rdatos de precipitación aleatorios

library(tidyverse)
library(sp) # for coordinates, CRS, proj4string, etc
library(gstat)
library(maptools)

# Coordinates of gridded precipitation cells
precGridPts <- ("ID lat long
                1 46.78125 -121.46875
                2 46.84375 -121.53125
                3 46.84375 -121.46875
                4 46.84375 -121.40625
                5 46.84375 -121.34375
                6 46.90625 -121.53125
                7 46.90625 -121.46875
                8 46.90625 -121.40625
                9 46.90625 -121.34375
                10 46.90625 -121.28125
                11 46.96875 -121.46875
                12 46.96875 -121.40625
                13 46.96875 -121.34375
                14 46.96875 -121.28125
                15 46.96875 -121.21875
                16 46.96875 -121.15625
                ")

# Read precipitation cells
precGridPtsdf <- read.table(text = precGridPts, header = TRUE)

Convertir a un objeto sp

sp::coordinates(precGridPtsdf) <- ~long + lat # longitude first

Agregue un sistema de referencia espacial (SRS) o un sistema de referencia de coordenadas (CRS). 

# CRS database: http://spatialreference.org/ref/epsg/
sp::proj4string(precGridPtsdf) <- sp::CRS("+proj=longlat +ellps=WGS84 +datum=WGS84")
str(precGridPtsdf)
#> Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
#>   ..@ data       :'data.frame':  16 obs. of  1 variable:
#>   .. ..$ ID: int [1:16] 1 2 3 4 5 6 7 8 9 10 ...
#>   ..@ coords.nrs : int [1:2] 3 2
#>   ..@ coords     : num [1:16, 1:2] -121 -122 -121 -121 -121 ...
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : chr [1:16] "1" "2" "3" "4" ...
#>   .. .. ..$ : chr [1:2] "long" "lat"
#>   ..@ bbox       : num [1:2, 1:2] -121.5 46.8 -121.2 47
#>   .. ..- attr(*, "dimnames")=List of 2
#>   .. .. ..$ : chr [1:2] "long" "lat"
#>   .. .. ..$ : chr [1:2] "min" "max"
#>   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
#>   .. .. ..@ projargs: chr "+proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0"

Convertir a UTM 10N

utm10n <- "+proj=utm +zone=10 ellps=WGS84"
precGridPtsdf_UTM <- spTransform(precGridPtsdf, CRS(utm10n))

Datos hipotéticos de precipitación anual generados mediante la distribución de Poisson. 

precDataTxt <- ("ID PRCP2016 PRCP2017 PRCP2018
                1 2125 2099 2203
                2 2075 2160 2119
                3 2170 2153 2180
                4 2130 2118 2153
                5 2170 2083 2179
                6 2109 2008 2107
                7 2109 2189 2093
                8 2058 2170 2067
                9 2154 2119 2139
                10 2056 2184 2120
                11 2080 2123 2107
                12 2110 2150 2175
                13 2176 2105 2126
                14 2088 2057 2199
                15 2032 2029 2100
                16 2133 2108 2006"
)

precData <- read_table2(precDataTxt, col_types = cols(ID = "i"))

Fusionar marco de datos Prec con archivo de forma Prec

precGridPtsdf <- merge(precGridPtsdf, precData, by.x = "ID", by.y = "ID")
precdf <- data.frame(precGridPtsdf)

Combinar el marco de datos de precipitación con el archivo de forma de precipitación (UTM)

precGridPtsdf_UTM <- merge(precGridPtsdf_UTM, precData, by.x = "ID", by.y = "ID")

# sample extent
region_extent <- structure(c(612566.169007975, 5185395.70942594, 639349.654465079, 
                             5205871.0782451), .Dim = c(2L, 2L), .Dimnames = list(c("x", "y"
                             ), c("min", "max")))

Defina la extensión de la interpolación espacial. Hazlo 4 km más grande en cada dirección

x.range <- c(region_extent[1] - 4000, region_extent[3] + 4000)
y.range <- c(region_extent[2] - 4000, region_extent[4] + 4000)

Cree la cuadrícula deseada a una resolución de 1 km

grd <- expand.grid(x = seq(from = x.range[1], to = x.range[2], by = 1000), 
                   y = seq(from = y.range[1], to = y.range[2], by = 1000))   

# Convert grid to spatial object
coordinates(grd) <- ~x + y
# Use the same projection as boundary_UTM
proj4string(grd) <- "+proj=utm +zone=10 ellps=WGS84 +ellps=WGS84"
gridded(grd) <- TRUE

Interpolar usando la distancia inversa ponderada (IDW)

idw <- idw(formula = PRCP2016 ~ 1, locations = precGridPtsdf_UTM, newdata = grd)  
#> [inverse distance weighted interpolation]

# Clean up
idw.output = as.data.frame(idw)
names(idw.output)[1:3] <- c("Longitude", "Latitude", "Precipitation")

precdf_UTM <- data.frame(precGridPtsdf_UTM)

Trazar resultados de interpolación

idwPlt1 <- ggplot() + 
  geom_tile(data = idw.output, aes(x = Longitude, y = Latitude, fill = Precipitation)) +
  geom_point(data = precdf_UTM, aes(x = long, y = lat, size = PRCP2016), shape = 21, colour = "red") +
  viridis::scale_fill_viridis() + 
  scale_size_continuous(name = "") +
  theme_bw() +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  theme(axis.text.y = element_text(angle = 90)) +
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 0))) 
idwPlt1

### Now looping through every year 
list.idw <- colnames(precData)[-1] %>% 
  set_names() %>% 
  map(., ~ idw(as.formula(paste(.x, "~ 1")), 
               locations = precGridPtsdf_UTM, newdata = grd)) 

#> [inverse distance weighted interpolation]
#> [inverse distance weighted interpolation]
#> [inverse distance weighted interpolation]

idw.output.df = as.data.frame(list.idw) %>% as.tibble()
idw.output.df

#> # A tibble: 1,015 x 12
#>    PRCP2016.x PRCP2016.y PRCP2016.var1.pred PRCP2016.var1.var PRCP2017.x
#>  *      <dbl>      <dbl>              <dbl>             <dbl>      <dbl>
#>  1    608566.   5181396.              2114.                NA    608566.
#>  2    609566.   5181396.              2115.                NA    609566.
#>  3    610566.   5181396.              2116.                NA    610566.
#>  4    611566.   5181396.              2117.                NA    611566.
#>  5    612566.   5181396.              2119.                NA    612566.
#>  6    613566.   5181396.              2121.                NA    613566.
#>  7    614566.   5181396.              2123.                NA    614566.
#>  8    615566.   5181396.              2124.                NA    615566.
#>  9    616566.   5181396.              2125.                NA    616566.
#> 10    617566.   5181396.              2125.                NA    617566.
#> # ... with 1,005 more rows, and 7 more variables: PRCP2017.y <dbl>,
#> #   PRCP2017.var1.pred <dbl>, PRCP2017.var1.var <dbl>, PRCP2018.x <dbl>,
#> #   PRCP2018.y <dbl>, PRCP2018.var1.pred <dbl>, PRCP2018.var1.var <dbl>
Tung
fuente