Interpolacja czasoprzestrzenna w R lub ArcGIS?

Próbuję obliczyć średnią wartość opadów z wielu punktów za pomocą narzędzia Odwrotna ważona odległość w ArcGIS 9.3.

Mój problem polega na tym, że: każdy punkt ma swoje własne szeregi czasowe, dlatego proces interpolacji powinien być możliwy przez wszystkie lata (rodzaj iteracji, że tak powiem).

Oto przykładowa tabela atrybutów:

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

Czy ktoś mógłby mi pokazać, jak to zrobić?

Edycja 1: W końcu zrobiłem to za pomocą kodu C ++, który wymagał siatki maski ArcGIS, plików danych i lokalizacji wszystkich punktów.

Edycja 2: Ostatnio użyłem R do wykonania tego zadania interpolacji. Można użyć hydroTSM, gstatlub spacetimepakietów. Kilka przykładowych linków poniżej:

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

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

Edycja 3: Dodano poniżej działający przykład dla przyszłych czytelników

Rozwiązałem to, wstawiając do modelu iterator „Wybór funkcji”. (W oknie ModelBuilder, w menu Wstaw-> Iteratory).

Użyj pola czasu jako zmiennej „grupuj według”. W ten sposób model będzie powtarzał się raz za każdym razem w klasie obiektów.

Następnie dołącz preferowane narzędzie do interpolacji (splajn, IDW, cokolwiek) do wyniku operacji z iteratora. Uruchom model, idź na wakacje na kilka tygodni, a kiedy wrócisz, będziesz mieć tyle siatek, ile punktów czasowych w klasie obiektów.

Zauważ, że to rozwiązanie zakłada, że ​​masz dyskretne punkty próbkowania z datą lub polem numerycznym, które wskazują pojedynczy punkt czasowy dla każdego rekordu w zestawie funkcji. Jeśli używasz formatu „czas rozpoczęcia” i „czas zakończenia”, może nie być to takie proste.

Wygląda na to, że na ten wątek odpowiada narzędzie IDW, ale jeśli poprosisz i wprowadzisz rok początkowy, a następnie przejdziesz przez pola roku za pomocą zmiennej wbudowanej w programie do tworzenia modeli, byłby to bardziej elegancki sposób obsługi modelowania .

PS: Zgadzam się z @AndyW, że jeśli rozwiązałeś to za pomocą IDW, opublikuj jako odpowiedź siebie, a następnie „zaznacz za pomocą kleszcza”

Dodaj własne rozwiązanie, używając Rdanych losowych opadów

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)

Konwertuj na obiekt sp

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

Dodaj system odniesienia przestrzennego (SRS) lub system odniesienia za pomocą współrzędnych (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"

Konwertuj na UTM 10N

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

Hipotetyczne dane o rocznych opadach wygenerowane przy użyciu rozkładu Poissona.

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"))

Scal ramkę danych Prec z plikiem kształtu Prec

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

Scal ramkę danych opadów z Preformitation shapefile (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")))

Określ zasięg interpolacji przestrzennej. Zwiększyć o 4 km w każdym kierunku

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

Utwórz żądaną siatkę w rozdzielczości 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

Interpoluj używając odwrotnej ważonej odległości (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)

Wykreśl wyniki interpolacji

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>