This tutorial demonstrates some simple spatial overlay analysis of polygon data using the IDEAS data model, as described in Robertson et al. 2020.
We will load some sample data from the stampr package, and pull out two polygons to demonstrate overlay operations.
library(stampr)
library(sp)
data(mpb)
P1 <- subset(mpb, TGROUP==1)[5,]
P2 <- subset(mpb, TGROUP==2)[7,]
plot(P2, border="green")
plot(P1, add=TRUE, border="blue")First we need to load some libraries;
We will use the con data connection to access a table called mpb which has the same data from the stampr package in IDEAS format.
mpb.i <- tbl(con,"MPB")
grid <- tbl(con,"FINALGRID2") %>% filter(RESOLUTION==19)
head(mpb.i)
#> # Source: lazy query [?? x 4]
#> # Database: NetezzaConnection
#> DGGID VALUE KEY TID
#> <dbl> <int> <chr> <int>
#> 1 4921607324 0 BOUNDARY 1264
#> 2 4921607324 0 ID 1264
#> 3 4921607324 1264 tid 1264
#> 4 4923988919 0 BOUNDARY 1264
#> 5 4924028283 0 BOUNDARY 1264
#> 6 4924185749 0 BOUNDARY 1264We want to pull out those same two polygons by identifying them by their ID values, as follows:
ID1 <- P1$ID
ID2 <- P2$ID
P1.i <- mpb.i %>% filter(KEY=="ID") %>% filter(VALUE==ID1) %>% inner_join(., grid, "DGGID") %>% mutate(WKT=inza..ST_AsText(GEOM)) %>% collect()
P2.i <- mpb.i %>% filter(KEY=="ID") %>% filter(VALUE==ID2) %>% inner_join(., grid, "DGGID") %>% mutate(WKT=inza..ST_AsText(GEOM)) %>% collect()
dbDisconnect(con)
plot(st_as_sf(P2.i, wkt='WKT', crs = 4326)['TID'], col='green', reset=FALSE)
plot(st_as_sf(P1.i, wkt='WKT', crs = 4326)['TID'], add=TRUE, col='blue')IDEAS data model
intersection <- P1.i %>% inner_join(., P2.i, "DGGID")
plot(st_as_sf(P2.i, wkt='WKT', crs = 4326)['TID'], col='green', reset=FALSE)
plot(st_as_sf(P1.i, wkt='WKT', crs = 4326)['TID'], add=TRUE, col='blue')
plot(st_as_sf(intersection, wkt='WKT.x', crs = 4326)['TID.x'], add=TRUE, col='red')