Пространственно-временная интерполяция в R или ArcGIS?

12

Я пытаюсь рассчитать среднее значение количества осадков по количеству точек, используя инструмент Обратное взвешенное расстояние в ArcGIS 9.3.

Моя проблема состоит в том, что: у каждой точки есть свой собственный временной ряд, поэтому процесс интерполяции должен быть в состоянии выполнять в течение всех лет (так сказать, итерация).

Ниже приведен пример таблицы атрибутов:

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

Кто-нибудь может показать мне, как это сделать?

Редактировать 1: Я наконец сделал это, используя код C ++, который требовал сетку маски ArcGIS, файлы данных и расположение всех точек.

Редактировать 2: я недавно использовал R, чтобы сделать эту задачу интерполяции. Вы можете использовать либо пакеты hydroTSM, gstatлибо spacetime. Несколько примеров ссылок ниже:

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

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

Изменить 3: добавлен рабочий пример ниже для будущих читателей

arcgis-desktop r interpolation spatial-statistics time-series Tung
источник
Это поможет? Временной ряд
Брэд Несом
Это можно сделать в R, но я думаю, что есть простой способ сделать это прямо в ArcMap. Все, что нужно OP - это перебирать отдельные переменные (годы) и вычислять интерполированный растр для каждой отдельной переменной. Тот факт, что значения в этом примере являются последовательными годами, не имеет значения.
Энди W
Спасибо за ваш ответ. На самом деле есть пакетный вариант, если щелкнуть правой кнопкой мыши по инструменту IDW, но все же это довольно утомительная работа, если у вас есть почасовые или ежедневные данные. KR
Тун
@thecatalyst - если пакетный инструмент IDW выполняет свою работу, вы должны опубликовать его в качестве ответа. Хотя это может быть утомительным, но если это случается нечасто (так как ежегодные оценки количества осадков нечасты), то есть мало причин для поиска других решений.
Энди W
@Andy: пакетный инструмент поможет, если у вас есть ограниченное количество, но у меня есть сотни данных, которые делают идею его использования немного нереальной. Я все еще ищу решение этой проблемы. KR
Тун

Ответы:

3

Я решил эту проблему, вставив в модель итератор «Выбор элементов». (В окне ModelBuilder в меню «Вставка» -> «Итераторы».)

Используйте ваше поле времени в качестве переменной group by. При этом модель будет повторяться по одному разу для каждого класса объектов.

Затем присоедините предпочитаемый инструмент интерполяции (сплайн, IDW и т. Д.) К выходу объекта из итератора. Запустите модель, отправьтесь в отпуск на несколько недель, и когда вы вернетесь, у вас будет столько сеток, сколько у вас будет временных точек в классе пространственных объектов.

Обратите внимание, что это решение предполагает наличие дискретных точек выборки по времени с полем даты или числовым полем, которое указывает одну временную точку для каждой записи в вашем наборе функций. Если вы используете формат «время начала» и «время окончания», это может быть не так просто.


источник
1
Кроме того, не забудьте использовать переменную «% n%» в имени выходного файла (или какой-либо другой способ создания уникального имени файла), иначе вы можете перезаписывать свой растр на каждой итерации. Для получения дополнительной информации см. Help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//… или просто Google «Примеры подстановки
TY. Приятно знать, что есть другой способ сделать это. Ура!
Тун
2

Кажется, что на эту ветку отвечает инструмент IDW, но если бы вы запрашивали и вводили начальный год, а затем перебирали поля года, используя встроенную переменную в построителе модели, тогда это был бы более элегантный способ обработки моделирования. ,

PS: Я согласен с @AndyW, что если вы решили это с помощью IDW, опубликуйте ответ самостоятельно, а затем "отметьте галочкой"

CDBrown
источник
1

Добавить собственное решение, используя Rданные о случайных осадках

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)

Преобразовать в объект sp

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

Добавьте систему пространственной привязки (SRS) или систему координат (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"

Конвертировать в UTM 10N

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

Гипотетические данные годовых осадков, полученные с использованием распределения Пуассона. 

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

Объединить Prec данные фрейма с Prec Shapefile

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

Объединить фрейм данных осадков с шейп-файлом осадков (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")))

Определите экстент для пространственной интерполяции. Увеличьте его на 4 км в каждом направлении

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

Создать желаемую сетку с разрешением 1 км

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

Интерполировать с использованием обратного взвешенного расстояния (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)

Результаты интерполяции участков

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
источник