Interesting Tim Peters interview...

Really enjoyed seeing the face behind the developer...

Tim Peters (core developer on CPython) discusses the early days of Python language development including meeting Guido, and porting Python to 64bit.  Additonal topics include:

• Why people are leaving compiler development.
• Constraints of free development time
• Sprint-based development
• Talks about new language features
• Parallelism and processes

Point Clustering w/ ArcGIS + Python

Here is a quick example of clustering an arcgis feature class of points.  WARNING: This is more of a "gist" than a production ready version.  I am releasing it to support my presentation at the Esri Developer Summit.  Come back soon and it should be updated with additional error handling, much much more efficient use of the tools.

import json
import multiprocessing
import os
import unittest

from collections import defaultdict

import arcpy

def cluster_worker(sourcePointFeatureClass, outputDirectory, levelInfos=None, whereClause=None, workerName=None, includeLevels=None):
    print '.rectify_worker({})'.format(workerName)
    arcpy.env.overwriteOutput = True
    outputGeodatabaseName = "temp{}.gdb".format(workerName)
    outputGeodatabase = os.path.join(outputDirectory, outputGeodatabaseName)
    if arcpy.Exists(outputGeodatabase):
        arcpy.Delete_management(outputGeodatabase)
    arcpy.CreateFileGDB_management(outputDirectory, outputGeodatabaseName)
    arcpy.env.workspace = outputGeodatabase

    points_feature_layer = 'temp_points'
    points_feature_class = 'source_points'
    
    arcpy.MakeFeatureLayer_management(sourcePointFeatureClass, points_feature_layer, whereClause)
    arcpy.CopyFeatures_management(points_feature_layer, points_feature_class)

    level_infos_json = """{"origin":{"x":-20037508.342787,"y":20037508.342787},"spatialReference":{"wkid":102100},"lods":[{"level":0,"resolution":156543.033928,"scale":591657527.591555},{"level":1,"resolution":78271.5169639999,"scale":295828763.795777},{"level":2,"resolution":39135.7584820001,"scale":147914381.897889},{"level":3,"resolution":19567.8792409999,"scale":73957190.948944},{"level":4,"resolution":9783.93962049996,"scale":36978595.474472},{"level":5,"resolution":4891.96981024998,"scale":18489297.737236},{"level":6,"resolution":2445.98490512499,"scale":9244648.868618},{"level":7,"resolution":1222.99245256249,"scale":4622324.434309},{"level":8,"resolution":611.49622628138,"scale":2311162.217155},{"level":9,"resolution":305.748113140558,"scale":1155581.108577},{"level":10,"resolution":152.874056570411,"scale":577790.554289},{"level":11,"resolution":76.4370282850732,"scale":288895.277144},{"level":12,"resolution":38.2185141425366,"scale":144447.638572},{"level":13,"resolution":19.1092570712683,"scale":72223.819286},{"level":14,"resolution":9.55462853563415,"scale":36111.909643},{"level":15,"resolution":4.77731426794937,"scale":18055.954822},{"level":16,"resolution":2.38865713397468,"scale":9027.977411},{"level":17,"resolution":1.19432856685505,"scale":4513.988705},{"level":18,"resolution":0.597164283559817,"scale":2256.994353},{"level":19,"resolution":0.298582141647617,"scale":1128.497176}],"units":"esriMeters"}"""
    level_infos = json.loads(unicode(level_infos_json, 'latin-1'))
    clustered_feature_classes = []
    for level_info in level_infos['lods']:
        if not includeLevels or int(level_info['level']) in includeLevels:
            cluster_feature_class = create_cluster_layer(points_feature_class, level_info)
            if cluster_feature_class:
                clustered_feature_classes.append(os.path.join(outputGeodatabase, cluster_feature_class))
    return (clustered_feature_classes, workerName, outputGeodatabase)

def create_cluster_layer(sourcePointFeatureClass, levelInfo):
    level, res = (levelInfo['level'], levelInfo['resolution'])

    integrate_name = 'integrate_{}'.format(level)
    cluster_name = 'clusters_{}'.format(level)
    output_name = 'clusters_{}_joined'.format(level)

    if not arcpy.Exists(output_name):
        arcpy.CopyFeatures_management(sourcePointFeatureClass, integrate_name)
        distance = "{} Meters".format(res * 10)
        try:
            arcpy.Integrate_management(integrate_name, distance)
        except:
            message = arcpy.GetMessage(0)
            if 'maximum tolerance exceeded' in message.lower():
                print 'maximum tolerance exceeded: {} - {}'.format(integrate_name, distance)
                return None
            else:
                raise

        cluster_features = arcpy.CollectEvents_stats(integrate_name, cluster_name)
        add_attributes(integrate_name, cluster_name, output_name, )

        arcpy.Delete_management(integrate_name)
        arcpy.Delete_management(cluster_name)

        arcpy.AddField_management(output_name, 'zoom_level', "SHORT")
        expression = "int({})".format(levelInfo['level'])
        arcpy.CalculateField_management(output_name, 'zoom_level', expression, "PYTHON")

    return output_name

def add_attributes(integratedFeatureClass, clustersFeatureClass, outputFeatureClass, mergeNumericFields=[], mergeTextFields=False):
    '''
    TODO: Add join option for string fields; 
    TOOD: Add customizable fields for summary

    summarizes attributes for clusters
    First - The first input value is used.
    Last -The last input value is used.
    Join -Merge the values together using the delimiter value to separate the values (only valid if the output field type is text).
    Min -The smallest input value is used (only valid if the output field type is numeric).
    Max -The largest input value is used (only valid if the output field type is numeric).
    Mean -The mean is calculated using the input values (only valid if the output field type is numeric).
    Median -The median is calculated using the input values (only valid if the output field type is numeric).
    Sum -The sum is calculated using the input values (only valid if the output field type is numeric).
    StDev -The standard deviation is calculated using the input values (only valid if the output field type is numeric).
    Count -The number of values included in statistical calculations. This counts each value except null values.
    '''
    fieldmappings = arcpy.FieldMappings()
    fieldmappings.addTable(clustersFeatureClass)

    integrated_fields = arcpy.ListFields(integratedFeatureClass)
    for f in integrated_fields:
        if f.type in ['Integer', 'SmallInteger', 'Double', 'Single']:
            template = f.name + '_{}'
            for m in ['Min', 'Max', 'Mean', 'Sum', 'Count', 'StDev', 'Median']:
                fieldmap = arcpy.FieldMap()
                fieldmap.addInputField(integratedFeatureClass, f.name)
                field = fieldmap.outputField
                field.name = template.format(m.lower())
                field.aliasName = template.format(m.lower())
                fieldmap.mergeRule = m
                fieldmap.outputField = field
                fieldmappings.addFieldMap(fieldmap)

    arcpy.SpatialJoin_analysis(clustersFeatureClass, integratedFeatureClass, outputFeatureClass, "#", "#", fieldmappings)

def create_membership_where(featureClass, fieldName, values):
    where_membership = ','.join(["'{}'".format(v) for v in values])
    return arcpy.AddFieldDelimiters(featureClass, fieldName) + " IN ({})".format(where_membership)

def get_unique_field_values(featureClass, fieldName, where=None, getValueFunction=None):
    unique_values = []
    with arcpy.da.SearchCursor(featureClass, [fieldName], where) as cursor:
        for r in cursor:
            if getValueFunction:
                unique_values.append(getValueFunction(r[0]))
            else:
                unique_values.append(str(r[0]))
    final_list = list(set(unique_values))
    return final_list

all_results = []
delete_files = []
def cluster_complete(results):
    clustered_feature_classes, workerName, delete_file = results
    print 'worker complete({})'.format(workerName)

    all_results += clustered_feature_classes
    delete_files.append(delete_file)

def run_clustering(featureClass, clusterField, outputDirectory, includeLevels=None):
    global all_results
    global delete_files
    arcpy.gp.overwriteOutput = True
    pool = multiprocessing.Pool(multiprocessing.cpu_count() - 1)
    for v in get_unique_field_values(featureClass, clusterField):
        where = create_membership_where(featureClass, clusterField, [v])
        worker_arguments = (featureClass, outputDirectory, None, where, v, includeLevels)
        pool.apply_async(cluster_worker, worker_arguments, callback=cluster_complete)
    pool.close()
    pool.join()

    final_output = featureClass + '_clustered'
    arcpy.CopyFeatures_management(all_results.pop(), final_output)
    arcpy.Append_management(all_results, final_output, "TEST")

    for d in delete_files:
        arcpy.Delete_management(d)

#=============================================================================================================
# TESTING
#=============================================================================================================
class TestClustering(unittest.TestCase):
    def test_run_clustering(self):
        feature_class = r'D:\data\EPA\FRS_INTERESTS.gdb\FACILITY_INTERESTS_WM'
        clusterField = 'STATE_CODE'
        outputDirectory = r'D:\data\EPA'
        levels = range(4, 16)
        run_clustering( feature_class, clusterField, outputDirectory, includeLevels=levels )

    def _add_attributes(self):
        arcpy.gp.overwriteOutput = True
        integrated_feature_class = r'C:\gis_working_directory\tempDE.gdb\integrate_10'
        clusters = r'C:\gis_working_directory\tempDE.gdb\clusters_10'
        outputDirectory = r'C:\gis_working_directory\tempDE.gdb\clusters_10_joined'
        add_attributes(integrated_feature_class, clusters, outputDirectory)

if __name__ == '__main__':
    unittest.main()

Great Sublime Text 2 Video

This a great resource for learning to Sublime Text 2:

Converting ArcGIS Feature Class to PyTables

Here is an example of converting a point geometry type ArcGIS Feature Class to a PyTables table. I decided to avoid using the arcpy.da.FeatureClassToNumpyArray, because I want more flexibility on how data is converted (e.g. geometry types, blobs, rasters, dates). Right now this only handles a point feature class, but I hope to add in other geometry types and make is into a module:

import arcpy
import tables
from tables import *

def feature_class_to_hdf5(featureClass, tableName, groupName='', indexFields=[]):
    '''
    TODO: add polygon and line geometry type support
    TODO: add blob field
    TODO: add raster support
    TODO: add Date logic
    '''
    try:
        #Validate feature class
        desc = arcpy.Describe(featureClass)
        if desc.dataType != 'FeatureClass':
            raise ValueError('featureClass must be a featureClass')
        if desc.featureType != 'Simple':
            raise ValueError('feature_class_to_table only support simple feature types')
        if desc.shapeType != 'Point':
            raise ValueError('sorry')

        #Construct output location for hdf5
        output_file_name = os.path.splitext(os.path.split(desc.catalogPath)[1])[0]
        if '.gdb' in desc.catalogPath:
            output_directory = os.path.split(os.path.split(desc.catalogPath)[0])[0]
        elif '.shp' in desc.catalogPath:
            output_directory = os.path.split(desc.catalogPath)[0]
        output_hd5_path = os.path.join(output_directory, output_file_name + '.h5')
        
        #Map ArcGIS Field Types to PyTable Field Types
        fields_dictionary = {}
        fields_dictionary['OID'] = Int32Col(dflt=1, pos = 0)
        fields_dictionary['X'] = Float64Col(dflt=1, pos = 1)
        fields_dictionary['Y'] = Float64Col(dflt=1, pos = 2)

        h5f_field_names = ['OID', 'X', 'Y']
        arc_field_names = ['OID@', 'SHAPE@X', 'SHAPE@Y']
        for f in desc.fields:
            if f.type in ['String','GUID']:
                fields_dictionary[str(f.name)] = StringCol(itemsize=f.length, dflt="n/a", pos=len(arc_field_names))
            elif f.type == 'Single':
                fields_dictionary[str(f.name)] = Float32Col(dflt=1, pos=len(arc_field_names))
            elif f.type == 'Double':
                fields_dictionary[str(f.name)] = Float64Col(dflt=1, pos=len(arc_field_names))
            elif f.type == 'SmallInteger':
                fields_dictionary[str(f.name)] = Int16Col(dflt=1, pos=len(arc_field_names))
            elif f.type == 'Integer':
                fields_dictionary[str(f.name)] = Int32Col(dflt=1, pos=len(arc_field_names))
            elif f.type == 'Date':
                fields_dictionary[str(f.name)] = Time64Col(dflt=1, pos=len(arc_field_names)) #Needs test...
            else:
                continue
            arc_field_names.append(str(f.name))
            h5f_field_names.append(str(f.name))

        #Create PyTables DescriptionClass
        DescriptionClass = type(tableName, (IsDescription,), fields_dictionary)

        with tables.openFile(output_hd5_path, 'w', title=featureClass) as h5_file:
            filters = Filters(complib='blosc', complevel=1)
            h5_feature_class = h5_file.createTable('/', tableName, DescriptionClass, filters=filters)
            with arcpy.da.SearchCursor(featureClass, field_names=arc_field_names) as cursor:
                for r in cursor:
                    new_row = h5_feature_class.row
                    for i, fieldName in enumerate(h5f_field_names):
                        try:
                            new_row[fieldName] = r[i]
                        except TypeError as e:
                            if 'unicode' in str(e).lower():
                                new_row[fieldName] = unicodedata.normalize('NFKD', r[i]).encode('ascii','ignore')
                            else:
                                raise
                    new_row.append()
            h5_file.flush()

    except Exception as e:
        print e
        raise

Funky Arcpy Cursors - functional programming examples

Python has great functional programming features with lambda expressions and first-class functions. There are many ways to cook a turkey, but here are some examples which rely heavily on a functional approach.

Map Functions to Feature Class Rows 

This example simply wraps a search cursor and executes functions. Not super helpful, but it shows the functional mindset and the use of functions as input parameters

def mapFunctionsToFeatureClassRows(featureClass, functionList):
    rows = arcpy.SearchCursor(inputFeatureClass, where)
    for r in rows:
        for f in functionList:
            f(r)
    del rows

Group Rows

def groupRows(inputFeatureClass, where, keyFunction, valueFunction):
    from collections import defaultdict

    groupings = defaultdict(list)
    rows = arcpy.SearchCursor(inputFeatureClass, where)
    for r in rows:
        groupings[keyFunction(r)].append(valueFunction(r))
    del rows
    return groupings

With 10.1+, you can utilize the arcpy.da library for significant performance increase:

def groupRows2(inputFeatureClass, fields, where, keyFunction, valueFunction):
    from collections import defaultdict
    groupings = defaultdict(list)
    with arcpy.da.SearchCursor(inputFeautreClass, fields, where) as cursor:
        for r in cursor:
            groupings[keyFunction(r)].append(valueFunction(r))
    return groupings

By convention, groupRows will call keyFunction(row) and valueFunction(row) supplying the current arcpy row object. The result will be a dictionary mapping keyFunction outputs to lists of valueFunction outputs:

if __name__ == '__main__':
    featureClass = arcpy.GetParameterAsText(0) or r'C:\python_working_directory\schools.gdb\school_points'
    schoolIdField = arcpy.GetParameterAsText(1) or r'schoolId'
    whereClause = None
    keyFunction = lambda r: r.getValue(stateSchoolId[:2])
    valueFunction = lambda r: r.getValue(stateSchoolId)

    groupsByDistrict = groupRows(features, None, keyFunction, valueFunction) 

>>> {'TX': ['TX201', 'TX202', 'TX203']}

 

Translate Append

Do you need to append to an existing feature class?  Basically the same data, but conflicting attribute tables?  Here's a functional solution:

def translateAppend(targetFeatureClass, appendFeatureClass, where, fieldTranslations, outputFeatureClass):
    import types
    
    arcpy.CopyFeatures_management(targetFeatureClass, outputFeatureClass)
    insertCursor = arcpy.InsertCursor(outputFeatureClass)
    appendRows = arcpy.SearchCursor(appendFeatureClass, where)
    
    for r in appendRows:
        newRow = insertCursor.newRow() 
        for k,v in fieldTranslations.items():
            value = None
            if isinstance(v, types.FunctionType):
                value = v(r)
            else:
                value = r.getValue(v)
            if value:
                newRow.setValue(k, value)
        insertCursor.insertRow(newRow)
    
    del newRow
    del appendRows
    del insertCursor

You supply the function with a dictionary whose keys are the fields your existing feature class. The values of the dict are functions which output the value for that field based on your new append data. The values can be lambda expressions, standard function (def statement). By convention, if you use a string value then that field from the append feature class is used:

def createGradeList(row):
    low_grade = row.getValue('low_grade')
    high_grade = row.getValue('high_grade')
    return range(low_grade, high_grade, 1)
   
#Map fields in input feature class to strings or functions
translations = {}

#Straight field-to-field name copy
translations['Shape'] = 'Shape'

#generate value using lambda expression
translations['state'] = lambda row: row.getValue(stateSchoolId)[:2] 

#lookup value in dictionary
translations['StateAbbrev'] = lambda row: state_abbrevs_index[row.getValue('FIPS')]

# generate value using function
translations['GradesCommaDelimited'] = createGradeList

translateAppend(feature_class, new_features, "1=1", translations, output_feature_class)

Get Unique Pixel Values By Channel (RGBA) using Python

The following code sample shows how to extract a list of unique pixel values for each channel in an image using the Python Imaging Library (PIL).

def getUniquePixelValuesByChannel(imagePath, hexBase=True):
    '''
    returns lists of unique pixel value for each channel in image(r,g,b,a)
    '''
    image = Image.open(imagePath)
    histogram = image.histogram()
    channels = [histogram[p : p + 256] for p in range(0, len(histogram), 256)]
    unique_pixels = []
    for c in channels:
        unique_pixels.append([hexBase and hex(p) or p for p in range(len(c)) if c[p] > 0])
    return unique_pixels

This function uses the image.histogram() method to get a list of the frequencies of each pixel value across each color channel (red, green, blue, alpha).  For a 4-band image, the resulting list is 1024 values long (256 for each channel) with the index of each item representing the color value in that channel.  

Below is an example of the function in-action returning the unique values for each channel in either hex or 256:

if __name__ == '__main__':

    imagePath = r'lookAtMyPixelValuesInHex.png'

    pixelValuesInHex = getUniquePixelValuesByChannel(imagePath, True)
    print 'hex result:', pixelValuesInHex
    '''
    hex result: [
        ['0x0', '0x7', '0xf', '0x17', '0x1f', '0x27', '0x37', '0xfd'], 
        ['0x7', '0xf', '0x17', '0x1f', '0x27', '0x2f', '0xfd'], 
        ['0x7', '0xf', '0x17', '0x1f', '0xfd'], 
        ['0x7', '0xf', '0x17', '0x1f', '0x27', '0x2f', '0x37', '0xfd']
    ]

    '''

    pixelValuesIn256 = getUniquePixelValuesByChannel(imagePath, False)
    print '256 color result:', pixelValuesIn256
    '''
    256 color result: [
        [0, 7, 15, 23, 31, 39, 55, 253], 
        [7, 15, 23, 31, 39, 47, 253], 
        [7, 15, 23, 31, 253], 
        [7, 15, 23, 31, 39, 47, 55, 253]
    ]
    '''

Normalize to USPS Street Abbreviations using Python

Here's a function i recently wrote for converting street suffixes using common street suffixes mapped to the usps standard.  You can find the list on the usps site.


Came in handy for some fuzzy matching of address strings...

def normalizeStreetSuffixes(inputValue):
        '''
        Use common abbreviations -> USPS standardized abbreviation to replace common street suffixes

        Obtains list from https://www.usps.com/send/official-abbreviations.htm
        '''
        usps_street_abbreviations = {'trpk': 'tpke', 'forges': 'frgs', 'bypas': 'byp', 'mnr': 'mnr', 'viaduct': 'via', 'mnt': 'mt',
         'lndng': 'lndg', 'vill': 'vlg', 'aly': 'aly', 'mill': 'ml', 'pts': 'pts', 'centers': 'ctrs', 'row': 'row', 'cnter': 'ctr',
          'hrbor': 'hbr', 'tr': 'trl', 'lndg': 'lndg', 'passage': 'psge', 'walks': 'walk', 'frks': 'frks', 'crest': 'crst', 'meadows': 'mdws',
           'freewy': 'fwy', 'garden': 'gdn', 'bluffs': 'blfs', 'vlg': 'vlg', 'vly': 'vly', 'fall': 'fall', 'trk': 'trak', 'squares': 'sqs',
            'trl': 'trl', 'harbor': 'hbr', 'frry': 'fry', 'div': 'dv', 'straven': 'stra', 'cmp': 'cp', 'grdns': 'gdns', 'villg': 'vlg',
             'meadow': 'mdw', 'trails': 'trl', 'streets': 'sts', 'prairie': 'pr', 'hts': 'hts', 'crescent': 'cres', 'pass': 'pass',
              'ter': 'ter', 'port': 'prt', 'bluf': 'blf', 'avnue': 'ave', 'lights': 'lgts', 'rpds': 'rpds', 'harbors': 'hbrs',
               'mews': 'mews', 'lodg': 'ldg', 'plz': 'plz', 'tracks': 'trak', 'path': 'path', 'pkway': 'pkwy', 'gln': 'gln',
                'bot': 'btm', 'drv': 'dr', 'rdg': 'rdg', 'fwy': 'fwy', 'hbr': 'hbr', 'via': 'via', 'divide': 'dv', 'inlt': 'inlt',
                 'fords': 'frds', 'avenu': 'ave', 'vis': 'vis', 'brk': 'brk', 'rivr': 'riv', 'oval': 'oval', 'gateway': 'gtwy',
                  'stream': 'strm', 'bayoo': 'byu', 'msn': 'msn', 'knoll': 'knl', 'expressway': 'expy', 'sprng': 'spg',
                   'flat': 'flt', 'holw': 'holw', 'grden': 'gdn', 'trail': 'trl', 'jctns': 'jcts', 'rdgs': 'rdgs',
                    'tunnel': 'tunl', 'ml': 'ml', 'fls': 'fls', 'flt': 'flt', 'lks': 'lks', 'mt': 'mt', 'groves': 'grvs',
                     'vally': 'vly', 'ferry': 'fry', 'parkway': 'pkwy', 'radiel': 'radl', 'strvnue': 'stra', 'fld': 'fld',
                      'overpass': 'opas', 'plaza': 'plz', 'estate': 'est', 'mntn': 'mtn', 'lock': 'lck', 'orchrd': 'orch',
                       'strvn': 'stra', 'locks': 'lcks', 'bend': 'bnd', 'kys': 'kys', 'junctions': 'jcts', 'mountin': 'mtn',
                        'burgs': 'bgs', 'pine': 'pne', 'ldge': 'ldg', 'causway': 'cswy', 'spg': 'spg', 'beach': 'bch', 'ft': 'ft',
                         'crse': 'crse', 'motorway': 'mtwy', 'bluff': 'blf', 'court': 'ct', 'grov': 'grv', 'sprngs': 'spgs',
                          'ovl': 'oval', 'villag': 'vlg', 'vdct': 'via', 'neck': 'nck', 'orchard': 'orch', 'light': 'lgt',
                           'sq': 'sq', 'pkwy': 'pkwy', 'shore': 'shr', 'green': 'grn', 'strm': 'strm', 'islnd': 'is',
                            'turnpike': 'tpke', 'stra': 'stra', 'mission': 'msn', 'spngs': 'spgs', 'course': 'crse',
                             'trafficway': 'trfy', 'terrace': 'ter', 'hway': 'hwy', 'avenue': 'ave', 'glen': 'gln',
                              'boul': 'blvd', 'inlet': 'inlt', 'la': 'ln', 'ln': 'ln', 'frst': 'frst', 'clf': 'clf',
                               'cres': 'cres', 'brook': 'brk', 'lk': 'lk', 'byp': 'byp', 'shoar': 'shr', 'bypass': 'byp',
                                'mtin': 'mtn', 'ally': 'aly', 'forest': 'frst', 'junction': 'jct', 'views': 'vws', 'wells': 'wls', 'cen': 'ctr',
                                 'exts': 'exts', 'crt': 'ct', 'corners': 'cors', 'trak': 'trak', 'frway': 'fwy', 'prarie': 'pr', 'crossing': 'xing',
                                  'extn': 'ext', 'cliffs': 'clfs', 'manors': 'mnrs', 'ports': 'prts', 'gatewy': 'gtwy', 'square': 'sq', 'hls': 'hls',
                                   'harb': 'hbr', 'loops': 'loop', 'mdw': 'mdw', 'smt': 'smt', 'rd': 'rd', 'hill': 'hl', 'blf': 'blf',
                                    'highway': 'hwy', 'walk': 'walk', 'clfs': 'clfs', 'brooks': 'brks', 'brnch': 'br', 'aven': 'ave',
                                     'shores': 'shrs', 'iss': 'iss', 'route': 'rte', 'wls': 'wls', 'place': 'pl', 'sumit': 'smt', 'pines': 'pnes',
                                      'trks': 'trak', 'shoal': 'shl', 'strt': 'st', 'frwy': 'fwy', 'heights': 'hts', 'ranches': 'rnch',
                                       'boulevard': 'blvd', 'extnsn': 'ext', 'mdws': 'mdws', 'hollows': 'holw', 'vsta': 'vis', 'plains': 'plns',
                                        'station': 'sta', 'circl': 'cir', 'mntns': 'mtns', 'prts': 'prts', 'shls': 'shls', 'villages': 'vlgs',
                                         'park': 'park', 'nck': 'nck', 'rst': 'rst', 'haven': 'hvn', 'turnpk': 'tpke', 'expy': 'expy', 'sta': 'sta',
                                          'expr': 'expy', 'stn': 'sta', 'expw': 'expy', 'street': 'st', 'str': 'st', 'spurs': 'spur', 'crecent': 'cres',
                                           'rad': 'radl', 'ranch': 'rnch', 'well': 'wl', 'shoals': 'shls', 'alley': 'aly', 'plza': 'plz', 'medows': 'mdws',
                                            'allee': 'aly', 'knls': 'knls', 'ests': 'ests', 'st': 'st', 'anx': 'anx', 'havn': 'hvn', 'paths': 'path', 'bypa': 'byp',
                                             'spgs': 'spgs', 'mills': 'mls', 'parks': 'park', 'byps': 'byp', 'flts': 'flts', 'tunnels': 'tunl', 'club': 'clb', 'sqrs': 'sqs',
                                              'hllw': 'holw', 'manor': 'mnr', 'centre': 'ctr', 'track': 'trak', 'hgts': 'hts', 'rnch': 'rnch', 'crcle': 'cir', 'falls': 'fls',
                                               'landing': 'lndg', 'plaines': 'plns', 'viadct': 'via', 'gdns': 'gdns', 'gtwy': 'gtwy', 'grove': 'grv', 'camp': 'cp', 'tpk': 'tpke',
                                                'drive': 'dr', 'freeway': 'fwy', 'ext': 'ext', 'points': 'pts', 'exp': 'expy', 'ky': 'ky', 'courts': 'cts', 'pky': 'pkwy', 'corner': 'cor',
                                                 'crssing': 'xing', 'mnrs': 'mnrs', 'unions': 'uns', 'cyn': 'cyn', 'lodge': 'ldg', 'trfy': 'trfy', 'circle': 'cir', 'bridge': 'brg',
                                                  'dl': 'dl', 'dm': 'dm', 'express': 'expy', 'tunls': 'tunl', 'dv': 'dv', 'dr': 'dr', 'shr': 'shr', 'knolls': 'knls', 'greens': 'grns',
                                                   'tunel': 'tunl', 'fields': 'flds', 'common': 'cmn', 'orch': 'orch', 'crk': 'crk', 'river': 'riv', 'shl': 'shl', 'view': 'vw',
                                                    'crsent': 'cres', 'rnchs': 'rnch', 'crscnt': 'cres', 'arc': 'arc', 'btm': 'btm', 'blvd': 'blvd', 'ways': 'ways', 'radl': 'radl',
                                                     'rdge': 'rdg', 'causeway': 'cswy', 'parkwy': 'pkwy', 'juncton': 'jct', 'statn': 'sta', 'gardn': 'gdn', 'mntain': 'mtn',
                                                      'crssng': 'xing', 'rapid': 'rpd', 'key': 'ky', 'plns': 'plns', 'wy': 'way', 'cor': 'cor', 'ramp': 'ramp', 'throughway': 'trwy',
                                                       'estates': 'ests', 'ck': 'crk', 'loaf': 'lf', 'hvn': 'hvn', 'wall': 'wall', 'hollow': 'holw', 'canyon': 'cyn', 'clb': 'clb',
                                                        'cswy': 'cswy', 'village': 'vlg', 'cr': 'crk', 'trce': 'trce', 'cp': 'cp', 'cv': 'cv', 'ct': 'cts', 'pr': 'pr', 'frg': 'frg',
                                                         'jction': 'jct', 'pt': 'pt', 'mssn': 'msn', 'frk': 'frk', 'brdge': 'brg', 'cent': 'ctr', 'spur': 'spur', 'frt': 'ft', 'pk': 'park',
                                                          'fry': 'fry', 'pl': 'pl', 'lanes': 'ln', 'gtway': 'gtwy', 'prk': 'park', 'vws': 'vws', 'stravenue': 'stra', 'lgt': 'lgt',
                                                           'hiway': 'hwy', 'ctr': 'ctr', 'prt': 'prt', 'ville': 'vl', 'plain': 'pln', 'mount': 'mt', 'mls': 'mls', 'loop': 'loop',
                                                            'riv': 'riv', 'centr': 'ctr', 'is': 'is', 'prr': 'pr', 'vl': 'vl', 'avn': 'ave', 'vw': 'vw', 'ave': 'ave', 'spng': 'spg',
                                                             'hiwy': 'hwy', 'dam': 'dm', 'isle': 'isle', 'crcl': 'cir', 'sqre': 'sq', 'jct': 'jct', 'jctn': 'jct', 'mountain': 'mtn',
                                                              'keys': 'kys', 'parkways': 'pkwy', 'drives': 'drs', 'tunl': 'tunl', 'jcts': 'jcts', 'knl': 'knl', 'center': 'ctr',
                                                               'driv': 'dr', 'tpke': 'tpke', 'sumitt': 'smt', 'canyn': 'cyn', 'ldg': 'ldg', 'harbr': 'hbr', 'rest': 'rst', 'shoars': 'shrs',
                                                                'vist': 'vis', 'gdn': 'gdn', 'islnds': 'iss', 'hills': 'hls', 'cresent': 'cres', 'point': 'pt', 'lake': 'lk', 'vlly': 'vly',
                                                                 'strav': 'stra', 'crossroad': 'xrd', 'bnd': 'bnd', 'strave': 'stra', 'stravn': 'stra', 'knol': 'knl', 'vlgs': 'vlgs',
                                                                  'forge': 'frg', 'cntr': 'ctr', 'cape': 'cpe', 'height': 'hts', 'lck': 'lck', 'highwy': 'hwy', 'trnpk': 'tpke', 'rpd': 'rpd',
                                                                   'boulv': 'blvd', 'circles': 'cirs', 'valleys': 'vlys', 'vst': 'vis', 'creek': 'crk', 'mall': 'mall', 'spring': 'spg',
                                                                    'brg': 'brg', 'holws': 'holw', 'lf': 'lf', 'est': 'est', 'xing': 'xing', 'trace': 'trce', 'bottom': 'btm',
                                                                     'streme': 'strm', 'isles': 'isle', 'circ': 'cir', 'forks': 'frks', 'burg': 'bg', 'run': 'run', 'trls': 'trl',
                                                                      'radial': 'radl', 'lakes': 'lks', 'rue': 'rue', 'vlys': 'vlys', 'br': 'br', 'cors': 'cors', 'pln': 'pln',
                                                                       'pike': 'pike', 'extension': 'ext', 'island': 'is', 'frd': 'frd', 'lcks': 'lcks', 'terr': 'ter',
                                                                        'union': 'un', 'extensions': 'exts', 'pkwys': 'pkwy', 'islands': 'iss', 'road': 'rd', 'shrs': 'shrs',
                                                                         'roads': 'rds', 'glens': 'glns', 'springs': 'spgs', 'missn': 'msn', 'ridge': 'rdg', 'arcade': 'arc',
                                                                          'bayou': 'byu', 'crsnt': 'cres', 'junctn': 'jct', 'way': 'way', 'valley': 'vly', 'fork': 'frk',
                                                                           'mountains': 'mtns', 'bottm': 'btm', 'forg': 'frg', 'ht': 'hts', 'ford': 'frd', 'hl': 'hl',
                                                                            'grdn': 'gdn', 'fort': 'ft', 'traces': 'trce', 'cnyn': 'cyn', 'cir': 'cir', 'un': 'un', 'mtn': 'mtn',
                                                                             'flats': 'flts', 'anex': 'anx', 'gatway': 'gtwy', 'rapids': 'rpds', 'villiage': 'vlg', 'flds': 'flds',
                                                                              'coves': 'cvs', 'rvr': 'riv', 'av': 'ave', 'pikes': 'pike', 'grv': 'grv', 'vista': 'vis', 'pnes': 'pnes',
                                                                               'forests': 'frst', 'field': 'fld', 'branch': 'br', 'grn': 'grn', 'dale': 'dl', 'rds': 'rds', 'annex': 'anx',
                                                                                'sqr': 'sq', 'cove': 'cv', 'squ': 'sq', 'skyway': 'skwy', 'ridges': 'rdgs', 'hwy': 'hwy', 'tunnl': 'tunl',
                                                                                 'underpass': 'upas', 'cliff': 'clf', 'lane': 'ln', 'land': 'land', 'bch': 'bch', 'dvd': 'dv', 'curve': 'curv',
                                                                                  'cpe': 'cpe', 'summit': 'smt', 'gardens': 'gdns'}
        words = inputValue.split()
        for w in words:
            if w in usps_street_abbreviations.keys():
                inputValue = inputValue.replace(w, usps_street_abbreviations[w])
        return inputValue

Online/Offline Mapping: ESRI Dev Summit 2012

I had a great time out in Palm Springs this year.  Saw some great presentations and will be definitely going back over the videos.  Looks like ESRI has released many of them in their re-vamped video site.

I presented this year on converting an existing Flex application to AIR and running it offline:

Here's a link to my presentation. 

Online/Offline Mapping: Finding tile urls based on a geographic extent

This is my second post in a series on developing web mapping applications which also function offline.  My last post covered map tile url schemas for various basemap providers.  While ESRI, Open Street Map, Google, Bing, and Yahoo / Nokia structure their urls differently, they all use the same column (x), row (y), and zoom level (z) values for a indicating tile position within a quad tree.  If you haven't skimmed over that post, or this sounds like Klingon, I recommend checking that post out. 

The here's the next challenge: take a geographic area of interest, and find urls for all the tiles which covering that area (at zoom levels of interest). 

My area of interest will be the fine City of New Orleans

    '''
    The City of New Orleans
    xmin: -90.283741
    ymin: 29.890626
    xmax: -89.912952
    ymax: 30.057766
    '''

Disclaimer:  Persisting map tiles locally breaks the terms of service for many popular basemap providers like Google, Bing, and ArcGIS Online.  I'm not a lawyer, nor do I aspire to be one, but if you start storing tiles offline without permission, then you do so at your own risk. 

Helper Methods for Finding Tiles

There are a bunch of functions you will need at your disposal to accomplish this task. These functions can be found in Joe Schwartz's Awesome White Paper on the Bing Maps Tiling System.  This is one of the best explanations of map tile systems I've found on the web.  Joe goes over the mercator projection, along with explanations of ground resolution, map scale, and pixel coordinates. Really really great stuff and I urge everybody to really dig into it.

At the end of the article, Joe provides sample code which contains the major functions you will need to dominate map tile systems.  While I've grown to really like C# development, I almost always translate interesting code back into warm and friendly Python.   Here's Joe's TileSystem class translated into Python with a couple extra methods for additional abstraction:

'''
Mostly direct port of awesome article by Joe Schwartz - http://msdn.microsoft.com/en-us/library/bb259689.aspx
'''
from __future__ import division
import math

class TileUtils(object):

    earth_radius = 6378137
    min_lat = -85.05112878
    max_lat = 85.05112878
    min_lng = -180
    max_lng = 180

    def clipValue(self, value, minValue, maxValue):
        '''
        Makes sure that value is within a specific range.
        If not, then the lower or upper bounds is returned
        '''
        return min(max(value, minValue), maxValue)

    def getMapDimensionsByZoomLevel(self, zoomLevel):
        '''
        Returns the width/height in pixels of the entire map
        based on the zoom level.
        '''
        return 256 << zoomLevel

    def getGroundResolution(self, latitude, level):
        '''
        returns the ground resolution for based on latitude and zoom level.
        '''
        latitude = self.clipValue(latitude, self.min_lat, self.max_lat);
        mapSize = self.getMapDimensionsByZoomLevel(level)
        return math.cos(latitude * math.pi / 180) * 2 * math.pi * self.earth_radius / mapSize

    def getMapScale(self, latitude, level, dpi=96):
        '''
        returns the map scale on the dpi of the screen
        '''
        dpm = dpi / 0.0254 #convert to dots per meter
        return self.getGroundResolution(latitude, level) * dpm

    def convertLatLngToPixelXY(self, lat, lng, level):
        '''
        returns the x and y values of the pixel corresponding to a latitude and longitude.
        '''
        mapSize = self.getMapDimensionsByZoomLevel(level)
       
        lat = self.clipValue(lat, self.min_lat, self.max_lat)
        lng = self.clipValue(lng, self.min_lng, self.max_lng)

        x = (lng + 180) / 360
        sinlat = math.sin(lat * math.pi / 180)
        y = 0.5 - math.log((1 + sinlat) / (1 - sinlat)) / (4 * math.pi)

        pixelX = int(self.clipValue(x * mapSize + 0.5, 0, mapSize - 1))
        pixelY = int(self.clipValue(y * mapSize + 0.5, 0, mapSize - 1))
        return (pixelX, pixelY)

    def convertPixelXYToLngLat(self, pixelX, pixelY, level):
        '''
        converts a pixel x, y to a latitude and longitude.
        '''
        mapSize = self.getMapDimensionsByZoomLevel(level)
        x = (self.clipValue(pixelX, 0, mapSize - 1) / mapSize) - 0.5
        y = 0.5 - (self.clipValue(pixelY, 0, mapSize - 1) / mapSize)

        lat = 90 - 360 * math.atan(math.exp(-y * 2 * math.pi)) / math.pi
        lng = 360 * x

        return (lng, lat)

    def convertPixelXYToTileXY(self, pixelX, pixelY):
        '''
        Converts pixel XY coordinates into tile XY coordinates of the tile containing
        '''
        return(int(pixelX / 256), int(pixelY / 256))

    def convertTileXYToPixelXY(self, tileX, tileY):
        '''
        Converts tile XY coordinates into pixel XY coordinates of the upper-left pixel
        '''
        return(tileX * 256, tileY * 256)

    def tileXYZToQuadKey(self, x, y, z):
        '''
        Computes quadKey value based on tile x, y and z values.
        '''
        quadKey = ''
        for i in range(z, 0, -1):
            digit = 0
            mask = 1 << (i - 1)
            if(x & mask) != 0:
                digit += 1
            if(y & mask) != 0:
                digit += 2
            quadKey += str(digit)
        return quadKey

    def quadKeyToTileXYZ(self, quadKey):
        '''
        Computes tile x, y and z values based on quadKey.
        '''
        tileX = 0
        tileY = 0
        tileZ = len(quadKey)

        for i in range(tileZ, 0, -1):
            mask = 1 << (i - 1)
            value = quadKey[tileZ - i]

            if value == '0':
                continue

            elif value == '1':
                tileX |= mask

            elif value == '2':
                tileY |= mask

            elif value == '3':
                tileX |= mask
                tileY |= mask

            else:
                raise Exception('Invalid QuadKey')

        return (tileX, tileY, tileZ)

    def convertLngLatToTileXY(self, lng, lat, level):
        pixelX, pixelY = self.convertLatLngToPixelXY(lat, lng, level)
        return self.convertPixelXYToTileXY(pixelX, pixelY)

    def getTileOrigin(self, tileX, tileY, level):
        '''
        Returns the upper-left hand corner lat/lng for a tile
        '''
        pixelX, pixelY = self.convertTileXYToPixelXY(tileX, tileY)
        lng, lat = self.convertPixelXYToLngLat(pixelX, pixelY, level)
        return (lat, lng)

The function convertLngLatToTileXY really sums up the class for our purposes. This function will help get the tile x/y values for the sides of the bounding box.  Using these functions, we can write an additional helper class to generate all the tile urls within our bounding box based off of a template tile url like the ones mentioned in my last post:

class TileFinder(object):

    def __init__(self, tileUtils, tileTemplate):
        '''
        Tile template is a template url which use {{x}}, {{y}}, {{z}} which will be replaced by
        their respective values.
        '''
        self.tileUtils = tileUtils
        self.tileTemplate = tileTemplate
    
    def getTileUrlsByLatLngExtent(self, xmin, ymin, xmax, ymax, level):
        '''
        Returns a list of tile urls by extent
        '''
        #Upper-Left Tile
        tileXMin, tileYMin = self.tileUtils.convertLngLatToTileXY(xmin, ymax, level)

        #Lower-Right Tile
        tileXMax, tileYMax = self.tileUtils.convertLngLatToTileXY(xmax, ymin, level)

        tileUrls = []

        for y in range(tileYMax, tileYMin - 1, -1):
            for x in range(tileXMin, tileXMax + 1, 1):
                tileUrls.append(self.createTileUrl(x, y, level)) 

        return tileUrls

    def createTileUrl(self, x, y, z):
        '''
        returns new tile url based on template
        '''
        return self.tileTemplate.replace('{{x}}', str(x)).replace('{{y}}', str(y)).replace('{{z}}', str(z))

Check out the constructor here.  It has two required arguments, one of which is a tile template.  This tile template uses {{x}}, {{y}}, and {{z}} values to indicate where substitutions should take place.  These substitutions happen in createTileUrl which is called for each x,y pair generated in getTileUrlsByLatLngExtent.

So now we have our helper classes for generating a list of tile urls based on an extent.  The last value to complete the trifecta is the zoom level.  I chose level 11 since my area of interest is fairly small (city level).  Here is a snippet for running the code:

    '''
    The City of New Orleans
    xmin: -90.283741
    ymin: 29.890626
    xmax: -89.912952
    ymax: 30.057766
    '''

    import urllib2
    
    #constructing the template url for arcgis online
    template = 'http://services.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer/tile/{{z}}/{{y}}/{{x}}.png'

    tileUtils = TileUtils()
    tileFinder = TileFinder(tileUtils, template)
    tile_urls = tileFinder.getTileUrlsByLatLngExtent(-90.283741, 29.890626, -89.912952, 30.057766, 11)

    for t in tile_urls: urllib2.urlopen(t).read()

So that's about all there is to it.  Below you can see the six tile urls that were produced in the code snippet above.  A special thanks to Stephen Ansari (@sansari22) who really got me started on learning tile schemas and provided me sample code to get going.

Online/Offline Mapping: Map Tiles and their URL Schemas

I've recently been developing online/offline mapping applications for use in places which simply don't have access to internet connections.  I'm currently preparing a presentation discussing topics related to bring online application to offline users, and in doing so will be blogging on related topics.  The first part of the series will be dedicated to working with tile map services.

One of the main challenges of bringing an online map to an offline environment is maintaining access to map tiles. From imagery to streets, tiled map services provide the necessary geographic context and semantic detail to visualize additional spatial layers.  While hard-drive space is still a major limiting factor in ability to store tiles, a significant number of tiles can be stored locally to allow use of an application for limited geographic area.

Disclaimer:  Persisting map tiles locally breaks the terms of service for many popular basemap providers like Google, Bing, and ArcGIS Online.  I'm not a lawyer, nor do I aspire to be one, but if you start storing tiles offline without permission, then you do so at your own risk.  With that said, let's look at how to get started in storing map tiles for offline use.


Map Tiles and their URL Schemas

If you want to save down map tiles, you need to first know the address of each tile.  Map tile addresses take the form of urls (uniform address locators) which web maps access via HTTP GET requests.  Map tile service providers structure their map tile urls differently, but they all use column (x), row (y), and zoom level (z) values somewhere in the url.  Most providers have settled on using the Web Mercator projection (sad but true), and most have used the same origin for tiling and tile dimensions which means that the zoom level, row and column ids are the same across providers.  Dimensions are almost all 256 by 256 pixels and are in png or jpeg format.

Below I have listed requests for a single tile covering the same geographic area (Pacific Ocean / California) for several basemap tile providers including: Open Street Map, ArcGIS Online, Google Maps, Bing Maps, Yahoo! Maps and Nokia.  Here are the details for the tile which is common to ALL providers and demonstrates how easy it can be to locate a tile for any of the providers below:

/*
Tile Details:
COLUMN (X): 2
ROW (Y): 6
ZOOM LEVEL (Z): 4
DIMENSIONS: 256 x 256
*/

Open Street Map:

O.S.M. uses a simple directory structure for storing its map tiles which we will see repeated in many of the other service providers.  While you may think that "Open" means you can freely scrape as many tiles as you want...remember that OSM is hosted by volunteers with limited server resources.  If you want to scrape a large number of tiles, check out hosting your own OSM tile server.  I used Michal Migurski's TileDrawer to quickly setup my own OSM EC2 instance.  Tile Drawer is great and doesn't even require you to SSH into the machine!  Just start it up and start scraping.

/*
Open Street Map Template:
http://{SERVER}/{ZOOM_LEVEL}/{COLUMN}/{ROW}.png

Open Street Map Example:
http://c.tile.openstreetmap.org/4/2/6.png
*/

ArcGIS Server Map Service:

ArcGIS Server is the most powerful, out-of-the-box mapping server currently available.  AGS uses a similiar tile directory structure as OSM, but reverses the column and row values in the url.  ESRI is very gracious with providing access to tons of free tile services.

/*
ArcGIS Server Tiled Map Service Template:
http://{SERVER}/tile/{ZOOM_LEVEL}/{ROW}/{COLUMN}.png

ArcGIS Server Tiled Map Service Example:
http://services.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer/tile/4/6/2.png
*/

Google Maps

Google uses a slightly different structure which involves url parameters.  They also include what appears to be a locale parameter which I'm guessing is used to change the language of label (although I have not yet confirmed).

/*
Google Maps Template:
http://{SERVER}/vt/lyrs={layer_id}&hl={locale}&x={COLUMN}&y={ROW}&z={ZOOM_LEVEL}

Google Maps Example: 
http://mt0.google.com/vt/lyrs=m@169000000&hl=en&x=2&y=6&z=4&s=Ga
*/

Bing Map (Formerly Microsoft Virtual Earth)

Bing has an interesting way to access tile which is using a quadkey.  A "quadkey" is basically the position of the tile within the "quadtree".  The real reason why Web Mercator become so successful is that it turns the earth into a square which can be recursively tiled by spliting tiles into 4 to obtain tiles for the next zoom level.  I will be covering how to generate quadkeys in coming posts, but for now, here is a python snippet which I used to get the tile below:

def tileXYZToQuadKey(x, y, z):
    quadKey = ''
    for i in range(z, 0, -1):
        digit = 0
        mask = 1 << (i - 1)
        if(x & mask) != 0:
            digit += 1
        if(y & mask) != 0:
            digit += 2
        quadKey += str(digit)
    return quadKey

>>> print tileXYZToQuadKey(2, 6, 4)
0230
>>> 

I wrote the function above based on some C# that was provided in the bing maps tiling specification which you can read more about here: Bing Maps tile specification. 

/*
Bing Maps Tile Template:
http://{SERVER}/tiles/a{QUAD_KEY}.jpeg?token={TOKEN}&mkt={LOCALE}&g={??}

Bing Maps Tile Example:
http://t0.tiles.virtualearth.net/tiles/a0230.jpeg?g=854&mkt=en-US&token=Anz84uRE1RULeLwuJ0qKu5amcu5rugRXy1vKc27wUaKVyIv1SVZrUjqaOfXJJoI0

Awesome Explanation: http://msdn.microsoft.com/en-us/library/bb259689.aspx
*/

Ovi Maps (Yahoo! / Nokia)
Yahoo! and Nokia use an API called Ovi Maps which I don't know much about (but I'm guessing Ovi was acquired by Nokia).  The tile structure is straight forward.

/*
Ovi Maps (Nokia and Yahoo! Maps) Template:
http://{SERVER}/maptiler/v2/maptile/???/{STYLE}/{ZOOM_LEVEL}/{COLUMN}/{ROW}/{SIZE}/{IMG_FORMAT}?token={TOKEN}&lg={LOCALE}

Ovi Maps (Nokia and Yahoo! Maps) Example:
http://4.maptile.lbs.ovi.com/maptiler/v2/maptile/279af375be/normal.day/4/2/6/256/png8?lg=ENG&token=TrLJuXVK62IQk0vuXFzaig%3D%3D&requestid=yahoo.prod&app_id=eAdkWGYRoc4RfxVo0Z4B
*/

Yandex Maps


Yandex has really nice maps and great to see soCal in cyrillic.

/*
Yandex Maps Template:
http://vec04.maps.yandex.net/tiles?l=map&v=2.26.0&x={{X}}&y={{y}}&z={{z}}&lang=ru-RU

Yandex Maps Example:
http://vec04.maps.yandex.net/tiles?l=map&v=2.26.0&x=2&y=6&z=4&lang=ru-RU
*/

Well that about concludes my post on map tile url schemas.  Look for my upcoming blog post where I will cover how to generate a list of tile x,y, z values for a specified zoom level based on a given extent.

Simple Source Code Count (Python)

Here is a simple way to count your lines of source code using python.

The command-line tool takes a source directory and a comma-delimited list of valid file extensions.

#Source Lines of Code ----------------------------------------------

import os
import sys

from optparse import OptionParser

def main(source_directory, extensions):
    files = get_file_path_list(source_directory, extensions)
    counts = [get_line_count(f) for f in files]
    
    print 'file count: %i' % len(files)
    print 'line count: %i' % sum(counts)

def get_file_path_list(root_directory, extensions):
    files = []
    for directory, dirnames, filenames in os.walk(root_directory):
        for f in filenames:
            ext = os.path.splitext(f)[1]
            if ext in extensions:
                files.append(os.path.join(directory, f))
    return files

def get_line_count(path):
    f = open(path, 'r')
    lineCount = len(f.readlines())
    f.close()
    return lineCount

def parseListOption(option, opt, value, parser):
    setattr(parser.values, option.dest, value.split(','))


if __name__ == '__main__':
    
    parser = OptionParser(usage="usage: %prog [options] filename",
                          version="%prog 1.0")
    
    parser.add_option("-d", "--dir",
                      dest="source_directory",
                      help="Path to directory of source code")

    parser.add_option("-e", "--exts",
                      type='string',
                      action='callback',
                      dest="extensions",
                      callback=parseListOption,
                      help="Path to directory of source code")
    
    (options, args) = parser.parse_args()

    main(options.source_directory, options.extensions)

I'm on a windows environment, so I run the tool with the following command on powershell:

PS > python countDirectoryLineCount.py -d source_directory -e extensions (.py,.as,.mxml,.cs,.java)

import this

>>> import this
The Zen of Python, by Tim Peters

Beautiful is better than ugly.
Explicit is better than implicit.
Simple is better than complex.
Complex is better than complicated.
Flat is better than nested.
Sparse is better than dense.
Readability counts.
Special cases aren't special enough to break the rules.
Although practicality beats purity.
Errors should never pass silently.
Unless explicitly silenced.
In the face of ambiguity, refuse the temptation to guess.
There should be one-- and preferably only one --obvious way to do it.
Although that way may not be obvious at first unless you're Dutch.
Now is better than never.
Although never is often better than *right* now.
If the implementation is hard to explain, it's a bad idea.
If the implementation is easy to explain, it may be a good idea.
Namespaces are one honking great idea -- let's do more of those!

ESRI Developer Summit Python Presentation

Back in March, I had the pleasure of attending the ESRI Developer Summit in Palm Springs where I gave the following presentation on Becoming a Python Developer...

I came to Python as a GIS Analyst without much experience in coding.  I stayed mostly within the arcgisscripting module and only later ventured into the Python Standard Library and third-party libraries.

In the video linked above, I talk a focused aspects of Python which GIS Analysts using the language should be aware of and how to get a development environment up and running. 

Python Geoprocessing Vs. C# Server Object Extensions (AGS 10)

The following deals with development of geoprocessing services specifically for ArcGIS Server 10

Introduction Video I made for Penn State's MGIS Program

Over the past weeks, I have written several server-side tools to support web mapping applications as part of my capstone project for Penn State's MGIS Program.  The following is a general overview of what I've found to be the strengths and weaknesses of  Geoprocessing Services in Python and Server Object Extensions in C#.

Over the coming weeks, I will be releasing the code and discussions of the individual tools in separate posts. 

Geoprocessing Services (Python)

For developers of web-based geoprocessing tools for ArcGIS Server 10,  the arcpy Python library is now a popular alternative to the more complex ArcObjects library.  With arcpy, you can access a decent percent of the ArcGIS tools with a clear syntax for linking them together. Knowledge of object-oriented programming is not required, and you can integrate GIS with an entire world of Python third-party libraries.  To web-enable arcpy, scripts are added to toolboxes which get published to the server.

Python Pros

• Accelerated development time
• ArcGIS includes many Python code samples
• Easy to debug business logic

In spite of their ease of development and clear syntax, aspects of python geoprocessing services limit their power.  First is that they execute slower that the equivalent code written in .NET ArcObjects. For desktop tools that may not be a problem, but for the web this can negatively impacts user experience. 

Even if speed was not an issue, geoprocessing services are more cumbersome to move from one environment to another (e.g. staging -> production) because access to required source data is more tightly coupled with the tool.  Take the example of a server-side point clusterer.  To access the point feature class containing un-clustered points, a python script may use either of the following logic:

This first example uses a hard-coded absolute path that will need to be changed if the tool moves to an environment with a different directory structure.  A better choice would be the relative path example, but there is still hard-coding of the dataset's geodatabase and name.  With either choice, the tool suffers from a lack of true encapsulation which ends up making it harder to move and reuse on other point layers.

An aspect of the ArcGIS workflow for creating scripting tools which hinders a tool's flexibility is the need for explicit registration of input and output parameters in a properties dialog.

Python Cons

• Slower execution speed when compared to C#
• Tighter coupling with source data
• Double managing of service input/output parameters

When using geoprocessing services, the developer gains development time by sacrificing execution speed, and the fine-grain control offered by ArcObjects.

C# Server Object Extensions (ArcObjects 10 .NET SDK)

C# Pros

• Faster code execution
• Loose coupling with tool source data
• Easy to deploy across environments

I began using .NET ArcObjects because python geoprocessing service felt sluggish.  The allure of greater speed justified learning a more complex language.  What I found was not only increased speed, but also easier deployment and lower maintenance.

The idea of the Server Object Extension is that instead of publishing a service on top of the server, you extend the server itself.  This extended functionality can then be enabled for any service along with other out-of-the-box capabilities (e.g. Export Map, Identify, Find, Query, Generate KML, WMS).

The beautiful thing about Server Object Extensions is that ArcObjects has hooks into the source data within a map service.  These hooks give C# access to datasets using map service layer indexes instead of directory paths allowing for looser coupling with the server environment. ESRI Web APIs also use map service layer indexes which makes client-server communication easier.










C# Cons

• Large and complex library
• Not Cross-Platform (won't run on AGS Java version).
• Steeper learning curve

Writing SOEs requires knowledge of a more complex library of tools accessed through a more complex programming language.  If a project does not have somebody with ArcObjects experience, diving into Server Object Extensions will be intimidating. In ArcObjects, finding the core methods for geoprocessing is easy, but getting data into a class which implements the correct interface for use as input in those methods can be difficult.

An example of this would be getting a geometry object, passed in as JSON from the client, into a feature class which can be used by a Spatial Analyst tool.  The solution I used took 80 lines of not fun code to write.  The image below is meant to give an impression of the complexity.

Another area of trouble is .NET and inability to run outside of a Window's environment.  Being that ESRI offer's a Java version of the ArcGIS Server, choosing C# means limiting a tool to Window's and the .NET version of ArcGIS Server.  This hasn't been a major issue for me yet, but I want as many options in hosting environments as possible.

Overall I have found the benefits of Server Object Extensions outweigh the ease of development offered by  Python when performance and reusablilty are important. Initially, I found the complexity of working in ArcObjects discouraging, but slowly I wrote function to deal this data type conversions and development got easier.  I now consider arcpy as a solution only on the desktop.

Where Camp DC

I had a bunch of fun at Where Camp DC's ignite talks.  The format was a 5 minute presentation of 20 auto-incrementing slides (15s per slide).  Somehow I gained a second chin and a lisp in the process...

Exception Handling in ArcGIS Python

When writing python scripts, it is very important to add some basic exception handling so you're not pulling your hair out when a script doesn't work.

Recently I was running a geometry simplification script I've been using for about two years. I had never had a problem with it, but suddenly I started ERROR 999999: Error executing function. This error is basically ArcGIS saying "I have no idea what happened". Man I was bummmin...no lines numbers...no idea why after so long I had a problem. Was it the data? Well, yes...it turned out the Shape field was named 'Shape_' instead of 'Shape'.

To find the error, I wrapped the script inside a try/except block and used two functions in my helper class:


def reportGPErrors(self, geoprocessor):
    '''
    reports any geoprocessing related errors
    '''
    gpmsgs = geoprocessor.GetMessage(0)
    gpmsgs += geoprocessor.GetMessages(2)
    geoprocessor.AddError(gpmsgs)
    print gpmsgs

def reportSysErrors(self, geoprocessor):
    '''
    reports python system related errors
    '''
    tb = sys.exc_info()[2]
    tbinfo = traceback.format_tb(tb)[0]
    pymsg = "PYTHON ERRORS:\nTraceback Info:\n"
    pymsg += tbinfo + "\n"
    pymsg += "Error Info:\n" + str(sys.exc_type) + ": "
    pymsg += str(sys.exc_value) + "\n"
    geoprocessor.AddError(pymsg)
    print pymsg


Both of the functions take a geoprocessor, which gives the ability to print to both the python console and the ArcGIS message window.

To implement these function in geoprocessing scripts, simply import helpers class, and place a call to these function in the try/except block:


from helpers import Helpers
oHelpers = Helpers()
gp = oHelpers.createGeoprocessor()

try:

    code...code...code

except:
    oHelpers.reportGPErrors(gp)
    oHelpers.reportSysErrors(gp)


Creating the Python Geoprocessor

I was recently hired to refactor about 20 geoprocessing services written in python. At the beginning of each script was about 5 lines of boiler plate to instantiate the geoprocessing object (gp), check out the appropriate extensions, and set the necessary environmental variables. Each script began with more or less:


import arcgisscripting
gp = arcgisscripting.create(9.3)
gp.overwriteoutput = 1
gp.CheckOutExtension("Spatial")
gp.cellSize = "0.0041666667"


Since this code was being repeated throughout all of the gp services, I simply turned it into a function in my helpers class. This removed about 100 lines of boiler plate and gave the ability to fine-tune the geoprocessor for all of the services. There are some additional debugging methods inside createGeoprocessor which also reside in the helpers class. I will do a future post on those basic error handling functions soon.


def createGeoprocessor(self):
    '''
    Returns a 9.3/9.3.1 Geoprocessor whose
    overwrite property is set to true
    and Spatial Analyst extension checked-out
    '''
    import arcgisscripting
    gp = arcgisscripting.create(9.3)
    gp.overwriteoutput = 1
    try:
        if(gp.CheckExtension("Spatial") == "Available"):
        gp.CheckOutExtension("Spatial")
    else:
        raise "LicenseError"
    except "LicenseError":
        self.toScreen(gp, Messages.SPATIAL_ANALYST_ERROR)
    except:
        self.reportGPErrors(gp)
        self.reportSysErrors(gp)
    return gp


Simple Weighted Average with Python

While I suspect there must already be a built in weighted average tool in some python library ( scipy / numpy), for a project I'm working on I just wrote my own.  I would love to hear what people think and if I can improve the function in any way.

def weightedAverage(valueWeightList):
        '''
        Takes a list of tuples (value, weight) and
        returns weighted average as calculated by
        Sum of all values * weights / Sum of all weights
        '''
        numerator = sum([v * w for v,w in valueWeightList])
        denominator = sum([w for v,w in valueWeightList])
        if(denominator != 0):
            return(float(numerator) / float(denominator))
        else:
            return None