Table to csv

On the ESRI support forum someone asked how to convert an attribute table to a csv file. This can easily be done with a Python script. First we import arcgisscripting and the CursorIterator and initialize the geoprocessing object.

import arcgisscripting, string
from cursoriterator import CursorIterator

gp = arcgisscripting.create()

The next step is to initialize the different variables for this script. The workspace is the path to the directory, file geodatabase, personal geodatabase, sde,... where the script can find the table or feature class that you want to convert. The table is the name of the table or feature class you want to export. The outputpath is the path to file where you want your export to be written to. Be aware that the output directory should exist. If the file doesn't exist it will be created otherwise it will be overwritten. The csvseparator is the column separator for the output file. It could be anything you want. If you want to use tabs you should set this variable to '\t'. The last variable is ignorefields. If a field name in the table or featureclass is equal to a string in this list then this field won't be exported. This can be used for excluding the shape field, the objectid field or other fields you don't want to appear in the output file.

gp.workspace = r'' ## workspace of the table or feature class
table = '' ## table or feature class from wich the attributes should be exported
outputpath = r'' ## path to the file where the output should be written to
csvseparator = ',' ## column separator field
ignorefields = [] ##list with fields to ignore

Then we define a function wich purpose is to create a list of fieldnames that are to be exported to the output file. This is done by iterating over the fields cursor that can be accessed by the gp function listfields. The fields that are in the ignorefields list won't be returned.

def get_fieldnames(fields, ignorefields=[]):
   ignorefieldsupper = map(string.upper, ignorefields)
   fields = CursorIterator(fields)
   fields_output = []
   for field in fields:
       if not field.name.upper() in ignorefieldsupper:
           fields_output.append(field.name)
   return fields_output


fields = gp.listfields(table)
fieldnames = get_fieldnames(fields)

Now we are going to iterate over the rows in the table or featureclass and create a list with a string representation of each row. For each field we add the value of a field to a list. Then we concactenate this values with the chosen separator and add this to the output list.

rows = gp.searchcursor(table)
rows = CursorIterator(rows)

output = []
output.append(csvseparator.join(fieldnames))
for row in rows:
   outputrow = []
   for fieldname in fieldnames:
       outputrow.append(str(row.getvalue(fieldname)))

   outputrow = csvseparator.join(outputrow)
   output.append(outputrow)

Finally we write the outputlist to the defined outputpath.

f = open(outputpath, 'w')
f.write('\n'.join(output))
f.close()

The full code can be found on the arcscripts site.

Projections and transformation

On his blog Richie Carmichael (First edition Second edition) showed how to use the ESRI Projection Engine from .NET. In this article I am going to do the same with Python.

My original problem was that I needed to convert some coordinates from ED50 Projected and ED50 geographic to ETRS89 projected for the UTM zone 31.

First of all we will need the pe.dll provided by ESRI which can be found in the BIN folder of the ArcGIS installation. Copy this dll in the same folder as your Python script.

To use this dll from Python you can use the ctypes package. This extension started has a separate project in Python 2.3 but it has been fully incorporated in Python 2.5. So if you're still using Python 2.4 that shipped with ArcGIS 9.2 You will have to download and install ctypes. This extensions creates the possibility to easily use C dll's from within Python.

Now lets start writing some code. First we import the ctypes package.

    import ctypes

The pe.dll is accessed in the following way.

    pe = ctypes.windll.pe

To transform from one projection to another we need to first create a projection object for the input spatial reference and the output spatial reference. Then we also need a geographic transformation object. We do this by calling the appropriate factory functions which can be found in the SDE C API documentation on EDN.

    input_proj = pe.pe_factory_projcs(input_proj_code)
    transformation = pe.pe_factory_geogtran(transformation_code)
    output_proj = pe.pe_factory_projcs(output_proj_code)

To lookup the code of your projection you will need to take a look at the following sites page 1 page 2 page 3 and page 4. The code of your geographic transformation can be found on the this sites page 1 page 2 and page 3.

To transform points with the projection engine we need to pass the points as a 2 dimensional array of points. In order to do this I created a small Python point class and methods to convert from a list of points to a 2 dimensional C array and back to a list of points. If you want more information about this I suggest you take a look at the ctypes documentation (tutorial reference) and also have a look at the needed datatypes for the pe functions.

    class Point:
        def __init__(self, x, y):
            self.x = x
            self.y = y
        def __str__(self):
            return 'x : ' + str(self.x) + ' y : ' + str(self.y)
    
    def _ctypes_to_points(point_array):
        points = []
        for point in point_array:
            point = Point(point[0], point[1])
            points.append(point)
        return points

    def _points_to_ctypes(points):
        array_type = ctypes.ARRAY(self.point_type, len(points))
        ctype_points = []
        for point in points:
            ctype_point = self.point_type(point.x, point.y)
            ctype_points.append(ctype_point)
        return array_type(*ctype_points)

Now we can create the function for doing the actual transformation from input projected reference system to the output projected reference system. As you can see we first need to convert the projected coordinates to geographic coordinates. Then transform these geographic coordinates with the give transformation. Finally we convert the geographic coordinates back to projected coordinates.

def transform(self, points):
    points_count = len(points)
    point_array = self._points_to_ctypes(points)
    pe.pe_proj_to_geog(self.input_proj, points_count, 
                       point_array)
    pe.pe_geog1_to_geog2(self.transformation, points_count, 
                         point_array, None)
    pe.pe_geog_to_proj(self.output_proj, points_count, 
                       point_array)
    return self._ctypes_to_points(point_array)

The complete code looks like this. I've put the projection related code in a class and at the bottom I added some test code.

import ctypes

class Point:
    def __init__(self, x, y):
        self.x = x
        self.y = y
    def __str__(self):
        return 'x : ' + str(self.x) + ' y : ' + str(self.y)

class Proj:
    def __init__(self, input_proj_code, transformation_code, 
        output_proj_code):
        self.pe = ctypes.windll.pe
        self.input_proj = self.pe.pe_factory_projcs(input_proj_code)
        self.output_proj = self.pe.pe_factory_projcs(output_proj_code)
        self.transformation = self.pe.pe_factory_geogtran(transformation_code)
        self.point_type = ctypes.c_double * 2

    def _points_to_ctypes(self, points):
        array_type = ctypes.ARRAY(self.point_type, len(points))
        ctype_points = []
        for point in points:
            ctype_point = self.point_type(point.x, point.y)
            ctype_points.append(ctype_point)
        return array_type(*ctype_points)

    def _ctypes_to_points(self, point_array):
        points = []
        for point in point_array:
            point = Point(point[0], point[1])
            points.append(point)
        return points

    def transform(self, points):
        points_count = len(points)
        point_array = self._points_to_ctypes(points)
        self.pe.pe_proj_to_geog(self.input_proj, points_count, 
                                point_array)
        self.pe.pe_geog1_to_geog2(self.transformation, 
                                  points_count, point_array, None)
        self.pe.pe_geog_to_proj(self.output_proj, points_count, 
                                point_array)
        return self._ctypes_to_points(point_array)

if __name__ == '__main__':
    ed50_proj_code = 23031
    ed50_etrs89_transformation = 1650
    etrs89_proj_code = 25831
    proj = Proj(ed50_proj_code, ed50_etrs89_transformation, etrs89_proj_code)
    points = [Point(585635, 5690770), Point(585654, 5690110)]
    points = proj.transform(points)
    for point in points:
        print point

I hope you liked this post. Have fun.

Using for-loops for cursors

When writing Python scripts to automate ArcGIS you often encounter all kind of cursors. Cursors are returned mostly returned by the gp object when calling functions like listtables, listfeatureclasses, listrasters, listfields, searchcursor, updatecursor, ... Usually cursors are used in the following way :

    import arcgisscripting
    gp = arcgisscripting.create()
    inputFC = r'path to a feature class'

    inRows = gp.searchcursor(inputFC)
    inRow = inRows.next()
    while inRow:
        ...do something
    inRow = inRows.next()

This is quiet a verbose way to loop through the rows. Fortunately there is a way to turn the cursor into an iterator that can be used in a for loop. This can be done by adding the following short class to your code.

    class CursorIterator(object):
        def __init__(self, cursor):
            self.cursor = cursor
            self.cursor.reset()
        def next(self):
            n = self.cursor.next()
            if n is None:
                self.cursor.reset()
                raise StopIteration
            return n
        def __iter__(self):
            return self

As it is nicely explained in Expert Python Programming by Tarek Ziadé we just created an iterator class. This is a container class with two methods namely next and __iter__. Usage of the CursorIterator class looks like this.

    rows = gp.searchcursor(inputFC)
    rowsIterator = CursorIterator(rows)
    for row in rowsIterator:
        ... do something

or shorter

    for row in CursorIterator(gp.searchcursor(inputFC)):
        ... do something

I hope you enjoyed my first post. Have fun.

Update : As Jason Schreider points out in his comment the CursorIterator class can be easily replaced by the following :

    rows = gp.searchcursor(inputFC)
    for row in iter(rows.next, None):
        ... do something