SoilPoint · blog soilpoint.dev
← SoilPoint

How to Get USDA Soil Data (SSURGO) from a Lat/Lon in Python, Without the T-SQL Pain

By Adrian · June 2026, updated September 2026 · ~8 min read

View as Markdown Open in Claude

Last updated: September 2026. I'm a solo developer, not a soil scientist. I built the JSON API at the end of this post because I got tired of the workflow in the first half. Everything before it is the free way to do it yourself, and it works.

You have a latitude and longitude. You want the soil at that point: the dominant soil series, texture, drainage class, hydrologic group, pH, available water capacity (AWC), organic matter. In the US, the authoritative source is SSURGO, the USDA-NRCS Soil Survey Geographic Database, and getting it programmatically is harder than it has any right to be.

This post walks the whole path:

  1. What SSURGO actually is, and why a single point is a multi-table join
  2. The official way: POST T-SQL to USDA's Soil Data Access (SDA), with a runnable Python example
  3. Why that's only step one: you still have to write the depth-weighted aggregation yourself
  4. The gotchas (query hangs, the SSURGO+STATSGO mix, row/byte limits, the JSON+COLUMNNAME quirk)
  5. A one-line JSON alternative
  6. Limits, and when you should use the free R/Python packages instead

If you only want a working answer right now, jump to “The one-line version”. Otherwise, let’s do it the real way first so the shortcut makes sense.

What SSURGO actually is

SSURGO isn’t a grid of values you can sample. It’s a vector polygon survey with a normalized relational schema bolted onto it. To get from a coordinate to a number, you traverse a hierarchy:

spatial polygon (a "map unit" on the ground)
   └─ mapunit        ← identified by mukey; e.g. "Diablo clay, 15 to 30 percent slopes"
        └─ component ← the soil series that make up that map unit, each with a % (comppct_r)
             └─ chorizon ← the soil horizons (layers) within a component, with depths
                  └─ pH, AWC, organic matter, sand/silt/clay... measured per horizon

So “the pH at this point” is really:

Most tutorials stop before that last step. More on it below.

SSURGO values come as _l / _r / _h (low / representative / high). Use _r (representative) unless you specifically need the range.

The official way: POST T-SQL to Soil Data Access

USDA exposes SSURGO through Soil Data Access (SDA), a web service that takes T-SQL in the body of an HTTP POST and returns the result set. The REST endpoint is:

https://SDMDataAccess.sc.egov.usda.gov/Tabular/post.rest

The docs are at sdmdataaccess.nrcs.usda.gov/WebServiceHelp.aspx. It’s authoritative, it’s free, and there’s no API key. Credit where due. It’s also an ASP.NET help page with no quickstart and no lat/lon recipe, which is why you’re reading this.

Step 1: coordinate → mukey

SDA ships a spatial helper, SDA_Get_Mukey_from_intersection_with_WktWgs84, that takes a WKT point in EPSG:4326. Note the WKT order is POINT(lon lat), longitude first.

import requests

SDA_URL = "https://SDMDataAccess.sc.egov.usda.gov/Tabular/post.rest"

def sda_query(sql: str) -> list[dict]:
    """POST a T-SQL query to SDA and return rows as dicts."""
    resp = requests.post(
        SDA_URL,
        json={"query": sql, "format": "JSON+COLUMNNAME"},
        timeout=120,  # SDA can hang ~a minute; don't set this low
    )
    resp.raise_for_status()
    table = resp.json().get("Table", [])
    if not table:
        return []
    header, *rows = table  # JSON+COLUMNNAME puts column names in row 0
    return [dict(zip(header, r)) for r in rows]

lat, lon = 37.368402, -121.771
mukey_sql = f"""
SELECT mukey
FROM SDA_Get_Mukey_from_intersection_with_WktWgs84('point({lon} {lat})')
"""
mukey = sda_query(mukey_sql)[0]["mukey"]
print(mukey)  # -> 1882921

That JSON+COLUMNNAME format is an SDA quirk: instead of an array of objects, you get a Table array whose first row is the column names and the rest are values. The helper above re-zips it back into dicts.

The most-cited SDA point query stops here, with a mukey and a map-unit name. No texture, no drainage, no pH. To get usable numbers you keep going.

Step 2: mukey → dominant component → horizon properties

Now the joins. You need to:

  1. Find the dominant component for the mukey (MAX(comppct_r)).
  2. Pull its horizons from chorizon, each with a top/bottom depth (hzdept_r, hzdepb_r).
  3. Depth-weight each property over your target window (here, 0–30 cm).

The depth-weighting is most of the work. For a property like pH over 0–30 cm, you weight each horizon’s value by how much of that horizon falls inside the window:

weighted_value = Σ (horizon_value × overlap_thickness) / Σ (overlap_thickness)

where overlap_thickness = max(0, min(horizon_bottom, 30) − max(horizon_top, 0)).

You can push some of this into SQL, but it gets ugly fast, and SDA’s T-SQL doesn’t make windowed depth math pleasant. A pragmatic split is: get the dominant component and its horizons from SDA, then do the weighting in Python where it’s readable and testable.

comp_sql = f"""
SELECT TOP 1 co.cokey, co.compname, co.comppct_r,
       co.drainagecl, co.hydgrp
FROM component AS co
WHERE co.mukey = '{mukey}'
ORDER BY co.comppct_r DESC
"""
comp = sda_query(comp_sql)[0]

hz_sql = f"""
SELECT ch.hzdept_r, ch.hzdepb_r,
       ch.ph1to1h2o_r AS ph,
       ch.awc_r       AS awc,
       ch.om_r        AS om
FROM chorizon AS ch
WHERE ch.cokey = '{comp["cokey"]}'
ORDER BY ch.hzdept_r
"""
horizons = sda_query(hz_sql)

def depth_weight(horizons, field, top=0, bottom=30):
    num = den = 0.0
    for h in horizons:
        val = h.get(field)
        if val in (None, ""):
            continue
        ht, hb = float(h["hzdept_r"]), float(h["hzdepb_r"])
        overlap = max(0.0, min(hb, bottom) - max(ht, top))
        if overlap <= 0:
            continue
        num += float(val) * overlap
        den += overlap
    return round(num / den, 2) if den else None

result = {
    "series":        comp["compname"],
    "componentPct":  comp["comppct_r"],
    "drainageClass": comp["drainagecl"],
    "hydrologicGroup": comp["hydgrp"],
    "ph_0_30":  depth_weight(horizons, "ph"),
    "awc_0_30": depth_weight(horizons, "awc"),
    "om_0_30":  depth_weight(horizons, "om"),
}
print(result)

This works. But notice what you just built: a brittle multi-query pipeline with hand-rolled depth math, no caching, and a hard dependency on a service that likes to hang.

The gotchas

If you go the SDA route in production, budget time for these:

None of this is a knock on USDA. SDA gives away an enormous, valuable dataset for free. It’s built for SQL-literate batch analysts, not for “give me JSON for this coordinate in a request handler.”

The clunky-but-good alternatives

Before the shortcut, here are the other options:

If one of those fits your workflow, use it. The rest of this post is for the case where you just want hosted JSON per coordinate, no install, no SQL, and cached answers that keep working through SDA outages.

The one-line version

That’s the API I built. It does the point-in-polygon lookup, picks the dominant component, depth-weights the 0–30 cm properties, and returns clean JSON. It keeps every answer in memory, so a point that’s been looked up since the last deploy still answers when SDA is down.

curl "https://soilpoint.dev/v1/soil?lat=37.368402&lon=-121.771"

Real response for Diablo clay near San Jose, CA (September 2026):

{
  "lat": 37.368402,
  "lon": -121.771,
  "mukey": "1882921",
  "mapUnitName": "Diablo clay, 15 to 30 percent slopes, MLRA 15",
  "soilSeries": "Diablo",
  "componentKind": "Series",
  "componentPct": 85,
  "texture": "Clay",
  "drainageClass": "Well drained",
  "hydrologicGroup": "C",
  "awc": 0.17,
  "ph": 6.8,
  "organicMatter": 2.5,
  "nccpi": 0.387,
  "nccpiClass": "Moderately low inherent productivity",
  "depthWindowCm": [0, 30],
  "dataVintage": "2025-08-29",
  "stale": false,
  "cachedAt": "2026-09-24T03:55:08.802Z",
  "dataSource": "USDA-NRCS SSURGO via Soil Data Access",
  "surveyDisclaimer": "Provided as-is; not a substitute for on-site investigation; not for regulatory determinations; US coverage only. Not a USDA-endorsed product."
}

Same fields, one GET, no T-SQL, no key. Here’s the Python equivalent of the whole first half of this article:

import requests

def soil_at(lat, lon):
    r = requests.get("https://soilpoint.dev/v1/soil",
                     params={"lat": lat, "lon": lon}, timeout=30)
    r.raise_for_status()
    return r.json()

print(soil_at(42.0308, -93.6319))
# Iowa: Spillville loam, drainage "Moderately well drained",
# hydrologic group "B/D", pH 6.5, AWC 0.2, organic matter 4.5%

How it works (the short version)

Here’s what happens behind that GET. It’s the same pipeline from the first half:

Limits

I’d rather you hit these in this post than in production:

So which should you use?

curl "https://soilpoint.dev/v1/soil?lat=YOUR_LAT&lon=YOUR_LON"

If you try it and something’s wrong or missing, tell me. I’m a solo dev and I read everything, and knowing where it breaks for you is how it gets better.

SoilPoint wraps USDA-NRCS Soil Data Access. It is not a USDA-endorsed product. Soil data is provided as-is and is not a substitute for on-site investigation.

Want to try it with your own coordinates?

Try it live API docs