Interpolação espaço-temporal em R ou ArcGIS?

12

Estou tentando calcular o valor médio da precipitação de vários pontos usando a ferramenta Distância inversa ponderada no ArcGIS 9.3.

Meu problema é o seguinte: cada ponto tem sua própria série temporal, portanto, o processo de interpolação deve ser capaz de realizar durante todos os anos (tipo de iteração, por assim dizer).

A seguir está uma tabela de atributos de amostra:

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 ......
......

Alguém poderia me mostrar como fazer isso?

Edit 1: Finalmente fiz isso usando o código C ++, que exigia a grade de máscara do ArcGIS, os arquivos de dados e a localização de todos os pontos.

Edição 2: Recentemente usei R para executar esta tarefa de interpolação. Você pode usar tanto hydroTSM, gstatou spacetimepacotes. Alguns exemplos de links abaixo:

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

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

Edit 3: Adicionado um exemplo de trabalho abaixo para futuros leitores

arcgis-desktop r interpolation spatial-statistics time-series Tung
fonte
Isso vai ajudar? Série temporal
Brad Nesom 22/11
Isso poderia ser feito no R, mas imagino que haja uma maneira simples de fazê-lo diretamente no ArcMap. Tudo o que o OP deseja é percorrer as variáveis ​​separadas (anos) e calcular a varredura interpolada para cada variável separada. O fato de os valores neste exemplo serem anos sequenciais não faz diferença.
Andy W
Obrigado pela sua resposta. Na verdade, existe uma opção de lote quando você clica com o botão direito do mouse na ferramenta IDW, mas ainda assim é um trabalho tedioso se você tiver dados horários ou diários. KR
Tung
@ thecatalyst - Se a ferramenta IDW em lote executar o trabalho, você deve publicá-la como resposta. Embora possa ser entediante, se não for frequente (como as estimativas anuais de precipitação são pouco frequentes), há poucas razões para procurar outras soluções.
Andy W
@ Andy: A ferramenta de lote ajudaria se você tiver um número limitado, mas eu tenho centenas de dados que tornam a idéia de usá-lo um pouco irrealista. Ainda estou procurando a solução desse problema. KR
Tung

Respostas:

3

Resolvi isso inserindo um iterador "Feature Selection" em um modelo. (Na janela ModelBuilder, no menu Inserir-> Iteradores.)

Use seu campo de tempo como sua variável "agrupar por". Ao fazer isso, o modelo irá repetir uma vez a cada vez na sua classe de recursos.

Em seguida, conecte sua ferramenta de interpolação preferida (spline, IDW, qualquer que seja) à saída do recurso do iterador. Execute o modelo, tire férias por algumas semanas e, quando voltar, você terá tantas grades quantas pontos de tempo na classe de recursos.

Observe que esta solução pressupõe que você tenha pontos de amostragem de tempo discretos com uma data ou campo numérico que indica um único ponto de tempo para cada registro no seu conjunto de recursos. Se você estiver usando o formato "horário de início" e "horário de término", talvez não seja tão simples.


fonte
1
Além disso, não se esqueça de usar a variável "% n%" no nome do arquivo de saída (ou alguma outra maneira de gerar um nome de arquivo exclusivo); caso contrário, você poderá sobrescrever a varredura a cada iteração. Para mais informações, consulte help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//... ou apenas google "Exemplos de substituição de variáveis em linha com as variáveis do sistema ModelBuilder"
TY. É bom saber que existe uma maneira diferente de fazer isso. Felicidades!
Tung
2

Parece que esse encadeamento é respondido pela ferramenta IDW, mas se você solicitar e inserir o ano inicial e iterar pelos campos do ano usando uma variável embutida no construtor de modelos, essa seria uma maneira mais elegante de lidar com a modelagem .

PS: Concordo com o @AndyW que se você o resolveu usando o IDW, poste como resposta e depois "marque com o tiquetaque"

CDBrown
fonte
1

Adicionar minha própria solução usando R& dados aleatórios de precipitação

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)

Converter em um objeto sp

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

Adicione um sistema de referência espacial (SRS) ou sistema de referência 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"

Converta para UTM 10N

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

Dados hipotéticos anuais de precipitação gerados usando a distribuição 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"))

Mesclar o quadro de dados Prec com o shapefile Prec

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

Mesclar quadro de dados de precipitação com o shapefile de precisão (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 a extensão da interpolação espacial. Torne-o 4 km maior em cada direção

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

Crie a grade desejada com resolução 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 distância 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)

Resultados da interpolação do gráfico

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
fonte