Soil Data Access

What is SDA?

SDA provides access to a SQL Server database containing a variety of soil information:

• Tabular data: map units, components, horizons, interpretations, ecological classes, laboratory characterization results, and more
• Spatial data: map unit polygons with geometry (available via WKT or for spatial queries), point locations (XY coordinates)

The SDA REST endpoint is available at:

https://sdmdataaccess.sc.egov.usda.gov/tabular/post.rest

All queries to SDA are executed through the SDA_query() function in soilDB, which handles the HTTP POST requests and JSON response parsing automatically.

SSURGO and STATSGO2 Data

SDA contains both the detailed Soil Survey Geographic Database (SSURGO) and the generalized U.S. General Soil Map (STATSGO2).

• SSURGO: Detailed mapping (scales 1:12,000 to 1:24,000). Identified by specific area symbols (e.g., CA630, IN005).
• STATSGO2: Generalized mapping (scale 1:250,000) for regional/national analysis. Identified by the single area symbol 'US'.

See the Data Filters section below for details on filtering out STATSGO2 data.

When to Use SDA vs. Local SSURGO

Use SDA for current data from NRCS without manual downloads, queries across multiple survey areas, quick one-time analyses, or access to interpretations and aggregated properties. Use local SSURGO databases for repeated queries to the same survey areas (faster), offline access, complete survey databases with all tables, or direct SQLite manipulation.

The same get_SDA_*() functions documented in this vignette can be applied to local SSURGO databases by passing the dsn parameter. See the Local SSURGO vignette for downloading and creating local SSURGO GeoPackage databases, and for examples of using these functions with local data.

Package Requirements

library(soilDB)

For spatial data examples, use sf and terra, which are both optional but recommended.

install.packages(c("sf", "terra"))

Core Syntax

The SDA_query() function executes SQL queries directly against the SDA database:

# Simple query to get the first 5 map units from an area
q <- "
SELECT TOP 5
  mapunit.mukey,
  mapunit.nationalmusym,
  legend.areasymbol,
  mapunit.musym,
  mapunit.muname
FROM legend
INNER JOIN mapunit ON mapunit.lkey = legend.lkey
WHERE legend.areasymbol = 'CA630'
ORDER BY mapunit.musym
"

result <- SDA_query(q)
head(result)

Handling Errors

If SDA_query() errors for any reason, the result will not be a data.frame with data. Instead it will be an object of class "try-error". You can handle the possibility of error using inherits(). For example:

result <- SDA_query(q)

# check if the result is an error response
if (!inherits(result, 'try-error')){
  head(result)
} else {
  # stop with error message extracted from result
  stop("Error occured! ", result[1])
}

You can set up a loop to retry a query that fails in the event of intermittent network or server issues. You should have a brief period of wait time (Sys.sleep()) between requests.

For repeated failures, it is best to increase that wait time in between each try up to, for example, a maximum of 3 attempts. The following example uses “exponential backoff” which increases the wait time in between requests incrementally from ~3 seconds to ~20 seconds over the course of several attempts.

n_tries <- 1
while (n_tries <= 3){
  result <- SDA_query(q)
  
  # check if the result is an error response
  if (!inherits(result, 'try-error')){
    head(result)
    
    # if no error, break out of the loop
    break
  } else {
    # on error, increment the number of tries
    # and sleep (time in seconds)
    Sys.sleep(exp(n_tries))
    
    n_tries <- n_tries + 1
    if (n_tries > 3) {
      stop("Error occured! ", result[1])
    }
  }
}

Understanding SDA SQL

SDA uses Microsoft SQL Server Transact-SQL (T-SQL) dialect.

Key differences of T-SQL from SQLite and other SQL dialects that are relevant to queries used in soilDB include:

When using SDA_query() with a local SQLite database (via dsn parameter), soilDB automatically converts some T-SQL syntax to SQLite equivalents to allow for both to work. This behavior is triggered when the dsn argument is specified. You can view the queries (and save them to run yourself) by passing the query_string=TRUE argument.

Basic Table Joins

Many custom SSURGO queries follow a standard pattern of relationships from legend, to mapunit, to component, to a related table to component such as component horizon.

In SSURGO data model primary and foreign key physical column names all end with the word "key". The primary key is a unique identifier of a particular record.

Foreign keys are references to primary keys of other tables and define the relationships between tables that can be used for joins.

Here we select a handful of values from all of the above mentioned tables corresponding to the first 10 component horizon (chorizon) in the SDA table.

# Get map unit, component, and horizon data
q <- "
SELECT TOP 10
  mu.mukey,
  leg.areasymbol,
  mu.musym,
  co.compname,
  co.comppct_r,
  ch.hzname,
  ch.hzdept_r,
  ch.hzdepb_r,
  ch.claytotal_r
FROM legend AS leg
INNER JOIN mapunit AS mu ON mu.lkey = leg.lkey
INNER JOIN component AS co ON co.mukey = mu.mukey
LEFT JOIN chorizon AS ch ON ch.cokey = co.cokey
WHERE leg.areasymbol = 'CA630'
  AND co.majcompflag = 'Yes'
  AND ch.claytotal_r IS NOT NULL
ORDER BY mu.musym, co.comppct_r DESC, ch.hzdept_r
"

result <- SDA_query(q)
head(result, 10)

We supplied all of the filtering constraints in a WHERE clause for clarity, with joins only using primary and foreign keys.

In general, the same constraints can be expressed as part of a join to relevant tables, which can be more efficient in complex queries.

Understanding LEFT JOIN vs INNER JOIN

One of the most important considerations in SSURGO queries is data completeness: not all relationships in the SSURGO database are populated for all soils or components. Understanding when to use LEFT JOIN vs. INNER JOIN is critical for correct analysis.

# These relationships are always present in SSURGO:
# - mapunit always belongs to a legend (survey area)
# - component always belongs to a mapunit
q <- "
SELECT TOP 10
  leg.areasymbol,
  mu.musym,
  co.compname,
  co.comppct_r
FROM legend AS leg
INNER JOIN mapunit AS mu ON mu.lkey = leg.lkey
INNER JOIN component AS co ON co.mukey = mu.mukey
WHERE leg.areasymbol = 'CA630'
ORDER BY mu.musym, co.comppct_r DESC
"

result <- SDA_query(q)

# Every row represents a valid component with all parent relationships
head(result)

# Some components have NO horizons defined
# INNER JOIN would exclude these components entirely
# LEFT JOIN preserves components even without horizon data
q <- "
SELECT
  leg.areasymbol,
  mu.mukey,
  mu.musym,
  co.compname,
  co.comppct_r,
  ch.chkey,
  ch.hzdept_r,
  ch.hzdepb_r
FROM legend AS leg
INNER JOIN mapunit AS mu ON mu.lkey = leg.lkey
INNER JOIN component AS co ON co.mukey = mu.mukey
LEFT JOIN chorizon AS ch ON ch.cokey = co.cokey
WHERE leg.areasymbol = 'CA630'
  AND co.majcompflag = 'Yes'
ORDER BY mu.musym, co.comppct_r DESC, ch.hzdept_r
"

result <- SDA_query(q)

# Check for components without horizons:
# Rows with NA in horizon columns indicate a component with no horizon data
if (is.data.frame(result) && nrow(result) > 0) {
  components_no_horizons <- result[is.na(result$chkey), ]
  if (nrow(components_no_horizons) > 0) {
    cat('Found', nrow(components_no_horizons), 'components without horizon data:\n')
    print(components_no_horizons[, c('musym', 'compname', 'comppct_r')])
  }
}

# Using INNER JOIN on horizons drops components that lack horizon data
q_inner <- "
SELECT
  leg.areasymbol,
  mu.mukey,
  mu.musym,
  co.compname,
  co.comppct_r,
  ch.desgnmaster,
  ch.hzdept_r,
  ch.hzdepb_r
FROM legend AS leg
INNER JOIN mapunit AS mu ON mu.lkey = leg.lkey
INNER JOIN component AS co ON co.mukey = mu.mukey
INNER JOIN chorizon AS ch ON ch.cokey = co.cokey
WHERE leg.areasymbol = 'CA630'
  AND co.majcompflag = 'Yes'
ORDER BY mu.musym, co.comppct_r DESC, ch.hzdept_r
"

result_inner <- SDA_query(q_inner)

# number of rows in LEFT JOIN result
nrow(result)

# number of rows in INNER JOIN result
nrow(result_inner)

# Some horizons have NO rock fragment data recorded
# LEFT JOIN preserves horizons even without fragment information
q <- "
SELECT
  mu.musym,
  co.compname,
  ch.desgnmaster,
  ch.hzdept_r,
  ch.hzdepb_r,
  rf.fragsize_h,
  rf.fragvol_r
FROM legend AS leg
INNER JOIN mapunit AS mu ON mu.lkey = leg.lkey
INNER JOIN component AS co ON co.mukey = mu.mukey
LEFT JOIN chorizon AS ch ON ch.cokey = co.cokey
LEFT JOIN chfrags AS rf ON rf.chkey = ch.chkey
WHERE leg.areasymbol = 'CA630'
  AND co.majcompflag = 'Yes'
ORDER BY mu.musym, co.comppct_r DESC, ch.hzdept_r
"

result <- SDA_query(q)

# Check for horizons without fragments:
if (is.data.frame(result) && nrow(result) > 0) {
  horizons_no_frags <- result[!duplicated(result$chkey) & is.na(result$fragsize_h), ]
  if (nrow(horizons_no_frags) > 0) {
    cat('Found', nrow(horizons_no_frags), 'horizons without rock fragment data:\n')
    print(horizons_no_frags[, c('musym', 'compname', 'desgnmaster', 'fragvol_r')])
  }
}

Aggregation Functions

SDA supports standard SQL aggregate functions (COUNT, SUM, AVG, MIN, MAX, etc.):

# Get component statistics per map unit
q <- "
SELECT
  mu.mukey,
  mu.musym,
  mu.muname,
  COUNT(co.cokey) AS n_components,
  SUM(co.comppct_r) AS total_comppct,
  AVG(co.comppct_r) AS avg_comppct
FROM legend AS leg
INNER JOIN mapunit AS mu ON mu.lkey = leg.lkey
INNER JOIN component AS co ON co.mukey = mu.mukey
WHERE leg.areasymbol = 'CA630'
GROUP BY mu.mukey, mu.musym, mu.muname
ORDER BY mu.musym
"

result <- SDA_query(q)
head(result)

Here we use SQL GROUP BY to define the columns used to determine unique groups. Other columns that are not part of the GROUP BY statement need to be included in some sort of aggregation expression to be used in SELECT clause.

When calculating map unit weighted averages using custom SQL, be aware that comppct_r (Component Percentage) does not always sum to 100% for a map unit due to minor component inclusions that may not be exported or missing data.

A robust weighted average should normalize weights by the sum of component percentages present in the data, rather than assuming a sum of 100. This is especially important if you are relying on INNER JOIN to chorizon table, as many mapunits have at least some components with no horizon data.

Data Filters

WHERE areasymbol != 'US'

/* Exclude Bedrock */
AND texturerv NOT IN ('BR', 'WB', 'UWB') 

/* Exclude "Duff" / Dry Organics */
AND ISNULL(horztype, '') != 'Dry Organic'
AND texturerv NOT IN ('SPM', 'MPM', 'HPM') 

AND desgnmaster != 'O'

q <- "
SELECT TOP 10
  component.compname,
  ch.hzname,
  ch.hzdept_r,
  ch.hzdepb_r,
  ch.texturerv,
  ch.horztype,
  ch.claytotal_r,
  ch.om_r
FROM chorizon AS ch
INNER JOIN component ON component.cokey = ch.cokey AND taxsubgrp LIKE '%Histic%'
ORDER BY component.cokey, hzdept_r
"

result <- SDA_query(q)
result

q <- "
SELECT TOP 10
  ch.hzname,
  ch.hzdept_r,
  ch.hzdepb_r,
  ch.texturerv,
  ch.horztype,
  ch.claytotal_r,
  ch.om_r
FROM chorizon AS ch
WHERE 
  -- 1. Exclude Bedrock
  ch.texturerv NOT IN ('BR', 'WB', 'UWB') 
  
  -- 2. Exclude Dry Organic / Duff
  AND ISNULL(ch.horztype, '') != 'Dry Organic'
  AND ch.texturerv NOT IN ('SPM', 'MPM', 'HPM')
ORDER BY cokey, hzdept_r
"

result <- SDA_query(q)
head(result)

q <- "
SELECT 
  COUNT(chkey) AS n, 
  kwfact, 
  kffact 
FROM chorizon 
WHERE desgnmaster = 'O' AND hzdept_r = 0
GROUP BY kwfact, kffact
ORDER BY n DESC
"
result <- SDA_query(q)
result$percent <- round(prop.table(result$n) * 100, 1)
head(result, 10)

Rock Fragment Data

Rock fragments are particles larger than 2 mm. In SSURGO, fragment data are stored in several places both aggregated up to the horizon level, and in child tables of horizon.

Fragments are classified by size class (e.g. gravel, cobbles, stones, boulders) and recorded as volume percentages. Derived weight fractions are also stored at the horizon level.

# Query horizons with their rock fragment content
q <- "
SELECT TOP 20
  mu.musym,
  co.compname,
  ch.hzname,
  ch.hzdept_r,
  ch.hzdepb_r,
  rf.fragkind,
  rf.fragsize_h,
  rf.fragvol_r
FROM legend AS leg
INNER JOIN mapunit AS mu ON mu.lkey = leg.lkey
INNER JOIN component AS co ON co.mukey = mu.mukey
LEFT JOIN chorizon AS ch ON ch.cokey = co.cokey
LEFT JOIN chfrags AS rf ON rf.chkey = ch.chkey
WHERE leg.areasymbol = 'CA630'
  AND co.majcompflag = 'Yes'
ORDER BY mu.musym, co.comppct_r DESC, ch.hzdept_r, rf.fragsize_h
"

result <- SDA_query(q)
head(result, 20)

# Total rock fragment content per horizon
q <- "
SELECT
  mu.musym,
  co.compname,
  co.comppct_r,
  ch.hzname,
  ch.hzdept_r,
  ch.hzdepb_r,
  SUM(rf.fragvol_r) AS calc_fragvoltot_r,
  COUNT(DISTINCT rf.fragsize_h) AS n_frag_classes
FROM legend AS leg
INNER JOIN mapunit AS mu ON mu.lkey = leg.lkey
INNER JOIN component AS co ON co.mukey = mu.mukey
LEFT JOIN chorizon AS ch ON ch.cokey = co.cokey
LEFT JOIN chfrags AS rf ON rf.chkey = ch.chkey
WHERE leg.areasymbol = 'CA630'
  AND co.majcompflag = 'Yes'
GROUP BY mu.musym, co.compname, co.comppct_r, ch.hzname,
         ch.hzdept_r, ch.hzdepb_r, ch.chkey
ORDER BY mu.musym, co.comppct_r DESC, ch.hzdept_r
"

result <- SDA_query(q)
head(result)

# Get rock fragment weight percentages (3-10 inch class)
frag_3_10 <- get_SDA_property(
  property = "Rock Fragments 3 - 10 inches - Rep Value",  # equivalent to 'frag3to10_r'
  method = "Dominant Component (Numeric)", # dominant component, weighted average by depth
  areasymbols = "CA630",
  top_depth = 0,
  bottom_depth = 100  # depth range for averaging
)

head(frag_3_10)

# get fragment volume (fragvoltot_l, fragvoltot_r, fragvoltot_h)
fragvol_total <- get_SDA_property(
  property = "fragvoltot_r",
  method = "Weighted Average", # weighted by component percentage and depth
  top_depth = 0,
  bottom_depth = 100,
  areasymbols = "CA630"
)

head(fragvol_total)

SDA Query Limitations

The SDA system has important hard limits:

1. Record Limit: Queries that return more than 100,000 records are rejected

2. Result Size: Response serialization limited to approximately 32 Mb

3. Memory Limit: Queries generally cannot exceed approximately 300 Mb memory usage

4. Timeouts and Complexity: Long-running queries may timeout or hit other internal limits during query planning

Note on Local vs. Remote: These limits primarily apply to the Remote SDA API. When querying a local SQLite database (e.g. created via createSSURGO() or downloadSSURGO()), you are generally unrestricted by record counts or timeouts. Complex tabular queries can often be run on entire states or the whole country locally without the need for chunking or simplification.

Strategies for Large Queries

1. Add WHERE clause to reduce result set:

This query will return >100k records (all components across all of SSURGO and STATSGO) (bad):

SELECT * FROM component

This query selects a specific component name pattern in WHERE clause (good):

SELECT * FROM component 
WHERE compname LIKE 'Miami%'

2. Use TOP clause to limit number of results (useful for testing):

Use a TOP clause to test with small result set first

SELECT TOP 100 * FROM chorizon

TOP can be used after sorting with ORDER BY to easily find the values associated with an ordered sequence. For instance, if you ORDER BY comppct_r DESC (order by component percentage in descending order), then take TOP 1 you will get the records associated with the dominant component. This is a very common aggregation scheme available as an option in most tools that use SSURGO data.

TOP is also generally useful for testing queries, or checking table result structure for a few resulting rows. If testing is successful, remove TOP, add filtering with WHERE or in JOIN condition (see below), or increase the TOP limit.

3. Aggregate and apply JOIN constraints early in the query:

Efficient queries often perform filtering operations as part of their JOIN statements. This can improve server-side query planning and reduce overhead that is not needed for the final result.

Return aggregated results (fewer rows) instead of raw data and constrain musym in join to mapunit instead of WHERE:

SELECT mu.mukey, 
  COUNT(*) AS n_hz,
  AVG(ch.claytotal_r) AS avg_clay
FROM component co
INNER JOIN mapunit mu ON mu.mukey = co.mukey AND
INNER JOIN chorizon ch ON ch.cokey = co.cokey
GROUP BY mu.mukey
WHERE compname LIKE 'Musick%'

4. Explicitly select the columns you need:

Efficient queries also are explicit about columns in the SELECT clause.

While SELECT * is fine to use for testing, it is prone to issues if schemas change, as well as returning more data than is necessary in most cases.

SELECT * will return all columns, usually including many columns you do not need:

SELECT TOP 1 *
FROM mapunit
INNER JOIN component ON mapunit.mukey = component.mukey
WHERE compname LIKE 'Musick%'

Use SELECT column1, column2, ... for concise results, with no duplicated columns in complex joins:

SELECT TOP 1 mapunit.mukey, musym, muname, 
             compname, localphase, comppct_r, majcompflag, compkind, cokey
FROM mapunit
INNER JOIN component ON mapunit.mukey = component.mukey
WHERE compname LIKE 'Musick%'

If a column name occurs multiple times in different tables, you need specify which one you want it to come from with <tablename>.<columname> syntax.

If you do need all columns from a specific table use <tablename>.* syntax, for example here we get all columns from chfrags, and specific columns from chorizon and component:

SELECT hzname, chfrags.*, component.cokey
FROM chfrags
INNER JOIN chorizon ON chorizon.chkey = chfrags.chkey
INNER JOIN component ON component.cokey = chorizon.cokey
WHERE compname LIKE 'Musick%'

5. Use makeChunks() for large vectors (see next section)

makeChunks() is a helper function in the `soilDB`` package. It takes a vector of values as input, and assigns each value a group (chunk) number. You can then calculate the unique chunk numbers and iterate over them to process data in bite-sized pieces. This is covered extensively in the next section.

T-SQL Syntax Specifics

SELECT TOP 10 * FROM chorizon WHERE claytotal_r IS NOT NULL

SELECT TOP 10 chkey, claytotal_r, ISNULL(claytotal_r, 0) AS clay 
FROM chorizon

SELECT TOP 10 
  claytotal_r,
  CAST(claytotal_r AS INT) AS clay_int,
  CAST(chkey AS CHAR(10)) AS chkey_str
FROM chorizon
WHERE claytotal_r != CAST(claytotal_r AS INT) 

SELECT * FROM mapunit WHERE muname LIKE '%clay%'

SELECT areasymbol + '-' + musym AS areamusym FROM mapunit

Why “Chunking”?

When you need data for many map units, components, soil survey areas, or other entities, a single query can exceed SDA’s limits.

SDA is a shared public resource, so it is generally encouraged to make requests that are as small and efficient as possible.

The makeChunks() function divides a large vector into smaller chunks for iterative querying.

Before implementing any chunked query approach, you should also be certain that your query can run properly on a small set, for example just single chunk, or a result that is truncated using a TOP clause.

Using makeChunks()

Suppose you have 5000 map unit keys to query. We can break that set up into chunks of 1000 map units, and iterate over each chunk using lapply()

large_mukey_list <- 461994:466995

# Divide into chunks of 1000
chunks <- soilDB::makeChunks(large_mukey_list, 1000)

# Now loop through chunks
results_list <- lapply(split(large_mukey_list, chunks), function(chunk_keys) {
  
  # Build IN clause
  in_clause <- soilDB::format_SQL_in_statement(chunk_keys) 
  
  # construct query 
  q <- sprintf(
    "SELECT mukey, musym, muname FROM mapunit WHERE mukey IN %s", in_clause
  )
    cat(sprintf("Chunk mukey %d to %d\n", min(chunk_keys), max(chunk_keys))) 
  SDA_query(q)
})

# Combine results
all_results <- do.call(rbind, results_list)

The above approach is readily converted to parallel processing, for example using future::future.lapply(), or other custom iterator functions that allow for progress reporting and other diagnostics.

If you do choose to use parallel processing, you should limit the number of concurrent connections made to the API to a reasonable number, say, 2 to 5 requests at a time, with a maximum of about 10. You should also implement proper error handling and retry mechanisms (see Handling Errors section above).

Example Chunked Property Query

Here we query soil properties for 6 soil survey areas in chunks of 2.

test_areasymbols <- c("CA067", "CA077", "CA632", "CA644", "CA630", "CA649")

# Chunk into groups of 5
chunks <- soilDB::makeChunks(test_areasymbols, 2)

results_list <- lapply(split(test_areasymbols, chunks), function(chunk_keys) {
  
  # Build IN clause
  in_clause <- soilDB::format_SQL_in_statement(chunk_keys) 
  
  q <- sprintf("
    SELECT
      mu.mukey,
      mu.musym,
      mu.muname,
      co.cokey,
      co.compname,
      SUM(hzdepb_r - hzdept_r) AS sum_hz_thickness
    FROM legend AS leg
    INNER JOIN mapunit AS mu ON mu.lkey = leg.lkey
    INNER JOIN component AS co ON co.mukey = mu.mukey
    INNER JOIN chorizon AS ch ON ch.cokey = co.cokey
    WHERE leg.areasymbol IN %s
    GROUP BY mu.mukey, mu.musym, mu.muname, co.cokey, co.compname
  ", in_clause)
  
  cat(sprintf("Chunk %s completed\n", paste0(chunk_keys, collapse = ", ")))
  SDA_query(q)
})

all_results <- do.call(rbind, results_list)
head(all_results)

Note: The above sample query is not a good way to determine soil thickness, as it does not exclude thickness of non-soil layers (e.g. bedrock) that should not be counted as part of the soil depth. As an exercise, use the logic from the previous sections to improve this example to filter out bedrock.

Recommended Chunk Sizes

• Typical queries: 1000-10000 items per chunk (safe for most queries)

• Simple lookups: 10000-100000 items per chunk (less data per item)

• Complex joins: 500-1000 items per chunk (more processing per item)

Monitor your query results and adjust based on SDA response times.

You may get an error:

Invalid query: The query processor ran out of internal resources and could not produce a query plan. This is a rare event and only expected for extremely complex queries or queries that reference a very large number of tables or partitions. Please simplify the query. If you believe you have received this message in error, contact Customer Support Services for more information.

If you see this, first confirm that your query syntax is correct for the result you are trying to obtain, for instance by trying it on a minimal small set of data. Critically evaluate the need for joins, subqueries, and complex filtering logic. Identify any areas of repeated logic and refactor to simplify redundant steps. If your query works as expected on a small data set, and you have optimized it, then consider breaking your big query up into (more) chunks.

Converting sf Objects to WKT

Important: SDA requires spatial inputs (WKT) to be in WGS84 (EPSG:4326) geographic coordinates.

Passing projected coordinates (like UTM or Albers) will result in query failures or zero results.

SDA does store an alternate projected version of the geometry data (e.g. mupolygonproj, rather than mupolygongeo). The projection is a Web Mercator projection ("EPSG:3857") so it is best avoided in analysis, though it is an option for visualization web map applications that lack capability to project the data.

You can use sf to define a geometry, such as a bounding recatangle, and convert it to WKT (Well-Known Text) for spatial queries:

library(sf)

# Define a bounding box for a region of interest
# (xmin, ymin, xmax, ymax)
bbox <- c(-120.9, 37.7, -120.8, 37.8)
bbox_sf <- st_as_sfc(st_bbox(c(
  xmin = bbox[1],
  ymin = bbox[2],
  xmax = bbox[3],
  ymax = bbox[4]
), crs = 4326))

wkt <- st_as_text(bbox_sf)

wkt

High-level Functions

The soilDB package provides two higher-level functions (SDA_spatialQuery() and fetchSDA_spatial()) that make obtaining spatial data easier. Both of these functions internally use SDA_query() and are often a better option for most use cases compared to crafting a custom spatial query.

• SDA_spatialQuery(): spatial inputs (e.g. point location, AOI rectangle) => spatial or tabular outputs

• fetchSDA_spatial(): tabular inputs (e.g. map unit key, area symbol, ecological site ID) => spatial outputs

The functions take different inputs, and support various outputs including SSURGO and STATSGO2 map unit polygon geometry, special feature points/lines/polygons, soil survey area boundaries, Major Land Resource Area (MLRA) boundaries, or simple tabular information using spatial data as input. These functions will handle projection to "EPSG:4326" internally as long as the input data has the correct CRS metadata.

# Example: Retrieve map unit polygons that intersect a bounding box
# This handles the WKT conversion and query construction automatically
# geomIntersection = TRUE clips the polygons to the bounding box

polygons <- SDA_spatialQuery(bbox_sf, what = "mupolygon", geomIntersection = TRUE)
polygons

Note that bbox_sf could be any geometry. Here is is a rectangular polygon, but it could be a point, or collection of points, or lines, or similar.

SDA_spatialQuery() allows you to get tabular data based on a spatial geometry.

For instance, to simply identify what soil survey area or map unit a single point overlaps with, you can do something like the following.

# Example: get some legend and mapunit-level tabular data at a point
point_sf <- sf::st_as_sf(data.frame(wkt = "POINT (-120.85 37.75)"),
                         wkt = "wkt",
                         crs = "EPSG:4236")
ssa <- SDA_spatialQuery(
  point_sf,
  what = "areasymbol"
)
ssa

mu <- SDA_spatialQuery(
  point_sf,
  what = "mukey"
)
mu

Above we pass a WKT string (the result of sf::st_centroid(bbox_sf)) within a data.frame to sf::st_as_sf(). We specify the wkt argument to specify which column the WKT string is stored in to construct our spatial query object.

This section will be expanded in the future, but for now an example application of the two higher-level functions is available in the “Dominant Ecological Site” vignette.

Spatial Queries: The Two-Step Pattern

For optimal performance, separate your spatial and tabular queries. Fetch the list of intersecting map unit keys (mukey) first, then use those keys to query attribute tables.

This avoids joining heavy geometry columns with complex attribute tables, which can trigger the 32Mb serialization limit. Then we can use the WKT in a custom SDA query to get map units intersecting the bounding box:

# Step 1: Get the mukeys that intersect the bounding box
q <- sprintf("
  SELECT DISTINCT mu.mukey
  FROM mapunit AS mu
  INNER JOIN mupolygon AS mp ON mp.mukey = mu.mukey
  WHERE mp.mupolygongeo.STIntersects(geometry::STGeomFromText('%s', 4326)) = 1
", wkt)

spatial_result <- SDA_query(q)
head(spatial_result, 5)

This query uses special geospatial functions defined in SQL Server to perform the geometric operations (in this case “intersection”)

SDA Spatial Helper Functions

In addition to standard T-SQL spatial functions, SDA provides pre-defined Table-Valued Functions (TVFs) optimized for intersections.

• SDA_Get_Mukey_from_intersection_with_WktWgs84(wkt): Returns mukey values intersecting the input WKT.
• SDA_Get_MupolygonWktWgs84_from_Mukey(mukey): Returns WKT geometry for a specific mukey.

This makes the logic much simpler for the common operations of finding mukey given a spatial location or geometry:

# Input MUST be WKT in WGS84 (EPSG:4326)
q <- "SELECT * FROM SDA_Get_Mukey_from_intersection_with_WktWgs84('POINT(-120.9 37.7)')"
result <- SDA_query(q)

Also, you can do the inverse (map unit key to geometry):

# Step 2: Get Geometry from Map Unit Key (mukey)
# Useful for retrieving the polygon boundary for a specific map unit
target_mukey <- 461994

q <- sprintf("SELECT * FROM SDA_Get_MupolygonWktWgs84_from_Mukey('%s')", target_mukey)

# Result contains the mukey and the geometry in WKT format
geometry_df <- SDA_query(q)

head(geometry_df)

Using the sf (or terra) package, you can convert the resulting text WKT column to a spatial object:

if (requireNamespace("sf", quietly = TRUE) &&
   !is.null(geometry_df)) {
  # Parse WKT column (usually named 'mupolygongeo' in mupolygon table, but 'MupolygonWktWgs84' in TVF result)
  poly_sf <- sf::st_as_sfc(geometry_df$MupolygonWktWgs84, crs = 4326)
  
  plot(poly_sf, main = paste("Geometry for mukey=", target_mukey))
}

Note: The TVF function name implies the CRS is WGS84 (authority:code EPSG:4326, SRS ID 4326; or, more correctly, OGC:CRS84).

See the SDA Advanced Query Guide for more details on TVF functions and other high-level features built into SDA.

Spatial and Tabular Property Integration

We can use the result of our spatial filtering with property queries.

Here we extract the unique map unit keys and pass as mukeys argument to get_SDA_property()

# Step 2: Use the mukeys to fetch tabular data
# First, get mukeys in bounding box
spatial_mukeys <- unique(spatial_result$mukey)

# Then query properties for those mukeys
if (length(spatial_mukeys) > 0) {
  clay_in_bbox <- get_SDA_property(
    property = "Total Clay - Rep Value",
    method = "Weighted Average",
    mukeys = spatial_mukeys,
    top_depth = 0,
    bottom_depth = 50
  )
  
  head(clay_in_bbox)
}

This shows how to get the data, see the Local SSURGO vignette for more spatial visualization examples using sf.

Available Properties

SDA has 80+ horizon properties and 40+ component properties.

Here are some common examples.

Component properties:

• "Taxonomic Order", "Taxonomic Suborder", "Taxonomic Temperature Regime"
• "Drainage Class", "Hydrologic Group"
• "Corrosion of Steel", "Corrosion of Concrete"
• "Range Production - Favorable/Normal/Unfavorable Year"

Horizon properties:

• Water: "Available Water Capacity - Rep Value", "0.1 bar H2O - Rep Value", "0.33 bar H2O - Rep Value", "15 bar H2O - Rep Value"

• Bulk Density: "Bulk Density 0.1 bar H2O - Rep Value", "Bulk Density 0.33 bar H2O - Rep Value", "Bulk Density 15 bar H2O - Rep Value", "Bulk Density oven dry - Rep Value"

• Texture: "Total Sand - Rep Value", "Total Silt - Rep Value", "Total Clay - Rep Value", and sand/silt sub-fractions

• Chemistry: "pH 1:1 water - Rep Value", "Electrical Conductivity - Rep Value", "Cation Exchange Capacity - Rep Value", "Sodium Adsorption Ratio - Rep Value", "Calcium Carbonate - Rep Value", "Gypsum - Rep Value"

• Other: "Saturated Hydraulic Conductivity - Rep Value", "Organic Matter - Rep Value", "Free Iron - Rep Value", "Oxalate Iron - Rep Value"

Alternately, you can directly use physical column names from the SSURGO schema like "claytotal_r", "ph1to1h2o_r", "ksat_r", etc.

Method: Dominant Component (Category)

Dominant component (category) get the categorical value from the dominant soil component. This query is a distinct method as categorical data need to be handled differently from numeric data; the numeric method is described below.

# Get taxonomic classification from dominant component
tax_order <- get_SDA_property(
  property = "Taxonomic Suborder",
  method = "Dominant Component (Category)",
  areasymbols = "CA630"
)

head(tax_order)

Method: Dominant Component (Numeric)

Dominant component (numeric) retrieves properties from the dominant component, optionally averaged over a depth range (in centimeters):

# Get bulk density at 0-25 cm depth
bulk_density <- get_SDA_property(
  property = c("dbthirdbar_l", "dbthirdbar_r", "dbthirdbar_h"),
  method = "Dominant Component (Numeric)",
  areasymbols = "CA630",
  top_depth = 0,
  bottom_depth = 25
)

head(bulk_density)

Method: Weighted Average

The weighted average method retrieves component percentage and depth-weighted average across all components in a map unit. In this example we explicitly exclude minor components by changing the default include_minors argument.

# Get weighted average clay content (0-50 cm)
clay_avg <- get_SDA_property(
  property = "Total Clay - Rep Value",
  method = "Weighted Average",
  areasymbols = "CA630",
  top_depth = 0,
  bottom_depth = 50,
  include_minors = FALSE
)

head(clay_avg)

Method: Min/Max

Method `“Min/Max” is used to get minimum or maximum values across all components in a mapunit.

# Get maximum sand content in any component
sand_max <- get_SDA_property(
  property = "Total Sand - Rep Value",
  method = "Min/Max",
  areasymbols = "CA630",
  FUN = "MAX",
  top_depth = 0,
  bottom_depth = 50
)

head(sand_max)

Method: Dominant Condition

Dominant condition aggregates components with the same category value for the properties of interest, then picks the dominant percentage. This gives extra weight to categories where multiple soils in the same map unit have the same value.

# Get dominant drainage class condition
drain_dominant <- get_SDA_property(
  property = "Drainage Class",
  method = "Dominant Condition",
  areasymbols = "CA630"
)

head(drain_dominant)

Method: None

Method "none" returns one row per component or component horizon, depending on the property selected, with no map unit aggregation. This generally returns much more data, but allows for custom aggregation outside of the SQL query, and deeper inspection of the source data.

You cannot mix component-level and horizon-level properties with this method. Run it twice, once for horizon-level and once for component-level, then join the component data to the horizon data using the cokey column to get a single data frame.

Also, the top_depth and bottom_depth arguments are deliberately ignored with this method. Apply your own filtering to the result if you only need a specific depth range.

# Get all pH values by component and horizon
ph_all <- get_SDA_property(
  property = "pH 1:1 water - Rep Value",
  method = "None",
  areasymbols = "CA630"
)

head(ph_all)

Inspecting Generated Queries

We can use query_string = TRUE to see the SQL generated by get_SDA_property():

q <- get_SDA_property(
  property = "Taxonomic Suborder",
  method = "Dominant Component (Category)",
  areasymbols = "CA630",
  query_string = TRUE
)

# View first 500 characters of query
cat(substr(q, 1, 500), "...\n")

Available Interpretations

SDA contains 600+ interpretations organized by category:

Method: Dominant Component

Get the rating from the dominant soil component:

# Get forestry ratings for dominant component
for_ratings <- get_SDA_interpretation(
  rulename = c("FOR - Potential Seedling Mortality",
               "FOR - Road Suitability (Natural Surface)"),
  method = "Dominant Component",
  areasymbols = "CA630"
)

head(for_ratings)

Method: Dominant Condition

Aggregate similar ratings, then return the dominant:

# Get dominant engineering interpretation
eng_ratings <- get_SDA_interpretation(
  rulename = "ENG - Septic Tank Absorption Fields",
  method = "Dominant Condition",
  areasymbols = "CA649"
)

head(eng_ratings)

Method: Weighted Average

Component-percentage-weighted average of ratings:

# Get weighted average agricultural rating
agr_weighted <- get_SDA_interpretation(
  rulename = "AGR - Winter Wheat Yield (MT)",
  method = "Weighted Average",
  areasymbols = "MT041",
  include_minors = TRUE
)

head(agr_weighted)

Note: in this example, a regional (Montana-specific) interpretation is used, so we used a Montana areasymbol (MT041). Some interpretations are only exported for specific states, regions, or soil survey areas.

Method: None

Return one row per component with no map unit aggregation:

# Get all component-level ratings
all_ratings <- get_SDA_interpretation(
  rulename = "FOR - Mechanical Planting Suitability",
  method = "None",
  areasymbols = "CA630"
)

head(all_ratings)

Understanding Output Columns

Use wide_reason = TRUE to pivot reason columns into separate sub-rule columns:

# Get detailed ratings with reasons pivoted
detailed <- get_SDA_interpretation(
  rulename = "FOR - Mechanical Planting Suitability",
  method = "Dominant Component",
  areasymbols = "CA630",
  wide_reason = TRUE
)

head(detailed)

Method: MAPUNIT (Default)

Returns an overall classification for each map unit:

hydric_class <- get_SDA_hydric(
  areasymbols = c("CA077", "CA630"),
  method = "MAPUNIT"
)

head(hydric_class)

# Check unique ratings
unique(hydric_class$HYDRIC_RATING)

Possible classifications:

• "Nonhydric" - No hydric components

• "Hydric" - All major components are hydric

• "Predominantly Hydric" - Hydric components >= 50%

• "Partially Hydric" - One or more major components are hydric

• "Predominantly Nonhydric" - Hydric components < 50%

The output includes:

• hydric_majors - Percentage of hydric major components

• hydric_inclusions - Percentage of hydric minor/inclusion components

• total_comppct - Total component percentage (should sum to 100)

Method: DOMINANT COMPONENT

Get hydric status of the dominant component only:

hydric_dom <- get_SDA_hydric(
  areasymbols = "CA630",
  method = "DOMINANT COMPONENT"
)

head(hydric_dom)

Method: DOMINANT CONDITION

Aggregate by hydric condition (Yes/No), then pick the dominant:

hydric_cond <- get_SDA_hydric(
  mukeys = c(461994, 461995, 462205),
  method = "DOMINANT CONDITION"
)

head(hydric_cond)

Method: NONE

Return one row per component (component-level hydric status):

hydric_all <- get_SDA_hydric(
  areasymbols = "CA630",
  method = "NONE"
)

head(hydric_all)

Filtering Results

You can filter by area symbol, map unit key, or custom WHERE clause:

# Get aggregate attributes for specific mukeys
muagg_filtered <- get_SDA_muaggatt(
  mukeys = c(461994, 461995, 463264)
)

head(muagg_filtered)

Method: DOMINANT COMPONENT

Get parent material from the dominant component:

pm_dom <- get_SDA_pmgroupname(
  areasymbols = "CA630",
  method = "DOMINANT COMPONENT",
  simplify = FALSE
)

head(pm_dom)

Method: DOMINANT CONDITION

Aggregate parent materials by group, then return the dominant:

pm_cond <- get_SDA_pmgroupname(
  mukeys = c(461994, 461995, 462205),
  method = "DOMINANT CONDITION"
)

head(pm_cond)

Simplify Parameter

By default, simplify = TRUE groups parent materials into broader categories. Set simplify = FALSE to get detailed group names:

# Simplified (broader groups)
pm_simple <- get_SDA_pmgroupname(
  mukeys = c(461994, 461995),
  simplify = TRUE
)

head(pm_simple)

# Detailed (specific groups)
pm_detailed <- get_SDA_pmgroupname(
  mukeys = c(461994, 461995),
  simplify = FALSE
)

head(pm_detailed)

Method: NONE (Default)

Return component-level ecological classes without aggregation:

eco_none <- get_SDA_coecoclass(
  method = "None",
  areasymbols = "CA630"
)

head(eco_none)

The default behavior targets the NRCS Ecological Site database-related entries, but there are many other ecological/vegetation class types in SSURGO.

Columns include:

• ecoclassid - Ecological Class ID

• ecoclassname - Full ecological class name

• ecoclasstypename - Type (e.g., "NRCS Rangeland Site", "NRCS Forestland Site")

• ecoclassref - Reference document

Method: DOMINANT COMPONENT

Get ecological site from the dominant component only:

eco_dom <- get_SDA_coecoclass(
  method = "Dominant Component",
  areasymbols = "CA630"
)

head(eco_dom)

Method: DOMINANT CONDITION

Aggregate ecological classes by ID, then return the dominant:

eco_cond <- get_SDA_coecoclass(
  method = "Dominant Condition",
  mukeys = c(461994, 461995, 462205),
  ecoclasstypename = "NRCS Forestland Site"
)

head(eco_cond)

Method: ALL

Return all ecological sites per map unit in wide format:

eco_all <- get_SDA_coecoclass(
  method = "All",
  mukeys = c(461994, 461995),
  threshold = 5
)

head(eco_all)

Filtering by ecological class type:

# Get rangeland sites only
eco_range <- get_SDA_coecoclass(
  method = "Dominant Condition",
  areasymbols = "CA630",
  ecoclasstypename = "NRCS Rangeland Site"
)

head(eco_range)

Available Tables

SDA contains four cosurfmorph tables:

Example: 3D Geomorphic Classification

# Get geomorphic position by component name
geom_3d <- get_SDA_cosurfmorph(
  table = "cosurfmorphgc",
  areasymbols = "CA630"
)

head(geom_3d)

Output includes columns like:

• compname - Component name

• geomposmntn, geomposmntn_n - Mountain position (count)

• p_geomposmntn - Mountain position proportion

Example: Hillslope Position by Area

# Get hillslope position aggregated by area symbol
geom_hill <- get_SDA_cosurfmorph(
  table = "cosurfmorphhpp",
  by = "areasymbol",
  areasymbols = "CA630"
)

head(geom_hill)

Example: Surface Shape

# Get surface shape classes
geom_shape <- get_SDA_cosurfmorph(
  table = "cosurfmorphss",
  areasymbols = "CA649"
)

head(geom_shape)

Example: Microrelief

# Get microrelief by component name
geom_micro <- get_SDA_cosurfmorph(
  table = "cosurfmorphmr",
  areasymbols = "CA630"
)

head(geom_micro)

Basic Usage

You can fetch all components and horizons for a specific area or set of map units using the WHERE argument:

library(aqp)

# Query soil components by areasymbol and musym
s <- fetchSDA(WHERE = "areasymbol = 'IN005' AND musym = 'MnpB2'")

# The result is a SoilProfileCollection
s

# Check the horizon data
head(horizons(s))

This function is very powerful for getting a quick snapshot of the soil profiles in an area.

It is worth noting that fetchSDA() returns data from both SSURGO and STATSGO2. To limit to SSURGO, use WHERE = "areasymbol != 'US'" (see the Data Filtering section), or specify a SSURGO area symbol.

Handling Duplicates

The duplicates argument (defaulting to FALSE) controls how components that appear in multiple survey areas are handled.

• duplicates = FALSE (Default): Returns only one instance of a component per nationalmusym. This eliminates “duplication” where multiple survey area (legend) map units are linked to the same source mapunit and component records. This is generally preferred for statistical analysis of soil properties.
• duplicates = TRUE: Returns a record for each unique map unit key (mukey). This is useful for visualization or when you need to account for every occurrence of a component across different survey area correlations.

Basic Usage

You can fetch lab data by specifying an identifier (like an area symbol) and the column it corresponds to (what argument):

# Fetch KSSL data for a specific soil survey area (CA630)
# 'what' argument specifies we are searching by 'area_code'
# 'tables' argument specifies which data tables to include (defaults to core tables)
ldm_data <- fetchLDM("CA630", what = "area_code")

# The result is a SoilProfileCollection
ldm_data

# Inspect site data
head(site(ldm_data))

Querying by Taxon Name

You can also search by correlated or sampled-as taxon names using a custom WHERE clause:

# Fetch lab data where the correlated or sampled name is 'Musick' or 'Holland'
# CASE statement handles differences between correlated and sampled names
ldm_taxon <- fetchLDM(WHERE = "(CASE WHEN corr_name IS NOT NULL 
                                THEN LOWER(corr_name) 
                                ELSE LOWER(samp_name) 
                            END) IN ('musick', 'holland')")

ldm_taxon

Advanced Usage: Specific Tables

By default, fetchLDM() retrieves a standard set of tables. You can request specific data, such as physical properties, using the tables argument:

# Fetch physical properties for soils correlated as "Typic Argialbolls"
ldm_phys <- fetchLDM(x = "Typic Argialbolls", 
                     what = "corr_taxsubgrp", 
                     tables = "lab_physical_properties")

# Inspect the available horizon data columns
names(horizons(ldm_phys))

fetchLDM supports retrieving data from various KSSL tables, including chemical properties, mineralogy, and x-ray analysis.

Example 1: Multi-Method Property Comparison

Compare the same property using different aggregation methods:

# Get clay content using three different methods
methods <- c("Dominant Component (Numeric)", "Weighted Average", "Max")

clay_results <- data.frame()

for (method in methods) {
  result <- get_SDA_property(
    property = "Total Clay - Rep Value",
    method = method,
    areasymbols = "CA630",
    top_depth = 0,
    bottom_depth = 50
  )
  
  result$method <- method
  clay_results <- rbind(clay_results, result)
}

# compare methods for a single map unit
subset(clay_results, mukey == 1865918)

Example 2: Properties and Interpretations for Land Use Planning

Combine soil properties and interpretations for a land use suitability assessment:

# Get drainage class and hydrologic group
drain_hydro <- get_SDA_property(
  property = c("Drainage Class", "Hydrologic Group"),
  method = "Dominant Condition",
  areasymbols = "CA630"
)

# Get an engineering interpretation
eng_interp <- get_SDA_interpretation(
  rulename = "ENG - Septic Tank Absorption Fields",
  method = "Dominant Condition",
  areasymbols = "CA630"
)

# Explicitly merge on mukey (and other shared columns to avoid duplication)
suitability <- merge(drain_hydro, eng_interp, 
                     by = c("mukey", "areasymbol", "musym", "muname"))

head(suitability)

Example 3: Hydric Status and Interpretations

Combine hydric classifications with wetland-related interpretations:

# Get a generalized mapunit-level hydric classification
# see ?get_SDA_hydric for details on method="mapunit"
hydric_stat <- get_SDA_hydric(
  areasymbols = "CA077",
  method = "MAPUNIT"
)

wet_interp <- get_SDA_interpretation(
  rulename = "WMS - Excavated Ponds (Aquifer-fed)",
  method = "Dominant Condition",
  areasymbols = "CA077"
)

wetland_assess <- merge(hydric_stat, wet_interp,
                        by = c("mukey", "areasymbol", "musym", "muname"),
                        all.x = TRUE)

subset(wetland_assess, rating_WMSExcavatedPondsAquiferfed < 0.6)

Example 4: Aggregating Across Multiple Areas

Here we use base R aggregate() (feel free to swap in your favorite R aggregation method e.g. dplyr, data.table, or collapse) to summarize across survey areas:

# Get properties for two areas
props_ca630 <- get_SDA_property(
  property = "Total Clay - Rep Value",
  method = "Dominant Component (Numeric)",
  areasymbols = "CA630"
)

props_ca649 <- get_SDA_property(
  property = "Total Clay - Rep Value",
  method = "Dominant Component (Numeric)",
  areasymbols = "CA649"
)

# Combine
all_props <- rbind(props_ca630, props_ca649)

# Aggregate by area symbol
area_summary <- aggregate(
  claytotal_r ~ areasymbol,
  data = all_props,
  FUN = function(x) {
    c(
      mean = mean(x, na.rm = TRUE),
      median = median(x, na.rm = TRUE),
      sd = sd(x, na.rm = TRUE),
      n = sum(!is.na(x))
    )
  }
)

area_summary

Example 5: Component-Level Analysis

Use method = "None" to perform custom component-level analysis:

# Get all component properties without aggregation
clay_by_comp <- get_SDA_property(
  property = "Total Clay - Rep Value",
  method = "None",
  areasymbols = "CA630",
  top_depth = 0,
  bottom_depth = 25,
  include_minors = TRUE
)

head(clay_by_comp)

# Calculate average clay by major vs. minor components
clay_summary <- aggregate(
  claytotal_r ~ majcompflag,
  data = clay_by_comp,
  FUN = mean
)

clay_summary

Key Takeaways

1. Use SDA_query() for custom SQL queries when predefined functions don’t fit your needs
2. T-SQL syntax: Use TOP, ISNULL(), CAST() for SDA queries
3. Query limits: Keep result sets < 100k records and < 32 Mb; use chunking for large datasets
4. Aggregation methods: Choose an appropriate method (Dominant Component, Weighted Average, None, etc.) for your analysis
5. Local vs. remote: Use dsn parameter to query local SSURGO databases (see Local SSURGO vignette)
6. Spatial queries: Combine WKT geometry with SDA_query() for location-based searches (or use SDA_spatialQuery())

Recommended Workflow

1. Start with small test queries using e.g. TOP 10, or a small number of inputs, to verify syntax and result table structure
2. Use query_string = TRUE to inspect generated SQL
3. Scale up with areasymbols or mukeys parameters
4. Use chunking for “large” data sets (greater than ~10,000 items, see makeChunks() section for suggested chunk sizes)
5. Cache results locally if repeating queries
6. Refer to official SDA documentation for table/column details

SDA Resources

• SDA Table Relationships Diagram: https://sdmdataaccess.sc.egov.usda.gov/documents/TableRelationshipsDiagram.pdf
• Column Descriptions: https://www.nrcs.usda.gov/sites/default/files/2022-08/SSURGO-Metadata-Table-Column-Descriptions-Report.pdf
• SDA Query Help
• SDA Query Samples
• SDA Table/Column Descriptions
• soilDB Reference Index
• Local SSURGO Vignette - For creating and querying local SSURGO databases
• Dominant Ecological Sites Vignette - For a fully worked example of combining spatial data with tabular (Ecological Site) data to make thematic maps

Core Tables

Common Columns