bcdc spatial

Interesting Tim Peters interview...

2013-05-24T18:59:00.001-07:00

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

2013-03-25T23:06:00.000-07:00

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

2013-03-20T14:38:00.001-07:00

This a great resource for learning to Sublime Text 2:

Converting ArcGIS Feature Class to PyTables

2013-02-27T08:35:00.000-08:00

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

2013-01-09T18:51:00.000-08:00

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

2012-09-23T20:58:00.001-07:00

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

2012-09-18T22:26:00.000-07:00

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

2012-05-07T20:33:00.000-07:00

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

2012-02-13T21:03:00.001-08:00

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

2012-01-29T14:36:00.001-08:00

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.

Some random ArcObjects that make writing Server Object Extensions easier...(C#)

2011-12-11T14:00:00.000-08:00

Server Object Extensions are one of the most powerful features of ArcGIS Server. To make writing them a bit easier for people, I am posting several of the functions I commonly use to deal with incoming data, execute an operation, and provide a response. As an example, I will take the scenario of doing a zonal statistics operation based on a user digitized polygon. The idea is that the user selects a variable for analysis (elevation, temperature, precipitation, etc.), draws a polygon in Flex/JS/Silverlight, the polygon is sent to the server, and zonal statistics for that polygon/variable are returned.

1. JSON to IGeometry: The server will receive a json object representing the polygon digitized by the user. The first step is to convert the JSONObject to IGeometry. I particularly liked NicoGIS's solution to the issue which is shown below:

public IGeometry ConvertAnyJsonGeometry(JsonObject jsonObjectGeometry)
        {
            object[] objArray;

            if (jsonObjectGeometry.TryGetArray("rings", out objArray))
            {
                return Conversion.ToGeometry(jsonObjectGeometry, esriGeometryType.esriGeometryPolygon);
            }

            if (jsonObjectGeometry.TryGetArray("paths", out objArray))
            {
                return Conversion.ToGeometry(jsonObjectGeometry, esriGeometryType.esriGeometryPolyline);
            }

            if (jsonObjectGeometry.TryGetArray("points", out objArray))
            {
                return Conversion.ToGeometry(jsonObjectGeometry, esriGeometryType.esriGeometryMultipoint);
            }

            try
            {
                return Conversion.ToGeometry(jsonObjectGeometry, esriGeometryType.esriGeometryPoint);
            }
            catch
            {
                try
                {
                    return Conversion.ToGeometry(jsonObjectGeometry, esriGeometryType.esriGeometryEnvelope);
                }
                catch
                {
                    return null;
                }
            }
        }

2. IGeometry to IFeatureClass: Once I have IGeometry, many times I want to convert it to an IFeatureClass so I can use it in an operation like zonal statistics. Zonal Statistics actually takes an IGeodataset, but IFeatureClass extends IGeodatatset, so all that is required is a cast. Here is a function to convert IGeometry to IFeatureClass:

public IFeatureClass CreateFeatureClassFromGeometry(IGeometry pGeometry, IFeatureWorkspace pOutFeatWorkspace, int wkid = 4236)
        {
            try
            {
                IFields pFields = new Fields() as IFields;
                {
                    // Set up the shape field for the feature class
                    IFieldsEdit pFieldsEdit = (IFieldsEdit)pFields;
                    IField pField = new Field();
                    IFieldEdit pFieldEdit = (IFieldEdit)pField;
                    pFieldEdit.Name_2 = "Shape";
                    pFieldEdit.Type_2 = esriFieldType.esriFieldTypeGeometry;

                    IGeometryDef pGeometryDef = new GeometryDef();
                    IGeometryDefEdit pGeometryDefEdit = (IGeometryDefEdit)pGeometryDef;
                    pGeometryDefEdit.GeometryType_2 = pGeometry.GeometryType;

                    ISpatialReference pSpatialReference;
                    if (wkid == 4326)
                    {
                        ISpatialReferenceFactory2 pSpaRefFact2 = new SpatialReferenceEnvironment() as ISpatialReferenceFactory2;
                        IGeographicCoordinateSystem pGeoCoordSys = pSpaRefFact2.CreateGeographicCoordinateSystem(wkid);
                        pSpatialReference = (ISpatialReference)pGeoCoordSys;
                    }

                    else if (wkid == 102100 || wkid == 3857)
                    {
                        ISpatialReferenceFactory2 pSpaRefFact2 = new SpatialReferenceEnvironment() as ISpatialReferenceFactory2;
                        IProjectedCoordinateSystem pProjCoordSys = pSpaRefFact2.CreateProjectedCoordinateSystem(wkid);
                        pSpatialReference = (ISpatialReference)pProjCoordSys;
                    }

                    else
                    {
                        throw new ArgumentNullException("Invalid Spatial Reference Well Known Id: Please use 4326 for Geographic or 102100 for Web Mercator");
                    }

                    ISpatialReferenceResolution pSpatialReferenceResolution = (ISpatialReferenceResolution)pSpatialReference;
                    pSpatialReferenceResolution.ConstructFromHorizon();

                    pGeometryDefEdit.SpatialReference_2 = pSpatialReference;
                    pFieldEdit.GeometryDef_2 = pGeometryDef;
                    pFieldsEdit.AddField(pField);

                    // Add other required fields to the feature class

                    IObjectClassDescription pObjectClassDescription = new FeatureClassDescription();

                    for (int i = 0; i < pObjectClassDescription.RequiredFields.FieldCount; i++)
                    {
                        pField = pObjectClassDescription.RequiredFields.get_Field(i);
                        if (pFieldsEdit.FindField(pField.Name) == -1)
                            pFieldsEdit.AddField(pField);
                    }
                }

                // Create the feature class
                string sFeatureClassName = "tmp" + Guid.NewGuid().ToString("N");
                IFeatureClass pFeatureClass = pOutFeatWorkspace.CreateFeatureClass(
                   sFeatureClassName, pFields, null, null, esriFeatureType.esriFTSimple, "Shape", null);

                // Add the input geometry to the feature class and return
                IFeatureCursor pFeatureCursor = pFeatureClass.Insert(true);
                IFeatureBuffer pFeatureBuffer = pFeatureClass.CreateFeatureBuffer();

                pFeatureBuffer.Shape = pGeometry;
                pFeatureCursor.InsertFeature(pFeatureBuffer);

                // Flush the feature cursor and release COM objects
                pFeatureCursor.Flush();
                Marshal.ReleaseComObject(pFeatureBuffer);
                Marshal.ReleaseComObject(pFeatureCursor);
                return pFeatureClass;
            }
            catch (Exception e)
            {
                Console.WriteLine(e.Message);
                throw e;
            }
        }

3. Create In-Memory Workspace: As you notice, the function above requires a workspace. I like to use in-memory workspaces and to create them I use the following function:

        public IWorkspace CreateInMemoryWorkspace()
        {
            try
            {
                // Create an InMemory workspace factory.
                IWorkspaceFactory workspaceFactory = new InMemoryWorkspaceFactory() as IWorkspaceFactory;

                // Create an InMemory geodatabase.
                IWorkspaceName workspaceName = workspaceFactory.Create("", "MyWorkspace", null, 0);

                // Cast for IName.
                IName name = (IName)workspaceName;

                //Open a reference to the InMemory workspace through the name object.
                IWorkspace workspace = (IWorkspace)name.Open();
                return workspace;
            }

            catch (Exception e)
            {
                Console.WriteLine(e.Message);
                throw e;
            }
        }

4. Get Raster Dataset from Map Service: Keep in mind that an SOE has access to any of the feature classes in the service which uses the SOE as a capability. This means that you shouldn't usually need to have any hardcoded data paths in the actual code. Rather you should try an only use data from the map service itself. If you send in the index of the value raster the user wants to run analysis on, look how easy it is to open:

public IGeoDataset GetGeoDatasetByMapServiceIndex(IMapServer3 mapServer, int layerID)
        {
            IMapServerDataAccess dataAccess = (IMapServerDataAccess)mapServer;
            return dataAccess.GetDataSource(mapServer.DefaultMapName, layerID) as IGeoDataset;
        }

5. Pre-loading Rasters during Init(): Loading and unloading rasters from memory can be an expensive operation. Fortunately with SOEs, you can load data into memory when the service starts up, and have it standing by for any request which come in. This is different than Python geoprocessing services which basically need to load data once the actual request is receive. One idea I've been using is to create a data dictionary when the service starts up which keeps my value rasters standing by:

public Dictionary<int, IGeoDataset> CreateGeodatasetDictionary(IMapServer3 mapServer)
        {
            Dictionary<int, IGeoDataset> data_dictionary = new Dictionary<int, IGeoDataset>();

            IMapServerInfo msInfo = mapServer.GetServerInfo(mapServer.DefaultMapName);
            IMapLayerInfos layerInfos = msInfo.MapLayerInfos;
            int c = layerInfos.Count;

            for (int i = 0; i < c; i++)
            {
                IGeoDataset data = GetGeoDatasetByMapServiceIndex(mapServer, c);
                data_dictionary.Add(c, data);
            }

            return data_dictionary;
        }

Basically this dictionary is a mapping between map service indexes and their IGeodatasets which I create in the SOE Init() function and store in memory for the life of the SOE.

6. Running Zonal Statistics: Now that we have both the zone dataset and the value raster at our figure tips, the zonal statistics operation is just a couple lines of code. There are many examples of how to use the zonalops class, but here's the idea:

IZonalOp oZonalOp = new RasterZonalOpClass();
ITable zonalStats = oZonalOp.ZonalStatisticsAsTable(zones, valueRaster, true);

7. ITable to IRecordSet: To return the results of the Zonal Stats from the SOE, you must first convert the ITable to an IRecordSet. Here an straightforward way to do it:

public static IRecordSet ConvertTableToRecordset(ITable table)
        {
            IRecordSetInit recordSetInit = new ESRI.ArcGIS.Geodatabase.RecordSetClass();
            recordSetInit.SetSourceTable(table, new QueryFilterClass());
            IRecordSet recordset = recordSetInit as IRecordSet;
            return recordset;
        }

8. Serialize IRecordset: To convert the IRecordset into return JSON, you can simply use the Conversion.ToJSON method as follows:

IRecordSet return_records = BSharpUtilities.ConvertTableToRecordset(areaTable);
byte[] jsonBytes = Conversion.ToJson(return_records);

return Encoding.UTF8.GetString(jsonBytes, 0, jsonBytes.Length);

That's about it for now. For information on getting started with SOEs, I would check out ESRI's samples that come with the .NET SDK and also NicoGIS posts

Simple Source Code Count (Python)

2011-11-14T20:06:00.000-08:00

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)

Managing MapMouseEvents: Double-click zoom vs. showing infowindow (ESRI Flex API 2.4)

2011-10-13T18:41:00.000-07:00

Apps built with the ESRI Flex API have the default behavior of zooming-in when you double-click on the map. This is a great feature, but the double-click and single-click dispatch the same event may confuse your handler functions.

When you click once on the map, MapMouseEvent.MAP_CLICK is dispatched.
When you double-click, MapMouseEvent.MAP_CLICK is dispatched twice before zooming begins.

Here is a great example of the events which are dispatched with the Flex API.

If you are listening to MapMouseEvent.MAP_CLICK to show an infowindow, the infowindow will open during the double-click zoom which I consider an undesirable behavior.

You can see the behavior I am trying to fix here.

The fix which I am using involves using a timer in the single-click handler to wait and see if a second MAP_CLICK event is dispatched. If so, I assume that a double-click zoom is going on and I don't let the infowindow open.

Here is the code from the example above modified so the infowindow does not appear during a double-click:

<?xml version="1.0" encoding="utf-8"?>
<s:Application xmlns:fx="http://ns.adobe.com/mxml/2009"
               xmlns:s="library://ns.adobe.com/flex/spark"
               xmlns:esri="http://www.esri.com/2008/ags"
               pageTitle="MapClick - click map to get current location">
    <fx:Script>
        <![CDATA[
            import com.esri.ags.events.MapMouseEvent;
            import com.esri.ags.geometry.MapPoint;
            import com.esri.ags.utils.WebMercatorUtil;
            
            /**
             * This whole timer business here is to deal with the
             * double-click zoom vs. the single click info window issue.
             * 
             * If another MAP_CLICK event is heard within the alloted time
             * then the infowindow will not open.
             */ 
            private const eventTimer:Timer = new Timer(300, 1);
            private function onMapClick(event:MapMouseEvent):void 
            {
                function onTimerComplete(timerEvent:TimerEvent):void
                {
                    this.eventTimer.removeEventListener(TimerEvent.TIMER_COMPLETE, onTimerComplete);
                    var latlong:MapPoint = WebMercatorUtil.webMercatorToGeographic(event.mapPoint) as MapPoint;
                    myMap.infoWindow.label = "You clicked at "
                        + event.mapPoint.x.toFixed(1) + " / " + event.mapPoint.y.toFixed(1)
                        + "\nLat/Long is: " + latlong.y.toFixed(6)
                        + " / " + latlong.x.toFixed(6);
                    myMap.infoWindow.show(event.mapPoint);
                }
                
                //If timer is already running then don't show the infowindow
                if(eventTimer.running)
                {
                    this.eventTimer.removeEventListener(TimerEvent.TIMER_COMPLETE, onTimerComplete);
                    this.eventTimer.stop();
                    return;
                }
                
                eventTimer.addEventListener(TimerEvent.TIMER_COMPLETE, onTimerComplete);
                eventTimer.start();
            }
        ]]>
    </fx:Script>
    
    <esri:Map id="myMap"
              mapClick="onMapClick(event)"
              scale="50000000">
        <esri:ArcGISTiledMapServiceLayer url="http://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer"/>
    </esri:Map>
</s:Application>

I would much rather not use a timer to achieve this. If anybody has a better way, please let me know with a comment.

Sorting an array on multiple fields of a nested object...

2011-09-15T19:46:00.000-07:00

I recently ran into the need to sort an array of AS3 objects by multiple fields of a nested object. In JSON, the AS3 would look sort of like this:

[

    {"feature":
        {"attributes":{"country":"Brazil", "state":"Bahia"}
        {"geometry":.....}
    },
    {"feature":
        {"attributes":{"country":"Brazil", "state":"Sao Paulo"}
        {"geometry":.....}
    }
]

Ok, so the example above is already sorted, but imagine this sort of array and needing to sort alphabetically descending first by country and then by system. I'm pretty sure that the AS3 sortOn() only handles first-level properties for each object in the array. I chose to write a custom function which would use array.sort():

private function nestedTwoColumnSort(array:Array, parentField:String, 
                                     sortField1:String, sortField2:String):void
{
    function multipleSort(a:Object, b:Object):Object
    {
        var w:String = a[parentField][sortField1].toLowerCase();
        var x:String = b[parentField][sortField1].toLowerCase();
        
        //sort first field
        if(w < x)
        {
            return -1;
        }
        
        if(w > x)
        {
            return 1;
        }
        
        //if first field is the same, then sort on second field
        if(w == x)
        {
            var y:String = a[parentField][sortField2].toLowerCase();
            var z:String = b[parentField][sortField2].toLowerCase();
            
            if(y < z)
            {
                return -1;
            }
            
            if(y > z)
            {
                return 1;
            }
        }
        
        return 0;
    }
    
    //Sort once for the first field
    array.sort(multipleSort);
    
    //Sort once for second field
    array.sort(multipleSort);
}

The nestedTwoColumnSort function takes the array, the name of the nested object, and the two fields you wish to sort by. It then runs the sort once to fully sort the first field, then again banking on all first fields to be the same.

Here's the actual function where I call the sort function:

private function onFeaturesLoadComplete(event:AppEvent):void
{
    var features:Array = event.data as Array;
    this.nestedTwoColumnSort(features, 'attributes', 'country', 'system');

    blah blah blah....

There very well might be a better way to do this, but this was pretty straightforward to implement and I'm not seeing any noticeable performance issues in the application (but it also only runs once in the app...).

I would love to hear how others accomplish the same task.

Singletonitis

2011-08-16T19:46:00.000-07:00

I am a self-loathing user of the Singleton design pattern. Global states make small scale development much less painful, especially when maintaining the state of a user interface in Flex.

Misko Hevery gives a great explanation of why the Singleton design pattern is really an anti-pattern. I am continuing to use Singletons for now, but I think it's time to learn Robotlegs or another dependency injection framework for ActionScript...

Cannot Log into REST Admin (AGS 10)

2011-08-03T11:19:00.000-07:00

I was dealing with an annoying issue today involving logging into the REST admin area of ArcGIS Server 10. Basically I would log in, nothing would happen and no error message would be displayed.

For anybody having this issue, it turns out it was related to the domain of the user's account. Basically if the user account is not local, but instead on a different domain, they will not be able to log into the ArcGIS Server 10 REST admin interface. This is a known bug #NIM059442.

Upgrade to the latest service pack and it should fix the problem.

Metro DC Flash User's Group

2011-07-24T19:00:00.000-07:00

This weekend I had the pleasure of attending a meeting of the Metro DC Flash User's Group. The meeting was held at Center for Digital Imaging Arts of Boston University (satellite location) in Georgetown, DC.

Doug Chaplow first presented on developing for mobile devices with Flash. One of the most important points I took away from his presentation was that AIR for iOS (not really AIR) and AIR for Android now have basically all the same capabilities. For a while, Adobe had been focused on Android during the political turmoil with Apple. Once Apple removed restricts on the cross-compiled apps, Adobe began work again on the cross-compiler for iOS and its now supposedly just as robust as AIR for Android.

One gentlemen in the group also mentioned that if you are writing an app for both iOS (Apple), do not include the word "Android" or "Blackberry" in your code or else it will get rejected by Apple.

Example of what NOT to do in your flash mobile app:

switch(OS_NAME)
{
    case 'Android':
        result = 'Get Rejected From Apple Store';
        break;
    default:
        result = 'Get Rejected From Apple Store';
        break;
}

Can anybody confirm if this is true? I wonder if this also applies to comments like:

//Developing for Android is better than for iOS

Don Anderson gave the second presentation on his MVC framework called Rivet. While the framework is still in development, it looks like the key goal was event listener management. In Rivet, it looked like all event listeners for a view are added and removed in a separate controller when a view is activated/deactivated. The architecture was definitely more pure ActionScript focused than Flex. Overall, the presentation was very good and I always welcome any talk on a variation of the Model-View-Controller architecture.

I would really urge anybody in the DC metro area to check out this meetup group. It was a great time and everybody was very friendly and knowledgeable.

import this

2011-07-21T08:24:00.001-07:00

>>> 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!

Setting Feature Layer maxAllowableOffset property

2011-07-18T14:32:00.000-07:00

Derek Swingley from ESRI recently enlightened me on the use of feature layers for the ArcGIS APIs for Flex, JavaScript, and Silverlight.

Max Allowable Offset basically determines the level of geometry generalization when returned from the mapping server. The larger the offset, the greater the simplification. To calculate the appropriate maxAllowableOffset, I am simply using the ground resolution of the current zoom level of the map (AS3):

private const EARTH_CIRCUM:Number = 2 * 3.14159265 * 6378137;
private function computeAllowableOffset():void
{
    fLayer.maxAllowableOffset = (EARTH_CIRCUM / (256 * Math.pow(2, oMap.level)));
}

If you are using a WGS84 Basemap, you can just divide the result above by the number of meters per degree:

private const EARTH_CIRCUM:Number = 2 * 3.14159265 * 6378137;
private function computeAllowableOffset():void
{
    fLayer.maxAllowableOffset = (EARTH_CIRCUM / (256 * Math.pow(2, oMap.level))) / (EARTH_CIRCUM / 360);
}

These functions make it so that at the equator there should only be one vertex per pixel. I've been wanting a little more stylized features, so I've been adding about 35 - 40% additional offset:

private const EARTH_CIRCUM:Number = 2 * 3.14159265 * 6378137;
private function computeAllowableOffset():void
{
    fLayer.maxAllowableOffset = (EARTH_CIRCUM / (256 * Math.pow(2, oMap.level))) * 1.35;
}

For more information on feature layers, check out Derek's blog posts

Great REST-based Server Object Extension posts from NicoGIS

2011-07-02T12:16:00.000-07:00

NicoGIS is an awesome GIS blog. Over the past year, I've been learning REST-based Server Object Extensions for ArcGIS Server 10. NicoGIS has some one of the best posts on the subject. I ended up reading the post in Italian for fun, but I bet Google can translate it.

ESRI Developer Summit Python Presentation

2011-06-26T20:21:00.000-07:00

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)

2011-06-26T20:13:00.000-07:00

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

2011-06-26T15:01:00.000-07:00

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...

Find Watershed SOE

2010-09-13T22:36:00.000-07:00

The Find Watershed Server Object Extension (SOE) plugs into ArcGIS Server 10 and adds functionality for generating watersheds over REST. The SOE takes two parameters: hydroshed id and a pour point location. Proper flow accumulation and flow direction rasters are retrieved by hydroshed id from a file geodatabase. The location parameter is the pour point. The goal of the tool is to generate a polygon representing the area of water which flows into a point.

I like the motto, "if you can't do in on desktop, you won't be able to do it on server." Keeping this idea in mind, I first wrote the tool as a add-in for ArcMap. This allowed me to debug core business logic for the tool without the complications of the server.

The main difference between a 'desktop tool' and Server Object Extension is the ability to deserialize requests and serialize responses over http. Since I was new to C# and Server Object Extensions, I first published a sample SOE called SimpleRESTSOE from the samples folder of the ArcGIS 10 install directory. This tool was simply an echo service (i.e. pass a string in, and get the string back). Requests declare 'text' and 'f' (output format) parameters as part of a url request:

/echo?text=helloWorld&f=json

and they SOE responds:

{"text":"helloWorld"}

I got the sample SOE registered with Windows using the command:

RegAsm SimpleRESTSOE.dll /codebase

I then registered the SOE with ArcCatalog and the Server Manager using the companion file located in the samples folder (RegisterSimpleRESTSOE.dll).

Once I verified the sample SOE was running, I simply started modifing the rest schema to take my custom parameters. To do this, I modified the CreateRestSchema method of the COM class for the tool. I also added some additional variables to support my business logic:

After modifying the CreateRestSchema method the tool was wired to take the following request:

/echo?hydroshed_id=sa7&location={"x"=-54,"y"=-24}&f=json

I wrapped my business logic in a function and called it from the REST operation's handler EchoInput. THe EchoInput method parses the incoming request, calculates the watershed, and serialize the output.

There is a block of 3 lines towards the bottom that comprise the business logic from my ArcMap tool. The first two are helper functions (CreateInMemoryWorkspace / CreateFCFromGeometry), and the third is the actual watershed operation:

Global Hydrosheds Lookup

I built a small client application to test Find Watershed over the web. In the app, the user first supplies the location by clicking on the map. The click location is passes to the server into a query against the lookup map to determine the hydroshed id then :

building a client

Site Selection with ArcObjects

2010-08-25T23:13:00.001-07:00

In continuing my application development class through Penn State, I wrote a small site selection tool using C# and ArcObjects. The interesting part of that the application was stand-alone in contrast to the ArcMap tools we wrote in weeks past. I made it through the practice exercises no problem, but ran into some issues when writing the deliverable for the week. Here are some of those problems and solutions:

Referencing Map Layers
One of the bugs in my app which gave me the most trouble was related the order in which layers are loaded into the application. The problem was caused by the following code:



            this.axMapControl1.AddLayerFromFile(dataPath + "counties.lyr");

            pCntyLayer = this.axMapControl1.get_Layer(0) as IFeatureLayer;

            pCntyFC = pCntyLayer.FeatureClass;



            this.axMapControl1.AddLayerFromFile(dataPath + "recareas.lyr");

            pRecLayer = this.axMapControl1.get_Layer(1) as IFeatureLayer;

            pRecFC = pRecLayer.FeatureClass;



            this.axMapControl1.AddLayerFromFile(dataPath + "interstates.lyr");

            pInterLayer = this.axMapControl1.get_Layer(2) as IFeatureLayer;

            pInterFC = pInterLayer.FeatureClass;



            this.axMapControl1.AddLayerFromFile(dataPath + "cities.lyr");

            pCityLayer = this.axMapControl1.get_Layer(3) as IFeatureLayer;

            pCityFC = pCityLayer.FeatureClass;

Notice the index values in the get_Layer function. I initially expected the layer indexes to be the same as the order in which they were added to the map. This turned out to be wrong. Actually each layer was added to the beginning of the array, meaning that for the code above, each of the indexes should be 0.

This bug through me on a wild goose chance because it turned out I got the first layer index correct, so I didn't get a error message until I was further into writing application logic. In the code above, each of the layers will be the county layer. By coincidence this was the first layer I was doing operations one which cause additional confusion when operations on the cities layer didn't work.

Creating the Geometry Bag
When creating a geometry bag of features from Query Filter, I first attempted writting my own functions but quick found that the professors solution was more elegant. I translated the VBA code into the following C# function:



private IGeometryBag getGeomBag(IQueryFilter pQueryFilter, IFeatureClass pFClass)

        {

            IEnumGeometry pEnumGeom = new EnumFeatureGeometryClass();

            IEnumGeometryBind pEnumGeomBind = pEnumGeom as IEnumGeometryBind;

            pEnumGeomBind.BindGeometrySource(pQueryFilter, pFClass);

            IGeometryFactory pGeomFactory = new GeometryEnvironmentClass();

            return pGeomFactory.CreateGeometryFromEnumerator(pEnumGeom) as IGeometryBag;

        }

This function applies the filter and return a nice geometry for additional site selection steps.

Conclusion
While I got my stand alone application working, I was unable to integrate it into excel and powerpoint. I am currently experimenting with VSTO (Visual Studio Tools for Office) to see if I find a solution.