SQL Fiddle - Your Database Scratchpad

Ever want to test something out on your database but did not want to mess everything up. Now is your chance, check out: sqlfiddle.  Just go to this site and select your database.  You can then create a quick schema and dummy data and bam! Test away.


Enjoy

Add Geometry Attributes - ArcGIS 10.2.1 Tool

New at 10.2.1 is a simplified ways of adding geometry properties to a feature class.  The Add Geometry Attributes geoprocessing tool is a short cut from using the calculate field's geometry calculations.

Let's use the old way first:

    fc = r"some datasource"
    arcpy.AddField_management(fc, fieldName1, "DOUBLE", 
                              fieldPrecision, fieldScale)
    arcpy.AddField_management(fc, fieldName2, "DOUBLE", 
                              fieldPrecision, fieldScale)
 
    # Calculate centroid
    arcpy.CalculateField_management(fc, fieldName1, 
                                    "!SHAPE.CENTROID.X!",
                                    "PYTHON_9.3")
    arcpy.CalculateField_management(fc, fieldName2, 
                                    "!SHAPE.CENTROID.Y!",
                                    "PYTHON_9.3")

So here we create the field and then perform a calculation on the data set to get the X/Y centroid of the data.

At 10.2.1 you can clean this up and just do the following:

fc = "some dataset"
arcpy.AddGeometryAttributes_management(fc, "CENTROID")

With just two lines of code, we did what took multiple line previously to 10.2.1.

The Add Geometry Attributes tool is a time saver, and makes life just a tad easier.  It's a convenient tool that allows for quick addition of geometry attributes.

Check out more about this tool here.

Enjoy

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!