-
Notifications
You must be signed in to change notification settings - Fork 0
/
SSURGO-grid-overlay-examples.R
117 lines (89 loc) · 3.89 KB
/
SSURGO-grid-overlay-examples.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
library(soilDB)
library(sf)
library(mapview)
# new interface to leaflet
# https://github.com/r-spatial/leafem
library(leafem)
source('local-functions.R')
## TODO
# sf methods are strange...
# make a shiny app that can generate nested grids at a given point and overlay on SSURGO data
# use a projected CRS for non-rotated grid cells
# generate clumps of grid cells at an idealized MMU
# generate nested grids by feature
coords <- c(-119.7936, 36.7426)
aea <- 5070
q <- sprintf("SELECT mupolygongeo.STIntersection( geometry::STPointFromText('POINT(%f %f)', 4326).STBuffer(0.1).STEnvelope() ) AS geom, mukey
FROM mupolygon
WHERE mupolygongeo.STIntersects( geometry::STPointFromText('POINT(%f %f)', 4326).STBuffer(0.1) ) = 1
;", coords[1], coords[2], coords[1], coords[2])
# query / convert to SPDF
x <- SDA_query(q)
x <- processSDA_WKT(x)
# transform to CONUS AEA
x <- st_transform(x, aea)
# convert points to sf and CONUS AEA
# coordinates from matrix -> sf
# that is annoying
p <- st_as_sf(data.frame(t(coords)), coords=c(1,2), crs = 4326)
p <- st_transform(p, aea)
# # random points for vizualization
# s <- st_sample(x, size=50)
### ... why doesn't this work as expected?
# # get index to intersecting polygons
# idx <- which(lengths(st_intersects(x, s)) > 1)
#
# # generate grid from these polygons
# # g.800 <- st_make_grid((x[idx,]), cellsize = 800)
#
# makeGridByFeature <- function(i, size) {
# st_make_grid(x[i, ], cellsize = size)
# }
#
# g <- lapply(idx, makeGridByFeature, size=800)
# g.800 <- do.call('c', g)
## note: 30, 90, 270, 810 cell sizes can be nested
# ~ 113 ac. polygon near Fresno State
idx <- 26
# ~ 34 ac. polygon
idx <- 629
ng <- makeNestedGrids(x[idx, ])
cols <- viridis::viridis(5)
par(bg=grey(0.65))
plot(st_geometry(x[idx,]), border='white')
plot(st_geometry(ng[['10m']]), add=TRUE, border=cols[1])
plot(st_geometry(ng[['30m']]), add=TRUE, border=cols[2])
plot(st_geometry(ng[['90m']]), add=TRUE, border=cols[3], lwd=1)
plot(st_geometry(ng[['270m']]), add=TRUE, border=cols[4], lwd=3)
plot(st_geometry(ng[['810m']]), add=TRUE, border=cols[5], lwd=3)
# ## TODO: pick a single delineation, or iterate over a selection of delineations
# # full grid from bbox of SSURGO
# # not very efficient for demonstrations
# g.30 <- st_make_grid(x, cellsize = 30)
# g.90 <- st_make_grid(x, cellsize = 90)
# g.270 <- st_make_grid(x, cellsize = 270)
# g.800 <- st_make_grid(x, cellsize = 810)
#
#
# ## indexing via spatial intersection, not very efficient nor intuitive
# idx.30 <- which(lengths(st_intersects(g.30, st_buffer(s, 50))) > 0)
# idx.90 <- which(lengths(st_intersects(g.90, st_buffer(s, 100))) > 0)
# idx.270 <- which(lengths(st_intersects(g.270, st_buffer(s, 100))) > 0)
# idx.800 <- which(lengths(st_intersects(g.800, st_buffer(s, 100))) > 0)
#
# par(mar=c(1,1,1,1))
#
# plot(st_geometry(x[idx, ]))
# plot(st_geometry(s), pch=16, col='firebrick', cex=0.5, add=TRUE)
# plot(st_geometry(g.800[idx.800]), col=NA, border='royalblue', lty=1, add=TRUE)
# plot(st_geometry(g.270[idx.270]), col=NA, border='orange', lty=1, add=TRUE)
# plot(st_geometry(g.90[idx.90]), col=NA, border='yellow', lty=1, add=TRUE)
# plot(st_geometry(g.30[idx.30]), col=NA, border='white', lty=1, add=TRUE)
# mapview(x)
mv <- mapview(x[idx, ], lwd = 3, color = 'black', legend = FALSE, highlight = FALSE)
mv <- addFeatures(map = mv, data = st_transform(ng[['810m']], '+proj=longlat +datum=WGS84'), color='royalblue', weight=4)
mv <- addFeatures(map = mv, data = st_transform(ng[['270m']], '+proj=longlat +datum=WGS84'), color='firebrick', weight=3)
mv <- addFeatures(map = mv, data = st_transform(ng[['90m']], '+proj=longlat +datum=WGS84'), color='orange', weight=2)
mv <- addFeatures(map = mv, data = st_transform(ng[['30m']], '+proj=longlat +datum=WGS84'), color='yellow', weight=1)
mv <- addFeatures(map = mv, data = st_transform(ng[['10m']], '+proj=longlat +datum=WGS84'), color='white', weight=1)
mv