More

Spatio-temporal interpolation in R or ArcGIS?

Spatio-temporal interpolation in R or ArcGIS?


I'm trying to calculate average rainfall value from a number of points using Inverse Weighted Distance tool in ArcGIS 9.3.

My problem is that: each point has it own time series, therefore the interpolation process should be able to carry out for all years (kind of iteration so to speak).

Following is a sample attribute table:

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…

Could anybody show me how to do that?


Edit 1: I finally did this by using C++ code which required ArcGIS mask grid, data files & locations of all the points.


Edit 2: I recently used R to do this interpolation task. You can use eitherhydroTSM,gstatorspacetimepackages. Few example links below:

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

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


Edit 3: Added a working example below for future readers


I solved this by inserting a "Feature Selection" iterator into a model. (In the ModelBuilder Window, under Insert->Iterators menu.)

Use your time field as your "group by" variable. By doing this, the model will iterate once for each time in your feature class.

Then attach your preferred interpolation tool (spline, IDW, whatever) to the feature output from the iterator. Run the model, go on vacation for a few weeks, and when you come back, you will have as many grids as you have time points in the feature class.

Note that this solution assumes you have discrete time sampling points with a date or numeric field that indicates a single time point for each record in your feature set. If you are using the "begin time" and "end time" format, it might not be so straight forward.


It seems that this thread is answered by the IDW tool, but if you were to request and input of the start year and then iterate through the year fields using an inline variable in model builder then this would be a more elegant way to handle the modelling.

PS: I agree with @AndyW that if you solved it using the IDW, post as an answer yourself and then "mark with the tick"


Add my own solution usingR& random precipitation data

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)

Convert to a sp object

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

Add a spatial reference system (SRS) or coordinate reference system (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"

Convert to UTM 10N

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

Hypothetical annual precipitation data generated using Poisson distribution.

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

Merge Prec data frame with Prec shapefile

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

Merge Precipitation data frame with Precipitation 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")))

Define the extent for spatial interpolation. Make it 4km larger on each direction

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

Create desired grid at 1km resolution

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

Interpolate using Inverse Distance Weighted (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)

Plot interpolation results

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 #> *      #> 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 , #> # PRCP2017.var1.pred , PRCP2017.var1.var , PRCP2018.x , #> # PRCP2018.y , PRCP2018.var1.pred , PRCP2018.var1.var 

Watch the video: TIMELAPSE OF THE FUTURE: A Journey to the End of Time 4K