Reading OSM Data Using Python

OpenStreetMap is a free editable community driven map.  The data can be downloaded in full from planet.openstreetmap.org, but you can also pull smaller sections from various web services.

For world pulls, I wouldn't recommend this code because you will run out of memory.  For smaller areas, same between 1-10 MB this should work depending on your desktop's specifications.

import xml
from xml.dom.minidom import parse
 
data = {}
keys = []
 
ways = {}
wkeys = []

relations = {}
rkeys = []

xmls = [r"c:\temp\small.osm"]
for f in xmls:
    dom = xml.dom.minidom.parse(f)
    #  parse out osm nodes
    #
    for n in dom.getElementsByTagName("node"):
        nid = n.getAttribute("id")
        data[nid] = {}
        data[nid]["lat"] = n.getAttribute("lat")
        data[nid]["lon"] = n.getAttribute("lon")
        for tag in n.getElementsByTagName("tag"):
            if(tag.hasAttribute("k")):
                k = tag.getAttribute("k")
                if(k not in keys):
                    keys.append(k)
                if(tag.hasAttribute("v")):
                    data[nid][k] = tag.getAttribute("v")
    # parse out osm ways/polygons
    #
    for n in dom.getElementsByTagName("way"):
        wid = n.getAttribute("id")
        ways[wid] = {}
        ways[wid]['ref'] = []
        ways[wid]['geomType'] = ""
        for nd in n.getElementsByTagName('nd'):
            if nd.hasAttribute('ref'):
                ref = nd.getAttribute('ref')
                ways[wid]['ref'].append(ref)
            del nd
        for tag in n.getElementsByTagName("tag"):
            if tag.hasAttribute("k") and \
               tag.hasAttribute('v'):
                k = tag.getAttribute("k")
                if k not in wkeys:
                    wkeys.append(k)
                ways[wid][k] = tag.getAttribute('v')
            del tag
        first = ways[wid]['ref'][0]
        last = ways[wid]['ref'][len(ways[wid]['ref'])-1]
        if ways[wid]['ref'][0] == ways[wid]['ref'][len(ways[wid]['ref'])-1]:
            ways[wid]['geomType'] = "polygon"
        else:
            ways[wid]['geomType'] = "polyline"
    for n in dom.getElementsByTagName('relation'):
        rid = n.getAttribute('id')
        relations[rid] = {}
        relations[rid]['member'] = []
        for mem in n.getElementsByTagName('member'):
            member = {}
            if mem.hasAttribtue('type'):
                member['type'] = mem.getAttribute('type')
            if mem.hasAttribute('ref'):
                member['ref'] = mem.getAttribute('ref')
            if mem.hasAttribute('role'):
                member['role'] = mem.getAttribute('role')
            relations[rid]['member'].append(member)
            del mem
        for tag in mem.getElementsByTagName("tag"):
            if tag.hasAttribute("k") and \
               tag.hasAttribute('v'):
                k = tag.getAttribute("k")
                relations[rid][k] = tag.getAttribute('v')            
            del tag
# do stuff with nodes, ways (lines), and polygons like converting the data from dictionaries to feature classes or shapefiles!

OSM files are just XML files that can be parsed with various xml processing tools in python. Some common ones are the minidom (core to python), lxml, and beautiful soup.

It is not the most efficient piece of code, but it demonstrates the power of python.

Enjoy, and feel free to make tweaks!

ArcGIS 10.2 SP 1 Announcement

Just in case you missed this with all the holiday madness, here is the 10.2.1 announcement.

http://support.esri.com/en/downloads/patches-servicepacks/view/productid/160/metaid/2055

Support Python GUI Support in ArcGIS

Please support my idea on having Python GUI support in ArcGIS.

http://ideas.arcgis.com/ideaView?id=087E00000004SmHIAU

I want to get this idea to over 1000 points before the beginning of 2014, and I need your help!

Thank you

GPX to Feature Class

A handy tool created at 10.1, but refined at 10.2 is the GPX to Feature Class tool, which allows users to convert GPX files from GPS units to feature classes.  Waypoints and Tracks can then be extracted from the converted GPX file, and rendered as desired.

The ArcGIS online help for the GPX to Feature Class can be found here.

import arcpy
import os
if __name__ == "__main__":
    input_gpx = r"C:\data.gpx"
    convert_fc = arcpy.env.scratchGDB + os.sep + "gpxfc"
    gpx_fc = arcpy.GPXtoFeatures_conversion(Input_GPX_File=input_gpx, 
                                            Output_Feature_class=convert_fc
                                           )[0]
    track_points_lyr = arcpy.MakeFeatureLayer_management(in_features=convert_fc, 
                                                        out_layer="TRKLYR", 
                                                        where_clause="Type = 'TRKPT'", 
                                                        )[0]
    way_point_lyr = arcpy.MakeFeatureLayer_management(in_features=convert_fc, 
                                                      out_layer="WPT", 
                                                        where_clause="Type = 'WPT'", 
                                                        )[0]    

    if int(arcpy.GetCount_management(track_points_lyr)[0]) > 0:
        trkpt = arcpy.CopyFeatures_management(way_point_lyr,
                                              env.scratchGDB + os.sep + "trkpt")[0]
        trkLine = arcpy.PointsToLine_management(trkpt,
                                                env.scratchGDB + os.sep + "lines",
                                                "Name")[0]
    if int(arcpy.GetCount_management(way_point_lyr)[0]) > 0:
        waypts = arcpy.CopyFeatures_management(way_point_lyr,
                                               env.scratchGDB + os.sep + "WAYPTS")[0]
        

Here the code converted the GPX file to a feature class, then further refined the results from just a set of points to tracks and waypoints as mentioned earlier in this post.

There is a difference between 10.1 and 10.2 that should be noted.  In 10.2, more of the GPX properties will be converted where as in 10.1, I have noticed that only a small set of the data is actually pulled in from the GPX files.  So if you need all the GPX properties, ArcGIS 10.2 is probably what you need.  You can always parse the XML in the GPX manually using the minidom library in python if this doesn't suit you needs.


Enjoy

Geometry Objects Make Life Easier

Sometimes you just want to work with a geometry, but you do not always want to go through all the steps of creating a cursor object.  In the 10.x framework, you can output almost any geoprocessing tool to a arcpy.Geometry object.

To gain access to all feature's geometries in a feature class, just do the following:

import arcpy
from arcpy import env
fc = r"c:\temp\data.shp"
geoms = arcpy.CopyFeatures_management(fc, arcpy.Geometry())
for g in geoms:
   print g.extent

This sample allows anyone to directly access the geometries of the input feature class without having to use the Cursor objects.

The same idea can be applied to other functions to the analysis function as well:

import arcpy
from arcpy import env
fc = r"c:\temp\data.shp"
geom = arcpy.Buffer_analysis(fc, arcpy.Geometry(), "100 Feet", "FULL", "ROUND")[0]
print geom.extent

Here the buffer tool outputs a single geometry to the geom object and the extent is displayed.

Where this becomes really powerful is when you need to perform geometry operations on your data, and want to put the results back into that row.

import arcpy
from arcpy import env
fc = r"c:\temp\data.shp"
with arcpy.da.UpdateCursor(fc, ["SHAPE@"]) as urows:
   for urow in urows:
      geom = arcpy.Buffer_analysis(urow[0], arcpy.Geometry(), "100 Feet", "FULL", "ROUND")[0]
      row[0] = geom
      urows.updateRow(urow)
      del urow
      del geom

Assuming that the input is a polygon, this snippet shows how geometries can be used as inputs and outputs thus allowing for easy insertion back into the original row.

Hope this helps!

Happy Halloween!

Happy Halloween.  Everyone Have Fun And Be Safe!

Table to Excel (10.2 Arcpy)

Table to excel is a great tool for converting your spatial data into an excel spreadsheet.  In the past, before 10.2, you would have to use a 3rd party module like Python Excel in order to convert your data to and from ArcGIS tables to excel spreadsheets.

It should be noted that this only support the MS Excel 5.0/95 format (.xls) file type.

The syntax and use is pretty simple for this tool, it has 4 inputs, the table path (string), out path (string), use field alias (Boolean and optional), and use domain and sub-type description (Boolean and optional).  Let's dive into an example:

import arcpy
import os
if __name__ == '__main__':
    table = r"c:\temp\scratch.gdb\StarbucksAddresses"
    out_excel = r"c:\temp\starbucks.xls"
    if os.path.isfile(out_excel):
        os.remove(out_excel)
    arcpy.TableToExcel_conversion(table, out_excel)

All we have done here is taken the simplest example and convert my list of Starbucks locations and converted to an excel spreadsheet so I can now use it in excel.

Enjoy

ArcGIS Tool help link - http://resources.arcgis.com/en/help/main/10.2/index.html#/Table_To_Excel/001200000054000000/