5

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.

Improve this question
1
  • 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
    Commented Apr 17, 2013 at 15:35

2 Answers 2

Reset to default
8

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.

Improve this answer
3

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')
Improve this answer
2
  • 3
    And if you want to avoid Shapely's GEOS dependency, there's github.com/larsbutler/geomet, an emerging WKB/WKT serializer.
    – sgillies
    Commented Apr 17, 2013 at 21:26
  • print LinearRing(a).wkt will probably fail, and not just because of the Python2 print.
    – Eric Duminil
    Commented Mar 18, 2023 at 21:27

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.