Is there a downloadable file of CONUS mapunit polygons?

[kevinwolz]

Currently I download mapunit polygons using SDA_spatialQuery(what = "mupolygon"), but I'm wondering if, similar to the downloadable file of the gSSURGO mapunit raster, there is a downloadable file of CONUS mapunit polygons.

[brownag]

@kevinwolz You can use the MUPOLYGON layer from either the ESRI File Geodatabase or the SSURGO Portal Geopackage for CONUS. e.g. st_read("your/path/ssurgoportal/database.gpkg", "MUPOLYGON")

For survey area boundaries see layer SAPOLYGON.

You can get the CONUS and state-level gSSURGO and SSURGO Portal data on Box. https://nrcs.app.box.com/s/w0gtf7ooxqfcd0cgfht5zprcmc8qvu0i & https://nrcs.app.box.com/v/soils/folder/17971946225.

At this time there is no standalone download for just the vector data analogous to the GeoTIFF MUKEY grid, you need to download the associated tabular info with the spatial. But I believe you have already downloaded the "all SSURGO" dataset, so it should be as simple as extracting and saving the mapunit polygons as a separate layer if needed.

I hope the new 2023 file geodatabases should be up soon fingers crossed. We have state-level gSSURGO but not CONUS or MUKEY grids.

[kevinwolz]

@brownag I have a quick follow up question for you on this.

Per your suggestion, I downloaded the 10/1/2024 ALL_SSURGO Geopackage and pulled out the mupolygon layer into its own GPKG file. However, even this file is 56GB! As such, it is impractical and unnecessary to read the whole file in. So, I'm using the wkt_filter argument of sf::st_read():

sf::st_read("mupolygon.gpkg", layer = "mupolygon", wkt_filter = sf::st_as_text(aoi))

This works great except for one thing: it returns much more data than the SDA alternative.

Here is my aoi in yellow:

Here are the polygons returned by soilDB::SDA_spatialQuery(aoi, what = "mupolygon", db = "SSURGO"). Exactly as expected - just the survey polygons that intersect aoi.

And here are the polygons returned by sf::st_read("mupolygon.gpkg", layer = "mupolygon", wkt_filter = sf::st_as_text(aoi)). It appears that the polygons in "mupolygon.gpkg" are grouped by survey area and so all polygons in the survey area are being returned.

Now I can easily process this data after it is read in to intersect it with aoi and clean things up, but this is still not ideal. The shear amount of data being returned when using the local file is enough to really slow things down. It's not the 56GB of the full CONUS dataset, but it can still be sizeable, especially for AOIs that happen to be in large survey areas (e.g. out West).

Do you have a suggestion for a better way to read in from the local version of the database in a way that returns something analagous to soilDB::SDA_spatialQuery(aoi, what = "mupolygon", db = "SSURGO")? Thanks!

[kevinwolz]

I'm guessing that the answer here is to utilize a more "direct" SQL interface with the GPKG file, but I'm not sure the best way to do that.

[brownag]

I can confirm that ALL_SSURGO_10_01_2024.gpkg has MULTIPOLYGON geometries that contain the full extent of each unique MUKEY, which I believe is causing this behavior that gives you these large SSA-bounded results.

This, I believe, is the result of the default behavior in the SSURGO Portal tool which "dissolves" polygons for MUKEY to create a multipolygon geometry:

Disable Disolve
By default, the polygons in the soilmu_a_XXXX.shp file are dissolved on the “MUKEY” column during import. Selecting this option removes the dissolving function.

Dissolving significantly reduces the number of records in the SQLite SSURGO Template database mupolygon table and improves performance of this layer in GIS software. Removing the dissolve will increase the database sizes. It will also make draw times, joins, relates, and selections take longer in GIS.

Dissolving does not result in any data loss.

Your use of the wkt_filter argument, or custom OGR query argument, should be fine when working with geopackages that were created with the "Disable dissolve" option checked on. Unfortunately that means you would need to download all the SSAs and build the geopackage yourself. I can check with our national database manager regarding the benchmarks they used for choosing this as the default behavior, and the issues it can pose.

I compared the results I get for a GPKG I created for a recent request I had for Humboldt Toiyabe National Forest, where I did disable dissolve because this is how I am used to seeing the data, and you get something much more like what you expect for a given AOI polygon:

library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE

dsn <- "C:/Geodata/ALL_SSURGO_10_01_2024/ALL_SSURGO_gpkg/ALL_SSURGO_10_01_2024.gpkg"
dsn2 <- "D:/RGISS/HumboldtToiyabeNF/HumboldtToiyabeNF_gpkg/HumboldtToiyabeNF.gpkg"
x2 <- sf::st_read(dsn2, query = "SELECT * FROM mupolygon LIMIT 1")
#> Reading query `SELECT * FROM mupolygon LIMIT 1'
#> from data source `D:\RGISS\HumboldtToiyabeNF\HumboldtToiyabeNF_gpkg\HumboldtToiyabeNF.gpkg' 
#>   using driver `GPKG'
#> Simple feature collection with 1 feature and 4 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -120.1404 ymin: 39.28037 xmax: -120.1353 ymax: 39.28624
#> Geodetic CRS:  WGS 84

# create an AOI polygon that is a bbox around a single POLYGON geometry
p <- sf::st_as_sf(wk::as_rct(st_bbox(x2[1,])))
# plot(st_geometry(p))

# MULTIPOLYGON, full extent of 6 mukey; pretty fast
res <- st_read(dsn, "mupolygon", wkt_filter = wk::as_wkt(p))
#> Reading layer `mupolygon' from data source 
#>   `C:\Geodata\ALL_SSURGO_10_01_2024\ALL_SSURGO_gpkg\ALL_SSURGO_10_01_2024.gpkg' 
#>   using driver `GPKG'
#> Simple feature collection with 6 features and 4 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -120.8304 ymin: 39.03236 xmax: -120.0053 ymax: 39.73071
#> Geodetic CRS:  WGS 84

plot(st_geometry(res))
plot(st_geometry(p), col = "RED", add = TRUE)

# MULTIPOLYGON, but only the portion of 6 mukeys
dsn2 <- "D:/RGISS/HumboldtToiyabeNF/HumboldtToiyabeNF_gpkg/HumboldtToiyabeNF.gpkg"
res2 <- sf::st_read(dsn2, "mupolygon", wkt_filter = wk::as_wkt(p))
#> Reading layer `mupolygon' from data source 
#>   `D:\RGISS\HumboldtToiyabeNF\HumboldtToiyabeNF_gpkg\HumboldtToiyabeNF.gpkg' 
#>   using driver `GPKG'
#> Simple feature collection with 7 features and 4 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -120.1969 ymin: 39.24355 xmax: -120.1072 ymax: 39.31529
#> Geodetic CRS:  WGS 84

plot(st_geometry(res2))
plot(st_geometry(p), col = "RED", add = TRUE)

Note also that st_read() has a couple default behaviors that can lead to MULTIPOLYGON results that is: type and promote_to_multi arguments. While it might be tempting to set type=3 to force POLYGON geometry, I would need to test that a bit more comprehensively to make sure that is a good idea for SSURGO in general such that it doesn't lead to omissions of more complicated geometries that may exist. In general SSURGO mupolygon shouldn't have MULTIPOLYGON geometries, but there are lots of complex geometries with "holes" and similar phenomena, and possibly some actual MULTIPOLYGON geometries in the source that escaped QA checks (not sure on this last point). Per the st_read() manual page "[with type=3] secondary rings from MULTIPOLYGONS may be dropped without warnings" so it is worth being a little careful about this approach.

[dylanbeaudette]

Yes. I'd suggest using ST_Intersects(geom, AOI) in your query to the local file. Depending on the size and frequency of queries, you are going to exceed the limits of file-based storage.

Have a look at the latest update to the SSURGO Portal product, there is an "all SSURGO" option.

https://www.nrcs.usda.gov/resources/data-and-reports/ssurgo-portal

You will of course run into similar problems querying with the wkt_filter argument to st_read(). I've never tested, but I'd wager (a small amount) that this approach to loading small AOIs will not effectively use the spatial index.

As a point of reference, a properly indexed copy of SSURGO spatial + tabular is on the order of 70GB.

Suggestions:

• use SDA for small AOIs
• customize the GeoPKG with additional spatial indexing, if missing
• query the geometry via SQL interface, using a combination of && (overlap) and ST_Intersection()
• ensure that the SQL interface is using WKB vs. WKT, i.e. prevent unnecessary conversion to/from
• use a RDBMS to store the data

[kevinwolz]

Thanks, @dylanbeaudette. For this application I cannot use SDA - I need an offline option for various reason, including intermittent downtime of the SDA API and an AOI size limit through that interface.

When using something this...

  SELECT * 
  FROM mupolygon 
  WHERE ST_Intersects(geometry, ST_GeomFromText('{aoi_wkt}', 4326))

I still get the whole, intersected survey area because (I think) the individual soil POLYGONs are grouped into MULTIPOLYGONs at the soil survey level.

I could ST_Cast the layer to POLYGON before doing the ST_Intersects, but that would cast all of CONUS and probably take days. Perhaps I first crop to the bounding box of the AOI, then cast, and finally to the intersection?

I'm not familiar with "using a combination of && (overlap) and ST_Intersection()".

[dylanbeaudette]

That is strange, how did you make or obtain the mupolygon table? I've never encountered this particular combination of factors.

I don't use GeoPKG all that often, is there an ST_GeometryType() function? In PostGIS, I'd use:

SELECT ST_GeometryType(geometry) FROM mupolygon LIMIT 1;

I just tried this in a SQLite DB and it didn't work. Maybe the client wasn't enabling the appropriate modules.

I agree that CASTing anything in this database is a bad idea. Long-term, I'd recommend building or downloading a version of these data that do not store the mupolygon geometry as a MULTI- type.

[kevinwolz]

Long-term, I'd recommend building or downloading a version of these data that do not store the mupolygon geometry as a MULTI- type.

Hah, yes that would be ideal. I just want USDA to do that for me. ;.) I'm not sure I have enough memory on any machine I have access to. An efficient, mukey vector polygon version analogous to the GeoTIFF MUKEY grid would be dreamy.

To make the mupolygon.gpkg, I first unzipped the 7z file on a big Linux machine:

sudo apt update
sudo apt install p7zip-full

7z x ALL_SSURGO_10_01_2024.7z

Unzipping the 7z file produces a file called ALL_SSURGO_10_01_2024.gpkg. To extract just the mupolygon layer from this, I did:

sudo apt update
sudo apt install gdal-bin
ogr2ogr -f "GPKG" mupolygon.gpkg ALL_SSURGO_10_01_2024.gpkg  mupolygon

I just assumed GPKG was the safest route because that's what was provided by USDA for the whole database.

[dylanbeaudette]

Interesting approach. The data are massive and best imported in-pieces, usually via iteration over SSA. This is how I typically load the individual SSURGO archives into an external database.

What does ogrinfo have to say about the original file? Are those MULTI geometries? If so, then we've got some splainin' to do!

If the original file does not encode as MULTI, then I'd suggest querying that directly.

You are probably one of the first people to download and attempt anything with one of these behemoth, all-in-on SSURGO products. Sorry for the inconvenience.

Efficient storage of a vector database this large is really only possible via something like PostGIS.

[kevinwolz]

What does ogrinfo have to say about the original file? Are those MULTI geometries? If so, then we've got some splainin' to do!

Yes, ogrinfo says the following about the mupolygon layer in the original ALL_SSURGO_10_01_2024.gpkg. They are Multipolygons from the beginning.

Layer name: mupolygon
Geometry: Multi Polygon
Feature Count: 326993
Extent: (-179.213580, -14.373855) - (179.859650, 71.441059)
Layer SRS WKT:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]
Data axis to CRS axis mapping: 2,1
FID Column = objectid
Geometry Column = shape
areasymbol: String (0.0)
spatialver: Integer64 (0.0)
musym: String (0.0)
mukey: Integer64 (0.0)

You are probably one of the first people to download and attempt anything with one of these behemoth, all-in-on SSURGO products. Sorry for the inconvenience.

That actually makes me feel better! I always feel like I'm being a pain with these questions because all of the SQL stuff is on the edge of my competence.

I definitely understand, though, that I am trying to do a very specific thing with what is ultimately not the best dataset for it: I am trying to create reports about individual farms (small AOIs) anywhere in CONUS (but must have easy, instant access to all of CONUS) that look good (at small-scale, that means vector data) and doesn't get tripped up when the SDA server is randomly out of service (happens regularly). All those factors lead me to a local source of mukey vector polygons. If any one factor were different, then I would switch to either SDA polygons or local mukey rasters.

[dylanbeaudette]

Helpful, thanks for checking. The MULTI-feature encoding was news to me, and I just confirmed with a colleague that it was intentional. I'll follow-up as to the particular reasons.

In the meantime, you might consider performing the intersection:

SELECT mukey, ST_Intersection(geom, AOI) as geom
FROM
mupolygon
WHERE ST_Intersects(geom, AOI);

What does that return? There are some additional optimizations that are worth trying--I use these on a regular basis in PostGIS, unsure if they work in SQLite.

I completely understand the decision to work with a local copy. SDA down-time is a problem that regularly creates headaches for me when testing things late at night. I'll talk with our SDA devs about that.

[dylanbeaudette]

I did some follow-up testing of single vs. multi-polygons in SSURGO. In general, expect intersection-based queries and subsets to be at least 10x slower with MULTIPOLYGON vs. POLYGON. File sizes of single vs. multi- geometry are almost the same. Seems the main rationale for using a default MULTIPOLYGON representation was render speed in Arc Pro.

I'm working with the SSURGO Portal devs on more flexible options and better documentation. Due to the nature of BBOX-based spatial indexing, I can't see any reason to prefer MULTIPOLYGON representation when using the database for analytical purposes.