1

I'm using this python code to extract data from the SSURGO database. At the moment I'm doing for a single lat long point. Would it be possible to run the query for a list of coordinates extracted from a csv file and save the results of each coordinate into a separate data frame?

import requests
import json
import xmltodict
import pprint
import pandas as pd

lat = "37.54189"
lon = "-120.96683"
lonLat = lon + " " + lat
url="https://SDMDataAccess.nrcs.usda.gov/Tabular/SDMTabularService.asmx"
#headers = {'content-type': 'application/soap+xml'}
headers = {'content-type': 'text/xml'}
body = """<?xml version="1.0" encoding="utf-8"?>
          <soap:Envelope xmlns:soap="http://www.w3.org/2003/05/soap-envelope" xmlns:sdm="http://SDMDataAccess.nrcs.usda.gov/Tabular/SDMTabularService.asmx">
   <soap:Header/>
   <soap:Body>
      <sdm:RunQuery>
         <sdm:Query>SELECT co.cokey as cokey, ch.chkey as chkey, comppct_r as prcent, slope_r, slope_h as slope, hzname, hzdept_r as deptht, hzdepb_r as depthb, awc_r as awc, 
                    claytotal_r as clay, silttotal_r as silt, sandtotal_r as sand, om_r as OM, dbthirdbar_r as bulk_density, wthirdbar_r as th33, ph1to1h2o_r as pH, ksat_r as sat_hidric_cond,
                    (dbthirdbar_r-wthirdbar_r)/100 as bd FROM sacatalog sc
                    FULL OUTER JOIN legend lg  ON sc.areasymbol=lg.areasymbol
                    FULL OUTER JOIN mapunit mu ON lg.lkey=mu.lkey
                    FULL OUTER JOIN component co ON mu.mukey=co.mukey
                    FULL OUTER JOIN chorizon ch ON co.cokey=ch.cokey
                    FULL OUTER JOIN chtexturegrp ctg ON ch.chkey=ctg.chkey
                    FULL OUTER JOIN chtexture ct ON ctg.chtgkey=ct.chtgkey
                    FULL OUTER JOIN copmgrp pmg ON co.cokey=pmg.cokey
                    FULL OUTER JOIN corestrictions rt ON co.cokey=rt.cokey
                    WHERE mu.mukey IN (SELECT * from SDA_Get_Mukey_from_intersection_with_WktWgs84('point(""" + lonLat + """)')) order by co.cokey, ch.chkey, prcent, deptht
        </sdm:Query>
      </sdm:RunQuery>
   </soap:Body>
</soap:Envelope>"""

response = requests.post(url,data=body,headers=headers)
print("###############################################\n")
# Put query results in dictionary format
my_dict = xmltodict.parse(response.content)

# Convert from dictionary to dataframe format
soil_df = pd.DataFrame.from_dict(my_dict['soap:Envelope']['soap:Body']['RunQueryResponse']['RunQueryResult']['diffgr:diffgram']['NewDataSet']['Table'])

1 Answer 1

2

Create a function which take coordinates as input and return a dataframe:

import requests
import xmltodict
import pandas as pd

def soil(lat, lon):
    #lat = "37.54189"
    #lon = "-120.96683"
    lonLat = "{0} {1}".format(lon, lat)
    url="https://SDMDataAccess.nrcs.usda.gov/Tabular/SDMTabularService.asmx"
    #headers = {'content-type': 'application/soap+xml'}
    headers = {'content-type': 'text/xml'}
    body = """<?xml version="1.0" encoding="utf-8"?>
              <soap:Envelope xmlns:soap="http://www.w3.org/2003/05/soap-envelope" xmlns:sdm="http://SDMDataAccess.nrcs.usda.gov/Tabular/SDMTabularService.asmx">
       <soap:Header/>
       <soap:Body>
          <sdm:RunQuery>
             <sdm:Query>SELECT co.cokey as cokey, ch.chkey as chkey, comppct_r as prcent, slope_r, slope_h as slope, hzname, hzdept_r as deptht, hzdepb_r as depthb, awc_r as awc, 
                        claytotal_r as clay, silttotal_r as silt, sandtotal_r as sand, om_r as OM, dbthirdbar_r as bulk_density, wthirdbar_r as th33, ph1to1h2o_r as pH, ksat_r as sat_hidric_cond,
                        (dbthirdbar_r-wthirdbar_r)/100 as bd FROM sacatalog sc
                        FULL OUTER JOIN legend lg  ON sc.areasymbol=lg.areasymbol
                        FULL OUTER JOIN mapunit mu ON lg.lkey=mu.lkey
                        FULL OUTER JOIN component co ON mu.mukey=co.mukey
                        FULL OUTER JOIN chorizon ch ON co.cokey=ch.cokey
                        FULL OUTER JOIN chtexturegrp ctg ON ch.chkey=ctg.chkey
                        FULL OUTER JOIN chtexture ct ON ctg.chtgkey=ct.chtgkey
                        FULL OUTER JOIN copmgrp pmg ON co.cokey=pmg.cokey
                        FULL OUTER JOIN corestrictions rt ON co.cokey=rt.cokey
                        WHERE mu.mukey IN (SELECT * from SDA_Get_Mukey_from_intersection_with_WktWgs84('point(""" + lonLat + """)')) order by co.cokey, ch.chkey, prcent, deptht
            </sdm:Query>
          </sdm:RunQuery>
       </soap:Body>
    </soap:Envelope>"""
    
    response = requests.post(url,data=body,headers=headers)
    print("###############################################\n")
    # Put query results in dictionary format
    my_dict = xmltodict.parse(response.content)
    
    # Convert from dictionary to dataframe format
    
    soil_df = pd.DataFrame.from_dict(my_dict['soap:Envelope']['soap:Body']['RunQueryResponse']['RunQueryResult']['diffgr:diffgram']['NewDataSet']['Table'])
    return soil_df

frames = []
coords = [[37.54189, -120.96683], [37.9014, -121.1895]]
for coord in coords: #For each coordinate pair, create a df and append to frames list
    frames.append(soil(*coord))
1
  • 1
    Thank you so much for your answer, I will try it out and report back. Many Many thanks!
    – Ralt_Jol
    Commented Jan 30, 2022 at 15:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.