1

I am parsing an OpenStreetMap gpx file which I exported from here. I'm trying to put the data into arcgis geodatabase feature classes. Some of the features in the gpx are points (saved as wpt), while others are polygons (saved as trk consisting of trkpt's). Using the script below, I am succesful in copying over the points, their geometry, the names and attributes I want. However somehow creating the geometry does not work for the polygons, although the attributes are copied over without trouble (I left the bit of code for the attributes out of the example below). I found this question which prompted me to try to build the polygons with the points ordered the other way round. However this didn't work either. I built in quite some print statements to check where it goes wrong, see printout below. Also not that printing poly.area results in 0.0, printing poly.trueCentroid results in None, etc. So the polygon simply isn't generated from the array. To me it appears that poly = arcpy.Polygon(ar) simply does not work.

GPX example:

<wpt lat="55.7962684" lon="-4.6260225">
 <name>Ian Brown Motors</name> 
 <desc>name=Ian Brown Motors shop=car_repair</desc> 
 <link href="http://osm.org/browse/node/3452158061" /> 
</wpt>
<trk>
 <name>Key Store Brodie's Minimarket</name> 
 <desc>building=yes name=Key Store Brodie's Minimarket opening_hours=08:00-18:00 shop=supermarket source=survey</desc> 
 <link href="http://osm.org/browse/way/28476052" /> 
 <trkseg>
  <trkpt lat="56.4387108" lon="-4.7134965" /> 
  <trkpt lat="56.4385955" lon="-4.7135666" /> 
  <trkpt lat="56.438564" lon="-4.7133971" /> 
  <trkpt lat="56.4386793" lon="-4.7133269" /> 
  <trkpt lat="56.4387108" lon="-4.7134965" /> 
 </trkseg>
</trk>

Script:

from xml.etree import ElementTree as Et
import os
import arcpy

arcpy.env.overwriteOutput = 1

wgs = arcpy.SpatialReference(r'C:\Program Files\ArcGIS\Desktop10.0\Coordinate '
                             r'Systems\Geographic Coordinate Systems\World\WGS '
                             r'1984.prj')

gdb = PATH_TO_GDB
f = PATH_TO_GPX
all_records = []
p = Et.parse(f)
root = p.getroot()

arcpy.CreateFeatureclass_management(gdb, 'shops_pt', 'POINT',
                                    spatial_reference=wgs)
arcpy.CreateFeatureclass_management(gdb, 'shops_pl', 'POLYGON',
                                    spatial_reference=wgs)

shops_pt = os.path.join(gdb, 'shops_pt')
shops_pl = os.path.join(gdb, 'shops_pl')

for fc in [shops_pt, shops_pl]:
    arcpy.AddField_management(fc, 'name', 'TEXT', 255)
    arcpy.AddField_management(fc, 'type', 'TEXT', 50)
    arcpy.AddField_management(fc, 'opening_hours', 'TEXT', 255)

ic_pt = arcpy.InsertCursor(shops_pt)
ic_pl = arcpy.InsertCursor(shops_pl)

for mainElement in root:
    if mainElement.tag == 'metadata':
        continue
    if mainElement.tag == 'wpt':
        r = ic_pt.newRow()
        y = float(mainElement.attrib['lat'])
        x = float(mainElement.attrib['lon'])
        p = arcpy.Point(x, y)
        r.shape = p
        ic_pt.insertRow(r)

    elif mainElement.tag == 'trk':
        r = ic_pl.newRow()
        for subElement in mainElement:
            if subElement.tag == 'trkseg':
                ar = arcpy.Array()
                for trkpt in subElement:
                    y = float(trkpt.attrib['lat'])
                    x = float(trkpt.attrib['lon'])
                    p = arcpy.Point(x, y)
                    print p
                    ar.add(p)
                print 'array length', len(ar)

                print 'ar[0]', ar.getObject(0), 'ar[-1]', ar.getObject(len(ar)-1)
                poly = arcpy.Polygon(ar)
                r.shape = poly
                print 'poly.centroid try 1', poly.centroid
                if poly.centroid is None:
                    ar2 = arcpy.Array()
                    for p in ar:
                        ar2.insert(0, p)
                    for p in ar2:
                        print p
                    poly = arcpy.Polygon(ar2)
                    print 'poly.centroid try 2', poly.centroid
        ic_pl.insertRow(r)

Printout sample for 1 polygon:

-4.7134965 56.4387108 NaN NaN
-4.7135666 56.4385955 NaN NaN
-4.7133971 56.438564 NaN NaN
-4.7133269 56.4386793 NaN NaN
-4.7134965 56.4387108 NaN NaN
array length 5
ar[0] -4.7134965 56.4387108 NaN NaN ar[-1] -4.7134965 56.4387108 NaN NaN
poly.centroid try 1 None
-4.7134965 56.4387108 NaN NaN
-4.7133269 56.4386793 NaN NaN
-4.7133971 56.438564 NaN NaN
-4.7135666 56.4385955 NaN NaN
-4.7134965 56.4387108 NaN NaN
poly.centroid try 2 None

Solution: as suggested by @Vince below in the comments, when putting the Array in a Polygon, explicitly declare the spatial reference, so it would become poly = arcpy.Polygon(ar, wgs).


8
  • 1
    Try applying the spatial reference to the shape -- poly = arcpy.Polygon(ar,wgs)
    – Vince
    Jun 25, 2015 at 12:08
  • Yes, I would do as Vince says and add a spatial reference, otherwise the polygon has no idea where it is supposed to be. Also, you can set up the wgs by passing in a well-known ID (WKID) rather than referencing a .prj file: wgs = arcpy.SpatialReference(4326)
    – crmackey
    Jun 25, 2015 at 12:17
  • The default spatial reference likely has a vertex resolution of 0.001 (1mm), which collapses the vertices, making the centroid undefined.
    – Vince
    Jun 25, 2015 at 12:22
  • @Vince that works well, thanks a lot. Would you care writing it up in an answer?
    – Menno
    Jun 25, 2015 at 12:28
  • 1
    @crmackey it gives me an error if I use it like that, but if I do it in the way at the bottom of this page it works fine: help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//…
    – Menno
    Jun 25, 2015 at 12:29

1 Answer 1

1

I think the comments on this question are worth placing into an answer which I'll make Community Wiki:

Try applying the spatial reference to the shape -- poly = arcpy.Polygon(ar,wgs) – Vince Jun 25 at 12:08


Yes, I would do as Vince says and add a spatial reference, otherwise the polygon has no idea where it is supposed to be. Also, you can set up the wgs by passing in a well-known ID (WKID) rather than referencing a .prj file: wgs = arcpy.SpatialReference(4326) – crmackey Jun 25 at 12:17


The default spatial reference likely has a vertex resolution of 0.001 (1mm), which collapses the vertices, making the centroid undefined. – Vince Jun 25 at 12:22


@crmackey it gives me an error if I use it like that, but if I do it in the way at the bottom of this page it works fine: help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//… – Menno Jun 25 at 12:29


I should point out that my suggestion only works for version 10.1 and beyond. – crmackey Jun 25 at 12:37

Your Answer

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

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