Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I am trying to write a python script which pulls GML tags out of an XML file and formats them into WKT for inserting into a PostGIS database. I have been successful in doing so for an XML containing a single part polygon using the follow code:

    rootElement = ET.parse("GMLExample_Polygon.xml").getroot()
    wkt = ""
    for subelement in rootElement.getiterator():
        for subsub in subelement:
            if subsub.tag == "{http://www.opengis.net/gml}X":
                x = subsub.text
            if subsub.tag == "{http://www.opengis.net/gml}Y":
                y = subsub.text
                point_for_pol = "%s %s, " % (x, y)
                wkt += point_for_pol
    wkt = wkt[:-2]

This code clearly won't work for multipolygons. I am unsure how to access the geometry for each polygon tag ("gml:Polygon srsName="BNG") separately and pull only the geometry nested under it. I am trying to use ElementTree, not sure if this is the best module to use? The XML is structured as follows:

<Order xsi:noNamespaceSchemaLocation="http://lalala.com/xml_polygon_order.xsd">
    <OrderRequest>
        <CustomerReference>999998</CustomerReference>
        <SiteAddress>
            <Premise>456</Premise>
            <Street>long street</Street>
            <Locality/>
            <Town>London</Town>
            <County>London</County>
            <PostCode>PN1 1PN</PostCode>
        </SiteAddress>
        <SiteGeography>
            <gml:Polygon srsName="BNG">
                <gml:outerBoundaryIs>
                    <gml:LinearRing>
                        <gml:coord>
                            <gml:X>452847.6009</gml:X>
                            <gml:Y>18596.0496</gml:Y>
                        </gml:coord>
                        <gml:coord>
                            <gml:X>415847.6009</gml:X>
                            <gml:Y>184596.0496</gml:Y>
                        </gml:coord>
                        <gml:coord>
                            <gml:X>415847.6009</gml:X>
                            <gml:Y>184596.0496</gml:Y>
                        </gml:coord>
                    </gml:LinearRing>
                </gml:outerBoundaryIs>
            </gml:Polygon>
            <gml:Polygon srsName="BNG">
                <gml:outerBoundaryIs>
                    <gml:LinearRing>
                        <gml:coord>
                            <gml:X>452847.6009</gml:X>
                            <gml:Y>18596.0496</gml:Y>
                        </gml:coord>
                        <gml:coord>
                            <gml:X>415847.6009</gml:X>
                            <gml:Y>184596.0496</gml:Y>
                        </gml:coord>
                        <gml:coord>
                            <gml:X>415847.6009</gml:X>
                            <gml:Y>184596.0496</gml:Y>
                        </gml:coord>
                    </gml:LinearRing>
                </gml:outerBoundaryIs>
            </gml:Polygon>
        </SiteGeography>
    </OrderRequest>
</Order>

Thanks for any help.

share|improve this question
1  
Hi Matt - I have never done it myself but a company called Astun Technology released an open source GML loader for Ordnance Survey GML files. You can find it here github.com/AstunTechnology/Loader. In the code it uses a GFS file that states what the structure of the GML looks like so the python module LXML can understand it. There may be some resuable code in there for your situation. –  tjmgis Apr 17 '13 at 15:35

2 Answers 2

up vote 4 down vote accepted

I enjoy using ElementTree. It's standardized in Python since 2.5 as xml.etree.ElementTree. Forgive me for being blunt, but you're using it wrong. I suggest trying the find, findtext, and findall methods when you know the structure of the data. Is Order your root element? If so,

>>> geography = rootElement.find('OrderRequest/SiteGeography')
>>> for polygon in geography.findall('{http://www.opengis.net/gml}Polygon'):
...     for coord in polygon.findall(
...             "{http://www.opengis.net/gml}outerBoundaryIs/"
...             "{http://www.opengis.net/gml}LinearRing/"
...             "{http://www.opengis.net/gml}coord"):
...         print(
...             coord.findtext("{http://www.opengis.net/gml}X"),
...             coord.findtext("{http://www.opengis.net/gml}Y"))
... 
('452847.6009', '18596.0496')
('415847.6009', '184596.0496')
('415847.6009', '184596.0496')
('452847.6009', '18596.0496')
('415847.6009', '184596.0496')
('415847.6009', '184596.0496')

http://pymotw.com/2/xml/etree/ElementTree/parse.html#finding-nodes-in-a-document has more advice on using ElementTree.

share|improve this answer

And in addition to the reply of sgillies, if you want the wkt format (wkt = "" ?), use his Shapely module:

from shapely.geometry.polygon import LinearRing

linearing = []
for polygon in geography.findall('{http://www.opengis.net/gml}Polygon'):
     for coord in polygon.findall("{http://www.opengis.net/gml}outerBoundaryIs/{http://www.opengis.net/gml}LinearRing/{http://www.opengis.net/gml}coord"):
         linearing.append((float(coord.findtext("{http://www.opengis.net/gml}X")),float(coord.findtext("{http://www.opengis.net/gml}Y"))))

print LinearRing(a).wkt
'LINEARRING (452847.6009000000194646 18596.0495999999984633, 415847.6009000000194646 184596.0495999999984633, 415847.6009000000194646 184596.0495999999984633, 452847.6009000000194646 18596.0495999999984633, 415847.6009000000194646 184596.0495999999984633, 415847.6009000000194646 184596.0495999999984633, 452847.6009000000194646 18596.0495999999984633)'

or with the two outerBoundaryIs:

for outerBoundaryIs in geography.findall('{http://www.opengis.net/gml}Polygon/{http://www.opengis.net/gml}outerBoundaryIs'):
    linearing = []
    for coord in outerBoundaryIs.findall("{http://www.opengis.net/gml}LinearRing/{http://www.opengis.net/gml}coord"):
        linearing.append((float(coord.findtext("{http://www.opengis.net/gml}X")),float(coord.findtext("{http://www.opengis.net/gml}Y"))))         
    print LinearRing(linearing).wkt

'LINEARRING (452847.6009000000194646 18596.0495999999984633, 415847.6009000000194646 184596.0495999999984633, 415847.6009000000194646 184596.0495999999984633, 452847.6009000000194646 18596.0495999999984633)'
'LINEARRING (452847.6009000000194646 18596.0495999999984633, 415847.6009000000194646 184596.0495999999984633, 415847.6009000000194646 184596.0495999999984633, 452847.6009000000194646 18596.0495999999984633)'

and if you want the srsName of the polygon:

for polygon in geography.findall('{http://www.opengis.net/gml}Polygon'):
    polygon.attrib.get('srsName') 
share|improve this answer
3  
And if you want to avoid Shapely's GEOS dependency, there's github.com/larsbutler/geomet, an emerging WKB/WKT serializer. –  sgillies Apr 17 '13 at 21:26

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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