WCS GetCoverage for Grid returning CSV Example

From Earth Science Information Partners (ESIP)

Simple WCS Python Example

This example uses popular python scripting language to call Datafed WCS via SOAP.

Getting python

Windows users: go to http://www.activestate.com and download Python 2.4 or later. You can also go straight to download download python

The main site for python is http://www.python.org but I prefer to install via the activestate package. Especially Mark Hammonds PythonWin is a nice IDE.

Getting Demo Script

Download http://datafed.net/demo/soapdemo.txt and save it with extension .py

type from command line python soapdemo.py and watch it run.

There are no parameters, the script is hard coded for demo purposes.

How it works

After some declarations there is the WCS query wrapped in SOAP envelope. The triple quotes allow multiline strig constants in python. We have now variable example_soap_in


#   this WCS query queries USA temperatures 2005-12-06T12:00:00 at 1350 m height.
#   The Envelope and Body nodes are standard SOAP wrapping.
#   The actual query is the http-post query defined by WCS schema.

example_soap_in = """
<soap:Envelope xmlns:soap="http://schemas.xmlsoap.org/soap/envelope/">
    <soap:Body>
        <wcs:GetCoverage version="1.0.0" service="WCS"
                xmlns:gml="http://www.opengis.net/gml" xmlns:wcs="http://www.opengis.net/wcs">
            <!-- datafed dataset_abbr.param_abbr defines the coverage name -->
            <wcs:sourceCoverage>THREDDS_GFS.T</wcs:sourceCoverage> 
            <wcs:domainSubset>
                <wcs:spatialSubset>
                    <!--
                        This element queries the USA at 1350 m elevation.
                        If the dataset has no elevation, only lat and lon limits
                        are needed. 
                    -->
                    <gml:Envelope srsName="WGS84(DD)">
                        <!-- lon_min lat_min elev_min -->
                        <gml:pos>-130 24 1350</gml:pos>
                        <!-- lon_max lat_max elev_max -->
                        <gml:pos>-65 50 1350</gml:pos>
                    </gml:Envelope>
                    <gml:Grid dimension="3">
                        <gml:limits>
                            <!--
                                grid size. Datafed does not support resizing, so
                                these numbers have no meaning for grid datasets.
                            -->
                            <gml:GridEnvelope>
                                <gml:low>0 0 0</gml:low>
                                <gml:high>99 99 0</gml:high>
                            </gml:GridEnvelope>
                        </gml:limits>
                        <gml:axisName>lat</gml:axisName>
                        <gml:axisName>lon</gml:axisName>
                        <gml:axisName>elev</gml:axisName>
                    </gml:Grid>
                </wcs:spatialSubset>
                <wcs:temporalSubset>
                    <!--
                        query data for one time only.
                        Datafed currently supports also time ranges,
                        but only for stored time interval.
                        Enumerated times are not yet supported.                        
                    -->
                    <gml:timePosition>2005-12-06T12:00:00</gml:timePosition>
                </wcs:temporalSubset>
            </wcs:domainSubset>
            <wcs:output>
                <!--
                    For test purposes, use Comma Separated Values.
                    NetCDF and GML are also supported.
                -->
                <wcs:format>CSV</wcs:format>
            </wcs:output>
        </wcs:GetCoverage>
    </soap:Body>
</soap:Envelope>
"""


# This function performs WCS SOAP query to datafed
# INPUT: the soap query, similar to example_soap_in
# output: the data stream

def query_datafed_wcs(soap_in_text):
    # open connection to std port 80
    c = httplib.HTTPConnection("webapps.datafed.net", 80)
    c.connect()
    try:
        # the Web Service Definition Language definition is at
        # http://webapps.datafed.net/WCS.asmx?WSDL
        c.putrequest("POST", "/WCS.asmx")
        c.putheader("soapAction", "http://datafed.net/WCS/GetCoverage")
        c.putheader("content-type", "text/xml")
        c.putheader("content-length", repr(len(soap_in_text)))
        c.endheaders()
        c.send(soap_in_text)
        r = c.getresponse()
        # uncomment this if you want to see the status
        print 'response status ' + repr(r.status) + ' ' + repr(r.reason)

        # 200 means OK. Anything other is a failure.
        if r.status <> 200:
            raise 'http-error', repr(r.status) + ' ' + repr(r.reason)

        x = xml.dom.minidom.parse(r)
        # uncomment this if you want to see the output envelope
        print x.toxml()

        # We have the output. It does not contain data, just
        # metadata and an uri pointer to the real data.
        try:
            # The uri is written into 'data' node, so let's look for it.
            for u in x.getElementsByTagNameNS("http://datafed.net/xs/wcs", "data"):
                # uncomment next if you want to see the uri
                print 'fetching data ' + u.firstChild.nodeValue
                return urllib2.urlopen(urllib2.Request(u.firstChild.nodeValue))
            # We did not find the uri, so the result is an error            
            raise 'soap-error', x
        finally:
            # due to circular references, you have to have the next line.
            # Otherwise memory will leak. Other python implementations
            # like Jython or IronPython rely on platform provided
            # garbage collection, and don't need it. 
            x.unlink()
    finally:
        #closing expensive resources ASAP is a good practice
        c.close()

#
#   If you run this script from command line, the next will be executed
#
if __name__ == "__main__":
    print 'querying datafed'
    csv = query_datafed_wcs(example_soap_in)
    print "here's the data"
    print ''
    # we have the csv stream, so let's look at it
    print csv.read()