diff --git a/articles/wcs-ssurgo.html b/articles/wcs-ssurgo.html index fc091691..85cb6e03 100644 --- a/articles/wcs-ssurgo.html +++ b/articles/wcs-ssurgo.html @@ -311,6 +311,7 @@

Web Coverage Service # add original point, after transforming to mukey grid CRS plot(project(p, "EPSG:5070"), add = TRUE, pch = 16) +

Manual Creation of Bounding Boxes

@@ -355,6 +356,7 @@

Manual Creation of Bounding Boxes # add original BBOX, after transforming to mukey grid CRS plot(st_transform(x, 5070), add = TRUE)

+

@@ -388,6 +390,16 @@

Map Unit Key Grids# check: print(mu) +#> class : SpatRaster +#> dimensions : 147, 219, 1 (nrow, ncol, nlyr) +#> resolution : 30, 30 (x, y) +#> extent : -1365483, -1358913, 2869259, 2873669 (xmin, xmax, ymin, ymax) +#> coord. ref. : NAD83 / Conus Albers (EPSG:5070) +#> source(s) : memory +#> categories : mukey +#> name : mukey +#> min value : 144983 +#> max value : 1716001 plot( mu, @@ -396,6 +408,7 @@

Map Unit Key Grids= FALSE, legend = FALSE ) +

SSURGO Polygons from SDA

@@ -426,6 +439,7 @@

SSURGO Polygons from SDA legend = FALSE) plot(p, add = TRUE, border = 'white') mtext('CONUS Albers Equal Area Projection (EPSG:5070)', side = 1, line = 1)

+

Grid Resolution Specification @@ -461,6 +475,7 @@

Grid Resolution Specification= FALSE, legend = FALSE )

+

Raster Soil Survey Data @@ -482,12 +497,42 @@

Raster Soil Survey Data # gSSURGO grid: 30m resolution (x <- mukey.wcs(a, db = 'gSSURGO', res = 30)) +#> class : SpatRaster +#> dimensions : 267, 200, 1 (nrow, ncol, nlyr) +#> resolution : 30, 30 (x, y) +#> extent : 1129000, 1135000, 1403000, 1411010 (xmin, xmax, ymin, ymax) +#> coord. ref. : NAD83 / Conus Albers (EPSG:5070) +#> source(s) : memory +#> categories : mukey +#> name : mukey +#> min value : 545800 +#> max value : 545887 # gNATSGO grid: 30m resolution (y <- mukey.wcs(a, db = 'gNATSGO', res = 30)) +#> class : SpatRaster +#> dimensions : 267, 200, 1 (nrow, ncol, nlyr) +#> resolution : 30, 30 (x, y) +#> extent : 1129000, 1135000, 1403000, 1411010 (xmin, xmax, ymin, ymax) +#> coord. ref. : NAD83 / Conus Albers (EPSG:5070) +#> source(s) : memory +#> categories : mukey +#> name : mukey +#> min value : 545800 +#> max value : 3244759 # RSS grid: 10m resolution (z <- mukey.wcs(a, db = 'RSS', res = 10)) +#> class : SpatRaster +#> dimensions : 800, 600, 1 (nrow, ncol, nlyr) +#> resolution : 10, 10 (x, y) +#> extent : 1129000, 1135000, 1403000, 1411000 (xmin, xmax, ymin, ymax) +#> coord. ref. : NAD83 / Conus Albers (EPSG:5070) +#> source(s) : memory +#> categories : mukey +#> name : mukey +#> min value : 3244721 +#> max value : 3244759 # graphical comparison par(mfcol = c(1, 3)) @@ -520,6 +565,7 @@

Raster Soil Survey Data ext = x ) plot(a, add = TRUE)

+

STATSGO @@ -531,6 +577,16 @@

STATSGO
 (statsgo <- mukey.wcs(a, db = 'statsgo', res = 300))
+#> class       : SpatRaster 
+#> dimensions  : 27, 20, 1  (nrow, ncol, nlyr)
+#> resolution  : 300, 300  (x, y)
+#> extent      : 1129000, 1135000, 1403000, 1411100  (xmin, xmax, ymin, ymax)
+#> coord. ref. : NAD83 / Conus Albers (EPSG:5070) 
+#> source(s)   : memory
+#> categories  : mukey 
+#> name        :  mukey 
+#> min value   : 659074 
+#> max value   : 664845
 
 # graphical comparison
 par(mfcol = c(1, 2))
@@ -550,6 +606,7 @@ 
STATSGO= FALSE,
   main = attr(statsgo, 'layer name')
 )

+

Hawaii SSURGO @@ -578,6 +635,7 @@

Hawaii SSURGO= attr(mu, 'layer name'), colNA = 'royalblue' )

+

Puerto Rico SSURGO @@ -606,6 +664,7 @@

Puerto Rico SSURGO= attr(mu, 'layer name'), colNA = 'royalblue' )

+

@@ -641,6 +700,47 @@

Map Unit Aggregate Values # extract raster attribute table for thematic mapping (rat <- cats(mu2)[[1]]) +#> ID mukey +#> 1 144983 144983 +#> 2 144984 144984 +#> 3 144985 144985 +#> 4 144986 144986 +#> 5 145005 145005 +#> 6 145009 145009 +#> 7 145010 145010 +#> 8 145011 145011 +#> 9 145012 145012 +#> 10 145015 145015 +#> 11 145017 145017 +#> 12 145019 145019 +#> 13 145020 145020 +#> 14 145056 145056 +#> 15 145057 145057 +#> 16 145058 145058 +#> 17 145059 145059 +#> 18 145060 145060 +#> 19 145068 145068 +#> 20 145069 145069 +#> 21 145070 145070 +#> 22 145076 145076 +#> 23 145079 145079 +#> 24 145118 145118 +#> 25 145183 145183 +#> 26 145195 145195 +#> 27 145208 145208 +#> 28 145250 145250 +#> 29 145253 145253 +#> 30 145264 145264 +#> 31 145269 145269 +#> 32 145275 145275 +#> 33 145278 145278 +#> 34 145328 145328 +#> 35 145329 145329 +#> 36 145340 145340 +#> 37 145343 145343 +#> 38 145385 145385 +#> 39 1715935 1715935 +#> 40 1716001 1716001 # optionally use convenience function: # * returns all fields from muaggatt table @@ -657,6 +757,13 @@

Map Unit Aggregate Values # check head(tab) +#> mukey aws050wta aws0100wta +#> 1 144983 7.26 14.02 +#> 2 144984 7.31 14.14 +#> 3 144985 7.36 14.27 +#> 4 144986 7.06 13.19 +#> 5 145005 5.21 9.34 +#> 6 145009 5.49 9.16 # set raster categories levels(mu2) <- tab @@ -675,6 +782,7 @@

Map Unit Aggregate Values ), range = c(0, 20) )

+

Interpretations for Soil Suitability / Limitation @@ -702,6 +810,41 @@

Interpretations for Soi # check head(tab) +#> areasymbol musym muname mukey +#> 1 MT629 101 McCollum fine sandy loam, 0 to 2 percent slopes 144983 +#> 2 MT629 102 McCollum fine sandy loam, 2 to 4 percent slopes 144984 +#> 3 MT629 103 McCollum fine sandy loam, 4 to 8 percent slopes 144985 +#> 4 MT629 104 McCollum fine sandy loam, gravelly substratum, 0 to 2 percent slopes 144986 +#> 5 MT629 120 Niarada gravelly loam, 0 to 4 percent slopes 145005 +#> 6 MT629 123 Niarada gravelly loam, cool, 15 to 30 percent slopes 145009 +#> rating_ENGConstructionMaterialsRoadfill rating_AWMIrrigationDisposalofWastewater +#> 1 1.00 0.00 +#> 2 0.98 0.05 +#> 3 0.98 0.68 +#> 4 1.00 0.15 +#> 5 0.98 0.19 +#> 6 0.12 1.00 +#> class_ENGConstructionMaterialsRoadfill class_AWMIrrigationDisposalofWastewater +#> 1 Well suited Not limited +#> 2 Moderately well suited Slightly limited +#> 3 Moderately well suited Moderately limited +#> 4 Well suited Slightly limited +#> 5 Moderately well suited Slightly limited +#> 6 Poorly suited Very limited +#> reason_ENGConstructionMaterialsRoadfill +#> 1 <NA> +#> 2 Low strength; Dusty +#> 3 Low strength; Dusty +#> 4 <NA> +#> 5 Dusty; Dusty; Dusty +#> 6 Slope; Dusty; Depth to bedrock; Slope; Slope; Dusty +#> reason_AWMIrrigationDisposalofWastewater +#> 1 <NA> +#> 2 Filtering capacity; Too steep for surface application +#> 3 Too steep for surface application; Too steep for surface application +#> 4 Filtering capacity; Filtering capacity; Droughty +#> 5 Large stones on the surface; Droughty; Filtering capacity; Droughty; Slow water movement; Droughty +#> 6 Too steep for surface application; Too steep for sprinkler application; Large stones on the surface; Droughty; Droughty; Too steep for surface application; Too steep for sprinkler application; Depth to bedrock; Large stones on the surface; Filtering capacity; Too steep for surface application; Droughty; Too steep for sprinkler application; Too steep for surface application; Too steep for sprinkler application; Droughty # set ordered factor levels (for nice label/legend order) tab$class_ENGConstructionMaterialsRoadfill <- factor( @@ -727,6 +870,7 @@

Interpretations for Soi horizontal = TRUE, # fuzzy ratings on X axis las = 1 # rotate axis labels 90 degrees )

+

From above graph we can see that the different suitability rating classes class_ENGConstructionMaterialsRoadfill each correspond to a range of fuzzy values @@ -753,6 +897,7 @@

Interpretations for Soi 'Irrigation Disposal of Wastewater\nWeighted Mean over Components' ) ) +

Component-level Properties @@ -791,6 +936,7 @@

Component-level Properties axes = FALSE, legend = "topleft" )

+

Another example is thematic mapping of the “simplified component parent material group”. First, set up a new AOI for the following examples:

@@ -813,6 +959,7 @@

Component-level Properties cex.main = 0.7, main = 'gSSURGO Map Unit Key Grid' ) +

We use get_SDA_pmgroupname() to obtain the tabular parent material information to relate to map unit keys:

@@ -833,6 +980,7 @@ 
Component-level PropertiesactiveCat(mu2) <- 'pmgroupname'
 
 plot(mu2, legend = "topleft", axes = FALSE)
+

We can also inspect a mapunit-level hydric rating derived from the default aggregation method in get_SDA_hydric().

@@ -850,6 +998,7 @@ 
Component-level Properties# set active category 
 activeCat(mu2) <- 'HYDRIC_RATING'
 plot(mu2, legend = "topleft", axes = FALSE)
+

Several Horizon-level Soil Properties @@ -880,6 +1029,20 @@

Several Horizon-level Soil Proper # check head(tab) +#> areasymbol musym muname mukey +#> 1 OK113 1 Apperson silty clay loam, 1 to 3 percent slopes 623396 +#> 2 OK113 4 Coyle loam, 1 to 3 percent slopes 623399 +#> 3 OK113 7 Keokuk very fine sandy loam, 0 to 1 percent slopes, occasionally flooded 623402 +#> 4 OK113 10 Bethany silt loam, 1 to 3 percent slopes 623405 +#> 5 OK113 11 Bethany silt loam, 3 to 5 percent slopes 623406 +#> 6 OK113 12 Bethany-Pawhuska complex, 1 to 5 percent slopes 623407 +#> dbthirdbar_r awc_r ph1to1h2o_r +#> 1 1.45 0.18 6.15 +#> 2 1.40 0.18 5.70 +#> 3 1.43 0.17 7.30 +#> 4 1.35 0.20 6.20 +#> 5 1.30 0.20 6.30 +#> 6 1.34 0.20 6.30 # convert areasymbol into a factor easy plotting later tab$areasymbol <- factor(tab$areasymbol) @@ -889,6 +1052,7 @@

Several Horizon-level Soil Proper # list variables in the RAT names(cats(mu)[[1]]) +#> [1] "mukey" "dbthirdbar_r" "awc_r" "ph1to1h2o_r" # convert categories associated with keys to values mu2 <- catalyze(mu)

@@ -896,16 +1060,22 @@

Several Horizon-level Soil Proper 
 plot(mu2$awc_r)
 box()
+

Plot aggregate soil properties.

 plot(mu2[['dbthirdbar_r']], cex.main = 0.7,
-     main = '1/3 Bar Bulk Density (g/cm^3)\nDominant Component\n0-25cm')
-
+     main = '1/3 Bar Bulk Density (g/cm^3)\nDominant Component\n0-25cm')
+

+
+
 plot(mu2[['awc_r']], cex.main = 0.7,
-     main = 'AWC (cm/cm)\nDominant Component\n0-25cm')
-
+     main = 'AWC (cm/cm)\nDominant Component\n0-25cm')
+

+
+
 plot(mu2[['ph1to1h2o_r']], cex.main = 0.7,
      main = 'pH 1:1 H2O\nDominant Component\n0-25cm')
+

Sand, Silt, and Clay at a Soil Survey Area Boundary @@ -914,7 +1084,7 @@

Sand, Silt, and Clay areas. In this case the one soil survey was published in 1979 and the other in 2004.

First, we setup BBOX and query map unit key WCS.

-
+
 # extract a BBOX like this from SoilWeb by pressing "b"
 bb <- '-91.6853 36.4617,-91.6853 36.5281,-91.5475 36.5281,-91.5475 36.4617,-91.6853 36.4617'
 wkt <- sprintf('POLYGON((%s))', bb)
@@ -927,11 +1097,12 @@ 
Sand, Silt, and Clay
 
 # note SSA boundary
 plot(mu, legend = FALSE, axes = FALSE)

+

 
Then we derive aggregate sand, silt, clay (RV) values from the
 largest component, taking the weighted mean over 25-50cm depth interval.
 We also will take the sand and clay values to calculate the surface
 texture class for comparison.

-
+
 # extract RAT for thematic mapping
 rat <- cats(mu)[[1]]
 
@@ -950,6 +1121,20 @@ 
Sand, Silt, and Clay
 
 # check
 head(tab)
+#>   areasymbol musym                                                              muname   mukey
+#> 1      MO091 73306               Gressy-Gatewood complex, 3 to 8 percent slopes, rocky  691980
+#> 2      MO149 73321                       Alred-Gatewood complex, 1 to 8 percent slopes 2502332
+#> 3      MO149 73322                      Alred-Gatewood complex, 8 to 15 percent slopes 2502334
+#> 4      MO149 76002 Batcave-Farewell complex, 1 to 3 percent slopes, frequently flooded 2503322
+#> 5      MO149 76046             Secesh silt loam, 1 to 3 percent slopes, rarely flooded 2503473
+#> 6      MO149 76047    Secesh-Tilk complex, 1 to 3 percent slopes, occasionally flooded 2503476
+#>   sandtotal_r silttotal_r claytotal_r
+#> 1       23.00       59.00       18.00
+#> 2       28.69       51.11       20.20
+#> 3       28.69       51.11       20.20
+#> 4       40.00       40.00       20.00
+#> 5       25.65       49.00       25.35
+#> 6       26.91       48.08       25.02
 
 # set raster categories
 levels(mu) <- tab[, c('mukey', vars)]
@@ -977,6 +1162,7 @@ 
Sand, Silt, and Clay
   cex.main = 0.7,
   main = paste0(names(r), " - 25-50cm\nDominant Component")
 )

+