GIS Solved

Visualizing Geometries

2012-01-25T20:35:00.000+01:00

In this post I'm going to give some more information about the GeometryVisualizer that I created recently.

The goal of the GeometryVisualizer is to be able to quickly visualize a GeoJson, WKT or ESRI JSON geometry. I created this to be able to quickly see how some ESRI JSON and WKT geometries looked like and visually verify some geometry conversion code I'm writing for my GeometryService in node.js. You can also use this to visualize the WKT or GeoJSON output of a PostGIS geometry query.

The ESRI JSON geometries are rendered with the ArcGIS JavaScript API. The WKT and GeoJSON geometries are rendered with OpenLayers.

Some code

In this section I'll show some parts of the code that made this web application possible and compare the ArcGIS JavaScript API with the OpenLayers API. First we'll take a look at the initialization of the map. The most notable facts about the initialization are that we need to add a layer before we are able to add graphics and that disabling the zoomwheel is less straight forward with the OpenLayers API.

ESRI JavaScript API map initialization:

var map = new esri.Map("esriJsonMap");
map.disableScrollWheelZoom();

// add a layer and directly hide it because otherwise we can't create graphics
var basemapURL= "http://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer"

var basemap = new esri.layers.ArcGISDynamicMapServiceLayer(basemapURL);
map.addLayer(basemap);
basemap.visible = false;

OpenLayers map initialization:

var map = new OpenLayers.Map("openLayersMap");

var controls = map.getControlsByClass('OpenLayers.Control.Navigation'); 
for(var i = 0; i<controls.length; ++i){ 
    controls[i].disableZoomWheel(); 
}

var layer = new OpenLayers.Layer.ArcGIS93Rest( "ArcGIS World Street Map", 
                "http://server.arcgisonline.com/ArcGIS/rest/services/World_Street_Map/MapServer",
                {layers: '0'} );
map.addLayer(layer);
layer.setVisibility(false);

To be able to render an ESRI JSON geometry string with the ArcGIS JavaScript API you have to do some specific steps :

Clear the graphics
Parse the geometry string representation
Create a graphic and set its symbol according to the geometry type.
Update the extent of the map

Not that for updating the extent, I had to access an undocumented property of the graphics object and for point graphics the extent had to be enlarged.

function addGeometry(geometryText) {
    
    if(map.graphics === null) return; // map not loaded yet or something else is wrong
            
    map.graphics.clear();

    var geometryJson = JSON.parse(geometryText);
    var graphic = new esri.Graphic(geometryJson);

    if (!graphic.symbol){
      switch(graphic.geometry.type){
        case "point":
        case "multipoint":
          graphic.setSymbol(new esri.symbol.SimpleMarkerSymbol());
          break;
        case "polyline":
          graphic.setSymbol(new esri.symbol.SimpleLineSymbol());
          break;
        default:
          graphic.setSymbol(new esri.symbol.SimpleFillSymbol());
          break;
      }
    }

    map.graphics.add(graphic);
    // use undocumented object
    var ext = graphic._extent;
    // if point expand the extent
    if(ext.getWidth() < 0.00000001 || ext.getHeight() < 0.00000001){
      var factor = 1;
      ext.update(ext.xmin - factor, ext.ymin - factor, ext.xmax + factor, ext.ymax + factor, ext.spatialReference);
    }

    map.setExtent(ext.expand(2));
}

The process of adding WKT and GeoJSON geometries was very similar.

Remove all features
Parse the geometry string representation with a predefined formatter (OpenLayers.Format.WKT or OpenLayers.Format.GeoJSON)
Update the extent of the map.

Note that I didn't had to specify the symbology and zooming to the extent of the feature was easier. The WKT parsing was very straight forward. For the GeoJSON parsing I had to add two things. First I had to ensure that the parsed JSON had a property "type" with as value "Feature" and the parsed object returned a feature collection with one feature instead of directly returning the feature.

function addGeometry(geometryInput, parser, type){
  vectorLayer.removeAllFeatures();
  var feature = parser.read(geometryInput);
  vectorLayer.addFeatures(feature);
  if(dojo.isArray(feature)){
    feature = feature[0];
  }
  var bounds = feature.geometry.getBounds();
  bounds = bounds.scale(1.1);
  map.zoomToExtent(bounds);  
}

function addWkt(geometryText) {
  var wktFormat = new OpenLayers.Format.WKT();
  addGeometry(geometryText, wktFormat);
}

function addGeoJson(geometryText) {
  var geoJsonFormat = new OpenLayers.Format.GeoJSON();
  var geometry = JSON.parse(geometryText);
  geometry["type"] = "Feature";
  addGeometry(geometry, geoJsonFormat);
}

This are the most interesting parts of the code. The full source code is in the GeometryVisualizer.

Any questions/remarks/improvements ? Let me know !

Creating a node.js GeometryService : part 2

2011-12-15T21:38:00.016+01:00

This is my second post on what I learned from creating an implementation of the geometryservice specification in node.js. The first post can be found here.

By default node.js doesn't reload your files when they have changed but while developing this can be very handy. The most easy to use tool for node.js on windows that monitors your files for changes, with an easy way to set the node.js debugging flag, is nodemonw. You can download nodemonw at https://github.com/cenanozen/nodemonw. Once you've downloaded the executable I suggest to copy it to your %Appdata%\npm directory or another directory thats in your path. To start nodemonw I now run the following command nodemonw --debug index.js. In the next section it will become clear why I added the --debug flag.

To be able to debug my node.js application I installed node-inspector (npm install -g node-inspector). As you can read in the readme of node-inspector it is really easy to get started with node-inspector. You just have to start node-inspector and then open http://127.0.0.1:8080/debug?port=5858 in your favorite WebKit based browser. On Wikipedia I found this list of WebKit based browsers. The most well known ones for Windows are Google Chrome and Safari. A screencast on node-inspector can be found here. There is also a node-inspector playlist on YouTube.

To be able to test parts of my GeometryService I used vows (npm install -g vows). Vows is an asynchronous behavior driven development framework. More info about it can be found on http://vowsjs.org/. I am now going to show a small part of the code from my GeometryService and some tests I wrote for this code. My directory structure for the code I'll show looks like this: lib/ -- datatransformer.js test/ -- datatransformers.test.js In datatransformer.js I started a function to convert an ESRI geometry JSON object to wkt based on its geometry type.

exports.esriGeoJsonToWKT = function esriGeoJsonToWKT (geometryType, geometry) {
  if(geometryType === "esriGeometryPoint") {
    return "POINT(" + geometry.x + " " + geometry.y + ")";
  }
}

And the content of datatransformer.test.js is:

var dt = require("../lib/datatransformer");
var vows = require('vows');
var assert = require('assert');

vows.describe('Esri geometry JSON to WKT').addBatch({
  'when converting esriGeometryPoint {"x":-117,"y":34}':{
    topic: function(){ 
      var geomType = 'esriGeometryPoint';
      var p = JSON.parse('{"x":-117,"y":34}');
      return dt.esriGeoJsonToWKT(geomType, p);
    },
    'we get "POINT(-117 34)"': function(topic){
      assert.equal(topic, "POINT(-117 34)");
    }
  },
  'but when converting esriGeometryPoint {"x":-117.01,"y":34.02}':{
    topic: function(){ 
      var geomType = 'esriGeometryPoint';
      var p = JSON.parse('{"x":-117.01,"y":34.02}');
      return dt.esriGeoJsonToWKT(geomType, p);
    },
    'we get "POINT(-117.01 34.02)"': function(topic){
      assert.equal(topic, "POINT(-117.01 34.02)");
    }
  }
}).exportTo(module);

As you can see I created two tests for converting point geometries from ESRI JSON to wkt. One for integer coordinates and one for decimal coordinates. The easiest way to run this tests with vows on windows was opening a commandline in the root directory of my node.js project and type vows --spec. This will run the tests it finds in the test and spec directories of your project. My output looks like this :

  ? Esri geometry JSON to WKT



when converting esriGeometryPoint {"x":-117,"y":34}

V we get "POINT(-117 "4)"

but when converting esriGeometryPoint {"x":-117.01,"y":34.02}

V we get "POINT(-117.01 34.02)"



V OK » 2 honored (0.007s)

On a side note if you ever encounter that node can't find any of your globally installed packages then it might help to add a new User Variable called NODE_PATH with as value %AppData%\npm\node_modules. That was it for this post more posts on this project will follow.

Creating a node.js GeometryService : part 1

2011-12-15T21:00:00.000+01:00

Ever since I encountered the geoservices-rest-specification (pdf) I've been thinking about creating my own implementation. I also wanted wanted to try out node.js so I combined both ideas and started implementing a GeometryService in node.js. If you've never heard of node.js, this is the short introduction from the homepage of node.js:
Node.js is a platform built on Chrome's JavaScript runtime for easily building fast, scalable network applications. Node.js uses an event-driven, non-blocking I/O model that makes it lightweight and efficient, perfect for data-intensive real-time applications that run across distributed devices. The GeometryService is a web service that contains GIS related utility methods like project, intersect, buffer,... A sample server can be found here http://sampleserver1.arcgisonline.com/ArcGIS/rest/services/Geometry/GeometryServer.

I first downloaded and installed the windows installer of node.js version 0.65. After a reboot all locations where added to my path so that I can now start a new cmd window and type node. Note that when you install packages globally with the node package manager called npm (npm install <package> -g) npm (the node package manager) then the packages are installed in Windows 7 under %AppData%\npm (C:\Users\<username>\AppData\Roaming\npm). The full list of packages is located here search.npmjs.org.

To get started learning node I first read the node.js tutorial at nodebeginner.com. By following this tutorial I got a great insight in how to create your own webserver. At this moment I'm continuing to build my webservice based on this code but I plan to use some framework like express or journey in a later stage. This is not very urgent because I'll first focus on the JSON output of the GeometryService and thus won't need to output any html.

For my first version of the GeometryService I decided to use PostGIS as my geometry processor. I know this introduces an extra overhead. But I think this is the easiest way to get something up and running in a short period of time. If you don't want to bother with installing PostgreSQL and PostGIS I suggest you to download the Community edition of the Open Geo Suite. To be able to connect to the database I installed node-postgres with following command line: npm install pg -g. If you're on windows and get build errors you might try to add a file called true.cmd to the directory where node.exe resides (e.g. C:\Program Files (x86)\nodejs) with as content exit 0. Below is a short snippet on how I connect to the database. Make sure to replace all placeholders in the connection string.

var pg = require("pg");

function executeSQL(sql, parameters, resultCallback){
  var connectionString = "pg://username:password@host:port/databasename";
  pg.connect(connectionString, function(err, client) {
    client.query(sql,parameters, function(err, result) {
      // TODO add error handling
      console.log(result);
      resultCallback(result);
    });
  });
}

That was it for today. In the next part of the series I'll write about debugging node.js and about my progress on the GeometryService.

C# REPL

2011-12-03T06:07:00.004+01:00

Until recently when you wanted to test small snippets of C# code you had to create a small console application or use the immediate window while debugging or use something like LinqPad.

But a few weeks ago Microsoft finally announced a solution for this called Roslyn. The goal of Roslyn is to provide an API for the compiler. This is still just a CTP but it is very promising. One of the features of Roslyn is an interactive window for C# also called a REPL (read eval print loop). Note that some C# features like Linq query expressions, events and the dynamic and async keywords have not been implemented yet. More information on the Roslyn CTP can be found here Introducing the Roslyn CTP and here http://msdn.com/roslyn.

Tip of the day:
When you want to reevaluate or edit a previous entry then use ALT+Up and ALT+Down.

Bonus tip:
You might also wanna try out Jash. This is a javascript shell that can be opened on any website with a bookmarklet. It features some code completion and is very practical for trying out things you`re unsure about.

Performance of Checked vs Unchecked

2011-11-07T22:16:00.006+01:00

If you're ever in doubt on using the checked keyword in C# or the /checked compiler directive because it might hurt your performance then don't worry anymore.

I did some integer additions on my portable in a checked and unchecked context and the performance decrease was only about 3-5% on operations of about 5 seconds. Don't forget that checked and unchecked only have an effect on integers. Float and double never check for overflows and the decimal always does. More info on this feature can be found on msdn and in this Code Project article

This is the code I used to benchmark the performance of a checked vs unchecked context. I also tested subtraction, multiplication and division with similar results. Pay great attention when using the checked context. Only inline operations are checked/unchecked, code in function calls and anonymous functions isn't.

using System;
using System.Diagnostics;

namespace GISSolved.TestConsole
{
    class Program
    {
        static void Main(string[] args)
        {
            GC.Collect(3, GCCollectionMode.Forced);
            var checkedResult = TimeIt(true);
            Console.WriteLine("Checked : {0}", checkedResult);
            GC.Collect(3, GCCollectionMode.Forced);
            var uncheckedResult = TimeIt(false);
            Console.WriteLine("Unchecked : {0}", uncheckedResult);
            Console.WriteLine("Difference : {0:0}%", (1 - ((double)checkedResult / (double)uncheckedResult)) * 100);
            Console.ReadLine();
        }
        
        static long TimeIt(bool isCheckedCalculation)
        {
            var s = new Stopwatch();
            s.Start();

            Int32 a = 123, b = 123;
            if (isCheckedCalculation)
            {
                checked
                {
                    for (int i = 0; i < 900000000; i++)
                    {
                        a = i + b + a;
                        a = 125;
                    }
                    // extra check to make sure that overflow exceptions are thrown
                    try
                    {   
                        Int32 x = a + Int32.MaxValue;
                    }
                    catch (OverflowException)
                    {
                        Console.WriteLine("Int32 Overflow");
                    }
                }
            }
            else 
            {
                unchecked
                {
                    for (int i = 0; i < 900000000; i++)
                    {
                        a = i + b + a;
                        a = 125;
                    }
                    Int32 x = a + Int32.MaxValue;
                }
            }
            s.Stop();
            return s.ElapsedMilliseconds;
        }
    }
}

MongoDB : Geospatial Indexing

2010-06-16T12:51:00.003+02:00

Just a short note to let you know that this post has become outdated because MongoDB supports native two dimensional geospatial indexing. More information can be found here.

2 Possible Ways To Speed Up NetworkX (and Python)

2010-06-01T08:33:00.008+02:00

I'm in the process of looking for a way to run the shortest path algorithms of NetworkX in a faster way. At this moment I tried psyco and PyPy.

Psyco

Using psyco is really easy. After installing it you only have to write the following two lines.

import psyco psyco.full()

Psyco didn't speedup my shortest path calculations.

If you want to know more about psyco I suggest you to read the documentation.

PyPy

In order to test PyPy I downloaded the latest windows binary and ran my code against PyPy instead of CPython by invoking my script on the command line with pypy.exe.

The startup time is a lot longer for PyPy but if you don't take this startup cost into account then the speedup on my machine is in the range of 10-20% for shortest path routing with NetworkX.

So far my tests. Any suggestions ?

Using NetworkX with Jython

2010-05-31T21:56:00.004+02:00

Just a quick note to let you know that the Python graph library NetworkX (version 1.1), which by the way is really good, can be ported to Jython with some minor modifications. Sadly enough you loose some performance in the process.

The two changes I made to get it running where :

Because numpy doesn't exist for Jython you can't use the current_flow_closeness_centrality function in current_flow_closeness.py. I moved the import numpy as np statement at the top to the _compute_C function at the bottom.
In the __init__.py file under generators I commented from atlas import * out.

In the following weeks I will try to port a small subset off NetworkX to Java or Scala. Hopefully I'll be able to outperform the Jython and CPython version.

Arithmetic, Population, and Energy.

2010-05-27T21:36:00.001+02:00

Just watched Dr. Albert A. Bartlett's presentation on "Arithmetic, Population, and Energy." on youtube. It's really good and entertaining. It's in 8 parts and I encourage you to watch them all.

The url of part 1 is : http://www.youtube.com/watch?v=F-QA2rkpBSY

Using yield to create an iterator

2009-10-29T21:26:00.003+01:00

Not so long ago I blogged about using the foreach statement with an ICursor. I achieved this by inheriting from IEnumerator and IEnumerable.

But by using the yield statement we can achieve a similar effect with much less code. The generator function looks like this:

public static IEnumerable<IFeature> Iter(IFeatureCursor cursor)
{
    IFeature feat;
    while((feat = cursor.NextFeature()) != null)
    {
        yield return feat;
    }
    yield break;
}

You can use this method like below.

IFeatureClass featureClass; // initialize feature class
IQueryFilter queryFilter = new QueryFilterClass();
queryFilter.WhereClause = ... your where clause here ...
bool recycling = true;
IEnumerable<IFeature> features = Iter(featureClass.Search(queryFilter, recycling));
foreach (IFeature feature in features)
{
    ... your code here ...
}

If you're running .NET 3.0 or more you can create the following extension method in an extension class. This method adds the SearchIter method to the IFeatureClass interface which replaces the Search method.

public static class Extensions
{
    public static IEnumerable<IFeature> SearchIter(this IFeatureClass featureClass, IQueryFilter queryFilter, bool recycling)
    {
        IFeatureCursor cursor = featureClass.Search(queryFilter, recycling);
        IFeature feat;
        while ((feat = cursor.NextFeature()) != null)
        {
            yield return feat;
        }
        yield break;
    }
}

The usage of the extension method is similar to the usage of the Iter method but you can replace

IEnumerable<IFeature> features = Iter(featureClass.Search(queryFilter, recycling));

with this

IEnumerable<IFeature> features = featureClass.SearchIter(queryFilter, recycling);

So this was it. I hope you learned something and feel free to leave a comment.

Looping over a Workspace

2009-10-29T20:42:00.005+01:00

Ever wanted to execute a function on all the tables and/or featureclasses in a workspace then the following code is something for you. What it basically does is first loop over all the featureclasses that are in featuredatasets. Then loop over all the other featureclasses and finally loop over all the tables in the provided workspace.

To achieve this and still be able to retrieve the result of each function execution on the tables and featureclasses I created a generator function. I did this by using the yield statement. In order to avoid duplicate code I also created a nested function for looping over the featureclasses.

def executeiter(gp, workspace, featclassfunction=None, tablefunction=None):
    import os
    
    def loopfcs():
        fcs = gp.listfeatureclasses()
        for fc in iter(fcs.next, None):
            yield featclassfunction(fc)
    
    gp.workspace = workspace
    if featclassfunction is not None:
        for dataset in iter(gp.listdatasets().next, None):
            datasetworkspace = os.path.join(workspace, dataset)
            gp.workspace = datasetworkspace
            for result in loopfcs():
                yield result

        gp.workspace = workspace
        for result in loopfcs():
            yield result
            
    if tablefunction is not None:
        tables = gp.listtables()
        for table in iter(tables.next, None):
            yield tablefunction(table)

If you don't need the result of the function you can call the following function. It takes the same arguments as the executeiter function but its a regular function instead of a generator.

def execute(gp, workspace, featclassfunction=str, tablefunction=str):
    for x in executeiter(workspace, featclassfunction, tablefunction):
        pass

To demonstrate the usage of my code I first initialized a geoprocessing object and a workspace variable. Then I used the executeiter method to make it return the uppercase version of the name of the tables and featureclasses in my workspace. I could also have passed for example the describe or the listfields method of the geoprocessing object or a custom function. When you want to pass a function that doesn't return a result like the deleterows or deletefeatures function its more convenient to call the execute function.

import arcgisscripting, string

gp = arcgisscripting.create()
workspace = r'D:\temp\temp.gdb' # path to your workspace

# print the uppercase names of all the tables and featureclasses
for uppername in executeiter(gp, workspace, string.upper, string.upper):
    print uppername

# delete all rows of all the tables and featureclasses
execute(gp, workspace, gp.deletefeatures, gp.deleterows)

I've come to the end of this post. Did I miss something ? Know a Python idiom I really should start using ? Feel free to comment.

Editing with ArcObjects

2009-10-13T21:24:00.003+02:00

Recently I found myself rewriting a lot of code which started and stopped edit sessions to decorate it with try ... finally blocks to make sure that every started edit session closes even when an error occurs. The code I wrote looked like below.

bool saveEdits = false;
try
{
    StartEditing();
    ... edit a featureclass or table ...
    saveEdits = true;
}
finally
{
    StopEditing(saveEdits);
}

I thought that there had to be a cleaner way to achieve the same functionality. Here comes the using statement to the rescue. To be able to use the using statement with a class it has to implement the IDisposable interface. So I wrote the following wrapper class for starting and stopping edit sessions.

public class EditSession : IDisposable
{
    IWorkspaceEdit _workspaceEdit = null;

    public IWorkspaceEdit WorkspaceEdit
    {
        get { return _workspaceEdit; }
        set { _workspaceEdit = value; }
    }

    public EditSession(IWorkspace workspace, bool withUndoRedo)
    {
        if (workspace != null)
        {
            _workspaceEdit = (IWorkspaceEdit)workspace;
        }
    }

    public static EditSession Start(IWorkspace workspace, bool withUndoRedo)
    {
        EditSession editSession = new EditSession(workspace, withUndoRedo);
        editSession.Start(withUndoRedo);
        return editSession;
    }

    public void Start(bool withUndoRedo)
    {
        if (_workspaceEdit.IsBeingEdited() == false)
        {
            _workspaceEdit.StartEditing(withUndoRedo);
            _workspaceEdit.StartEditOperation();
        }
    }

    public void SaveAndStop()
    {
        Stop(true);
    }

    public void Stop(bool save)
    {
        if (_workspaceEdit.IsBeingEdited() == true)
        {
            _workspaceEdit.StopEditOperation();
            _workspaceEdit.StopEditing(save);
        }
    }

    #region IDisposable Members

    public void Dispose()
    {
        if (_workspaceEdit != null && _workspaceEdit.IsBeingEdited())
        {
            Stop(false);
        }
    }

    #endregion
}

You can use the EditSession class like this.

using (EditSession editSession = EditSession.Start(workspace, false))
{
    ... edit a featureclass or table ...
    editSession.SaveAndStop();
}

If you're using .NET 3.0 or a later version you can use the following class which creates an extension method for the IWorkspace.

public static class Extensions
{
    public static EditSession StartEditing(this IWorkspace workspace, bool withUndoRedo)
    {
        return EditSession.Start(workspace, withUndoRedo);
    }
}

With the extension method you can directly call StartEditing on a workspace.

using (EditSession editSession = workspace.StartEditing(false))
{
    ... edit a featureclass or table ...
    editSession.SaveAndStop();
}

Inserting Features and Rows

2009-09-27T11:59:00.006+02:00

Recently I was inserting features in a feature class but my code was rather slow so I looked around for another method and the following is what I found.

First you need to start an edit session and then you can use the code below. This is just a short body of code to give you an idea on how to use the different objects.

IFeatureCursor insertFeatCursor = outputFeatClass.Insert(true);

foreach (object featureToInsert in featuresToInsert)
{
    IFeatureBuffer outputFeatBuffer = outputFeatClass.CreateFeatureBuffer();
    // set the shape
    outputFeatBuffer.Shape = ...
    // set the different values
    outputFeatBuffer.set_Value(fieldIndex, value);
    // insert the feature buffer
    insertFeatCursor.InsertFeature(outputFeatBuffer);
}
insertFeatCursor.Flush();

If you now save and close your edit session the features are inserted. The code for inserting rows in a table is very similar.

ITable table;
ICursor insertCursor = table.Insert(true);
IRowBuffer rowBuffer = table.CreateRowBuffer();
rowBuffer.set_Value(fieldIndex, value);
insertCursor.InsertRow(rowBuffer);
insertCursor.Flush();

If you have any comments or questions, let me know !!!

Geoprocessing 1 : Intersect

2009-09-19T14:02:00.005+02:00

The intersect operation is one of the many overlay operations. When intersecting layer A with a layer B the result will include all those parts that occur in both A and B. In this post I'm going to compare the results of the ArcGIS Intersect tool with the corresponding operation in 3 other GIS packages. The reason why did this was that I noticed that the output from ArcGIS contains more features then I expected when dealing with overlapping geometries. This is usually not a problem but it becomes one when you have lots of overlapping geometries.

The GIS tools I used to compare the result of ArcGIS with were : uDig, Quantum GIS and gvSIG. All these and some more can be found on Portable GIS. Portable GIS brings open source GIS to your usb.

To test the code I prepared 2 shapefiles. One with a long small polygon and the other one with 3 overlapping rectangles (see image below).

When intersecting with ArcGIS (version 9.2 and 9.3 tested) the output looks like below. As you can see in the attribute table the result contains 9 polygons.

To be able to do the intersect operation with uDig I had to install the Axios Spatial Operations Extension. You can install this extension by clicking Help -> Find and install from the menubar. With Quantum GIS I needed to enable the ftools plugin to enable the geoprocessing functionality. The results of the 3 used open source GIS packages where the same. Only 3 features where created in the output shapefile. For completeness I added the screenshots of the results.

uDig :

Quantum GIS :

gvSIG :

Have any comments or questions ? Let me know !

Using foreach with the ICursor

2009-09-08T21:47:00.004+02:00

My first post on this blog was about how to loop over an ESRI cursor in Python with the for statement instead of the while statement. The same problem exists in .NET when using the various cursor objects like ICursor and IFeatureCursor. Normally you would code something like this :

ICursor cursor = table.Search(null, false);
IRow row;
while ((row = cursor.NextRow()) != null)
{
    // some code
}

Or this:

ICursor cursor = table.Search(null, false);
IRow row = cursor.NextRow();
while (row != null)
{
    // some code
    row = cursor.NextRow();
}

And I want to replace it with the following.

ITable table; // get table from somewhere
IQueryFilter queryFilter; // create QueryFilter or leave it null
CursorGS cursorGS = new CursorGS(table, queryFilter);
foreach (IRow row in cursorGS)
{
    // ... insert your code here
}

To make a class usable in a foreach statement you need to implement the IEnumerable and IEnumerator interfaces. I found a good introduction to making class usable in a foreach statement in this Microsoft article. So what I did was creating a class that inherited from IEnumerator<IRow> and IEnumerable<IRow> and implement all the needed properties and methods. For IEnumerator these where Current, MoveNext, Reset and Dispose and for IEnumerable only the method GetEnumerator was needed. The biggest problem was the Reset method because an ICursor doesn't have a reset method. I decided to set the cursor null and the recreated at the moment that MoveNext is called. This also makes sure the ICursor is only created when we really need it. Below you can find my full class

using System;
using System.Collections;
using System.Collections.Generic;
using System.Text;
using ESRI.ArcGIS.Geodatabase;

namespace GisSolved.GeoDiff.Esri
{
    public class CursorGS : IEnumerator<IRow>, IEnumerable<IRow>
    {
        ITable _table;
        IQueryFilter _queryFilter;

        ESRI.ArcGIS.ADF.ComReleaser _comReleaser;
        IRow _currentRow = null;
        ICursor _cursor = null;

        public CursorGS(ITable table, IQueryFilter queryFilter)
        {
            _table = table;
            _queryFilter = queryFilter;
            _comReleaser = new ESRI.ArcGIS.ADF.ComReleaser();
        }

        #region IEnumerator<IRow> Members

        public IRow Current
        {
            get { return _currentRow; }
        }

        #endregion

        #region IEnumerator Members

        object System.Collections.IEnumerator.Current
        {
            get { return _currentRow; }
        }

        public bool MoveNext()
        {
            if(_cursor == null) // initialize the cursor
            {
                _cursor = _table.Search(_queryFilter, false);
                _comReleaser.ManageLifetime(_cursor);
            }
            _currentRow = _cursor.NextRow();
            return _currentRow != null;
        }

        public void Reset()
        {
            _cursor = null;
        }

        #endregion

        #region IEnumerable<IRow> Members

        public IEnumerator<IRow> GetEnumerator()
        {
            return (IEnumerator<IRow>)this;
        }

        #endregion

        #region IEnumerable Members

        System.Collections.IEnumerator System.Collections.IEnumerable.GetEnumerator()
        {
            return (System.Collections.IEnumerator)this;
        }

        #endregion

        #region IDisposable Members

        public void Dispose()
        {
            _comReleaser.Dispose();
        }

        #endregion
    }
}

As you can see I used the ESRI ComReleaser object to make sure that the ICursor references get released properly.

So now you know how to create a wrapper around the ICursor you can start creating other wrappers for often used cursors like the IFeatureCursor, IFields, ... and I can go to my bed. Success and feel free to post any comments or your implementation !!!

DatagridView Tricks

2009-08-03T21:13:00.005+02:00

This post will entirely be about the .NET DatagridView control. I will show you how to synchronize the scrolling of two Datagridviews, disable column sorting, disable cell selection and focus cues.

To synchronize the scrolling of two DatagridViews you should subscribe to the Scroll event of the two Datagridviews and add the following code. What it does is equalizing the first row index and the scrolled horizontal offset. I did set the FirstDisplayedScrollingRowIndex because the VerticalScrollingOffset is a readonly property.

private void dataGridViewA_Scroll(object sender, ScrollEventArgs e)
{
    dataGridViewB.FirstDisplayedScrollingRowIndex = dataGridViewA.FirstDisplayedScrollingRowIndex;
    dataGridViewB.HorizontalScrollingOffset = dataGridViewA.HorizontalScrollingOffset;
}

private void dataGridViewB_Scroll(object sender, ScrollEventArgs e)
{
    dataGridViewA.FirstDisplayedScrollingRowIndex = dataGridViewB.FirstDisplayedScrollingRowIndex;
    dataGridViewA.HorizontalScrollingOffset = dataGridViewB.HorizontalScrollingOffset;
}

For the following tricks I created a custom control called DatagridViewGS that derives from the DatagridView control.

Disabling column sorting for all columns is done by setting the SortMode to NotSortable for every column of the DatagridView. In order to do that I subscribed the DataBindingComplete event and added the following code to it

void DatagridViewGS_DataBindingComplete(object sender, DataGridViewBindingCompleteEventArgs e)
{
    foreach (DataGridViewColumn c in Columns)
    {
        c.SortMode = DataGridViewColumnSortMode.NotSortable;
    }
}

Cell selection is disabled by subscribing to the CellStateChanged event and setting the Selected property of the cell whose state changed to false. I've also set the SelectionMode of the DataGridView to DataGridViewSelectionMode.CellSelect and the MultiSelect property to false.

void DatagridViewGS_CellStateChanged(object sender, DataGridViewCellStateChangedEventArgs e)
{
    if(e.StateChanged == DataGridViewElementStates.Selected)
        e.Cell.Selected = false;
}

After disabling cell selection there where still focus cues. These are the small dotted lines around the selected elements. To hide these I needed to override ShowFocusCues property of the DataGridView and make it always return false.

protected override bool ShowFocusCues
{
    get
    {
        return false;
    }
}

The full code of my custom DatagridView looked like this :

using System;
using System.Collections.Generic;
using System.Text;
using System.Windows.Forms;
using System.Drawing;

namespace GisSolved.GUI
{
    internal class DatagridViewGS:DataGridView
    {
        public DatagridViewGS()
        {
            MultiSelect = false;
            SelectionMode = DataGridViewSelectionMode.CellSelect;
            CellStateChanged += new DataGridViewCellStateChangedEventHandler(DatagridViewGS_CellStateChanged);
            
            DataBindingComplete += new DataGridViewBindingCompleteEventHandler(DatagridViewGS_DataBindingComplete);
        }
        
        // Disable column sorting
        void DatagridViewGS_DataBindingComplete(object sender, DataGridViewBindingCompleteEventArgs e)
        {
            foreach (DataGridViewColumn c in Columns)
            {
                c.SortMode = DataGridViewColumnSortMode.NotSortable;
            }
        }
        // Disable cell selection
        void DatagridViewGS_CellStateChanged(object sender, DataGridViewCellStateChangedEventArgs e)
        {
            if(e.StateChanged == DataGridViewElementStates.Selected)
                e.Cell.Selected = false;
        }
        // Disable focus cues
        protected override bool ShowFocusCues
        {
            get
            {
                return false;
            }
        }
    }
}

Do you know other useful tricks with DatagridViews or other .NET controls ? Feel free to write them down in the comments section !

Building and Installing Spatial Index on Ubuntu

2009-07-27T20:53:00.000+02:00

This is a short post about installing Spatial Index on Ubuntu. First download the latest release. At the moment it's version 3.2 and cd to the location of your download. Then issue the following commands :

tar xzvf spatialindex-1.3.2.tar.gz
cd spatialindex-1.3.2
./configure
sudo make install

If you want to configure your install you should take a look at the installation notes. You can find more information about Spatial Index at the projects trac.

Tokyo Cabinet 2 : Loading and querying point data

2009-07-16T21:48:00.006+02:00

After setting up Tokyo Cabinet and Ruby its time to use it. As with my post about MongoDB I'm going to load 500.000 POIs in a database and query them with a bounding box query. I will use the table database from Tokyo Cabinet because it supports the most querying facilities. With a table database you can query numbers with full matched and range queries and for strings you can do full matching, forward matching, regular expression matching,...

To load the data in my database I will need to read my shapefile with POIs with Ruby and write the attributes to a new database. First we create the database with the following code.

require 'tokyocabinet'
include TokyoCabinet

# create the object
tdb = TDB::new

# open or create  the database
if !tdb.open("poi_db.tct", TDB::OWRITER | TDB::OCREAT)
  STDERR.printf("open error: %s\n", tdb.errmsg(tdb.ecode))
end

To read the features in my shapefile I am going to use the Ruby bindings for GDAL/OGR. Because I installed Tokyo Cabinet on GISVM I already had FWTools installed but I still needed to install the Ruby bindings for it. I did this with the following command.

sudo apt-get install libgdal-ruby

Now we are going to read a shapefile with 500.000 point features and write the records to the database. First we open the shapefile and get the layer. Then we loop over the features, create a new record and fill the record with the x,y information and the other fields when they aren't empty. The values need to be converted to strings otherwise the record can't be saved. Then we put the record in the database.

require 'gdal/ogr'

# open my shapefile
dataset = Gdal::Ogr.open("poi_500000.shp")
layer = dataset.get_layer(0) 

feature_defn = layer.get_layer_defn

layer.get_feature_count.times do |i|
 record = Hash.new # create new record
    feature = layer.get_feature(i)
 geom = feature.get_geometry_ref()
 record['x'] = geom.get_x(0).to_s()
 record['y'] = geom.get_y(0).to_s()
 pkey = tdb.genuid # init primary key
 feature_defn.get_field_count.times do |i|
  field_defn = feature_defn.get_field_defn(i)
  fieldname = field_defn.get_name_ref
  value = feature.get_field_as_string(i);
  if not value.nil? and value != ""
   if field_defn.get_name_ref == "ID"
    pkey = value
   else
    record[fieldname] = value.to_s()
   end
  end
 end
 # store the record in Tokyo Cabinet
 tdb.put(pkey, record)
end

To add indexes on the x and y field we call the following code. This creates two supplementary files called poi_db.tct.idx.x.dec and poi_db.tct.idx.y.dec.

# add index on x and y
tdb.setindex('x', TDB::ITDECIMAL)
tdb.setindex('y', TDB::ITDECIMAL)

To query the POIs in the database I created a function to query the POIs for a given bounding box and then I benchmarked it. I used the same bounding box as in my previous posts about MongoDB, Rtree, Pythonnet and PostGIS.

# query POIs by bounding box
def query(tdb, minx, maxx, miny, maxy)
 qry = TDBQRY::new(tdb)
 qry.addcond("x", TDBQRY::QCNUMGE, minx.to_s())
 qry.addcond("x", TDBQRY::QCNUMLE, maxx.to_s())
 qry.addcond("y", TDBQRY::QCNUMGE, miny.to_s())
 qry.addcond("y", TDBQRY::QCNUMLE, maxy.to_s())
 qry.setorder("x", TDBQRY::QONUMASC)

 res = qry.search
 puts res.length # number of results found
 return res
end

require 'benchmark'
puts Benchmark.measure { query(tdb, 4.5, 5.0, 50.5, 51.0) }

The query returned 98000 POIs. I ran the benchmark 12 times and this where the results :

  1.620000   0.190000   1.810000 (  1.866339)
  1.570000   0.030000   1.600000 (  1.625303)
  1.640000   0.030000   1.670000 (  1.668573)
  1.650000   0.000000   1.650000 (  1.664806)
  1.650000   0.020000   1.670000 (  1.708228)
  1.730000   0.010000   1.740000 (  1.744645)
  1.410000   0.310000   1.720000 (  1.749268)
  1.620000   0.050000   1.670000 (  1.724199)
  1.610000   0.010000   1.620000 (  1.657794)
  1.660000   0.020000   1.680000 (  1.680383)
  1.710000   0.020000   1.730000 (  1.767141)
  1.720000   0.010000   1.730000 (  1.809114)

According to the Ruby documentation the benchmark outputs the user CPU time, the system CPU time, the sum of the user and system CPU times, and the elapsed real time. So this means that the query took between 1.65 and 1.87 seconds to get a list of 98000 POIs within the given bounding box. This is a nice indication of the speed of Tokyo Cabinet.

To demonstrate how you can access the attribute I created the following code. It loops over the first 100 found POIs and prints the ID and the x- and y-coordinate.

res = query(tdb, 4.5, 5.0, 50.5, 51.0)
# print the first hundred found POIs
i = 0
res.each do |rkey|
 rcols = tdb.get(rkey)
 puts rcols['id'].to_s() + " " + rcols['x'].to_s() + " " + rcols['y'].to_s()
 i += 1
 if i > 100
  break
 end
end

Now we are ready to close the database. I hope you enjoyed this post and as always I welcome any comments.

# close the database
if !tdb.close
 ecode = tdb.ecode
 STDERR.printf("close error: %s\n", tdb.errmsg(ecode))
end

Installing Tokyo Cabinet and Ruby on Ubuntu

2009-07-07T20:45:00.005+02:00

After MongoDB its time for another alternative to relational databases called Tokyo Cabinet. Tokyo Cabinet is a library of routines for managing a file based key-value store. It's a high performing database and it can be accessed over a network with Tokyo Tyrant. In this post I install Tokyo Cabinet, Ruby and the Ruby bindings for Tokyo Cabinet. But there will be a follow up post where I load and query POIs like I did with MongoDB and PostgreSQL/PostGIS.

Tokyo Cabinet only works on Linux so I installed it in on an Ubuntu virtual machine. It took me some time to figure everything out but if I can do it you can too. First you need to download the latest version of Tokyo Cabinet from the project site. Once downloaded you open a terminal window, navigate to the download location and issue the following commands.

tar xzvf tokyocabinet-1.4.20.tar.gz
cd tokyo-cabinet/
# install dependencies
sudo apt-get install checkinstall build-essential libbz2-dev
# now compile
./configure --prefix=/usr
make clean
make
# creates and installs a Debian package
sudo checkinstall -D

I decided to use the Ruby bindings so if you don't have it you can install it with the below command. This installs Ruby, an interactive shell, an interactive reference, the Ruby documentation and the dev part of the Ruby Standard Library. We will need the ruby-dev package for building the Ruby bindings of Tokyo Cabinet.

sudo apt-get install ruby irb ri rdoc ruby1.8-dev

After downloading the Ruby bindings for Tokyo Cabinet from the project site, I installed them with the following commands.

tar xzvf tokyocabinet-ruby-1.26.tar.gz
cd tokyocabinet-ruby-1.26
ruby extconf.rb
make
sudo make install

In the specifications document you can find a lot of information about Tokyo Cabinet and the underlying concepts. As you can read there are four types of databases : a hash database, a B+ tree database, a fixed-length database and a table database. In the examples directory of the Ruby bindings you will find some samples that create and use these four database. There is also a sample that uses the abstract database API with which you can communicate to the four database types. For more info on the Ruby bindings you can read the docs.

Problems with the installation ? Here are the sources that helped me with the installation process or post a comment and maybe I can help out.
http://openwferu.rubyforge.org/tokyo.html
http://oui.com.br/blog/nando-en/post/installing-tokyo-cabinet
http://www.ubuntugeek.com/how-to-install-ruby-on-rails-ror-in-ubuntu.html
http://blogs.law.harvard.edu/hoanga/2006/10/27/fixing-mkmf-load-error-ruby-in-ubuntu/

Drag Drop from ArcCatalog

2009-06-30T21:25:00.011+02:00

As you probably already noticed you can drag drop items within ArcGIS. For example you can drag a feature class from ArcCatalog to ArcMap or from ArcMap or ArcCatalog to a geoprocessing form. In this post I'm going to show you what has to be done to enable drag drop behavior from ArcCatalog and Windows Explorer to a textbox.

To enable drag drop on a textbox I added event handlers for DragEnter, DragOver and DragDrop for my textbox called textBoxPath. In the DragEnter and DragOver handlers I check whether the dragged object is valid. If this is the case the drag drop effect is set to All. This is a combination of the Copy, Move, and Scroll effect. But its in the DragDrop event handler and more precisely in the helper function GetPaths that the most import stuff happens. As you can see the Data property of the DragEventArgs is processed by the GetPaths method and if any paths to feature classes or tables are found the first path is shown. After that I added a small hack to put the cursor at the end of the text in the textbox.

private void TextBoxPath_DragEnter(object sender, DragEventArgs e)
{
 e.Effect = EsriDragDrop.IsValid(e.Data) ? DragDropEffects.All : DragDropEffects.None;
}

private void TextBoxPath_DragOver(object sender, DragEventArgs e)
{
 e.Effect = EsriDragDrop.IsValid(e.Data) ? DragDropEffects.All : DragDropEffects.None;
}

private void TextBoxPath_DragDrop(object sender, DragEventArgs e)
{
 List<string> paths = EsriDragDrop.GetPaths(e.Data);
 TextBox txtBoxPath = (TextBox)sender;
 if (paths.Count > 0)
 {
  // set value of textbox to the first found path
  txtBoxPath.Text = esriDragDrop.Paths[0];
  // place cursor at the end of the textbox
  txtBoxPath.SelectionStart = txtBoxPath.TextLength;
 }
}

In the below EsriDragDrop class I placed the IsValid and GetPaths methods. The IsValid method checks whether the IDataObject coming from the drag event contains any valid objects. The GetPaths method retrieves those valid objects and returns the paths to the found feature classes and tables. It uses the IDataObjectHelper interface and its GetNames and GetFiles methods to access the objects in the IDataObject. Note that only feature classes and tables will be returned by my code but this constraint can easily be removed by not checking the Type of the datasetName. I didn't add any functionality to check whether the file path dragged from a Windows Explorer to the textbox was valid but you can implement this by using the Geoprocessor object and its Exists and Describe methods or by trying to open the table or feature class.

using System.Collections.Generic;
using System.Windows.Forms;
using ESRI.ArcGIS.esriSystem;
using ESRI.ArcGIS.Geodatabase;
using ESRI.ArcGIS.SystemUI;

namespace GisSolved.DragDrop
{
 public class EsriDragDrop
 {
  const string DATAOBJECT_ESRINAMES = "ESRI Names";
  public static bool IsValid(IDataObject dataObject)
  {
   return dataObject.GetDataPresent(DATAOBJECT_ESRINAMES) ||
    dataObject.GetDataPresent(System.Windows.Forms.DataFormats.FileDrop);
  }
  public static List<string> GetPaths(IDataObject dataObject)
  {
   List<string> foundPaths = new List<string>();
   IDataObjectHelper dataObjectHelper = new DataObjectHelperClass();
   dataObjectHelper.InternalObject = (object)dataObject;
   
   if (dataObjectHelper.CanGetNames())
   {
    IEnumName enumNames = dataObjectHelper.GetNames();
    IName name;
    while ((name = enumNames.Next()) != null)
    {
     if (name is IDatasetName)
     {
      IDatasetName datasetName = (IDatasetName)name;
      // only accept feature classes and tables
      if (datasetName.Type == esriDatasetType.esriDTFeatureClass ||
       datasetName.Type == esriDatasetType.esriDTTable)
      {
       string path = System.IO.Path.Combine(datasetName.WorkspaceName.PathName, datasetName.Name);
       foundPaths.Add(path);
      }
     }
    }
   }
   else if (dataObjectHelper.CanGetFiles())
   {
    string[] paths = (string[])dataObjectHelper.GetFiles();
    foreach (string path in paths)
    {
     // TODO : Add code here to check if the file path is a valid path
     foundPaths.Add(path);
    }
   }
   return foundPaths;
  }
 }
}

This is all you need to implement drag drop behavior from ArcCatalog or Windows Explorer to your textbox. If you want to implement drag drop from ArcMap to your form I suggest you to read this and this. Any comments or suggestions ? Let me know.

Python Toolbox 4 : Clone Digger

2009-06-25T20:35:00.006+02:00

Looking for duplicate code or opportunities to refactor, let me introduce you to a great Python tool called Clone Digger. As the projects page says

Clone Digger aimed to detect similar code in Python and Java programs. The synonyms for the term "similar code" are "clone" and "duplicate code".

Once installed you call the clonedigger.py file with as arguments the path for the output html and the path to a folder or code file to analyze. If you call it with the parameter -h it outputs the different commandline options. To show the power of Clone Digger I used the following extract from actual code.

import arcgisscripting

gp = arcgisscripting.create()

def create_point(x, y):
    p = gp.createobject('Point')
    p.x = x
    p.y = y
    return p

def createPolygon(xMin, xMax, yMin, yMax):
    polygon = gp.createobject('array')
    ##Add the first point
    newPoint = createPoint(xMin, yMin)
    polygon.Add(newPoint)
    ##Add the second point
    newPoint = createPoint(xMin, yMax)
    polygon.Add(newPoint)
    ##Add the third point
    newPoint = createPoint(xMax, yMax)
    polygon.Add(newPoint)
    ##Add the fourth point
    newPoint = createPoint(xMax, yMin)
    polygon.Add(newPoint)
    ##Close the polygon
    newPoint = createPoint(xMin, yMin)
    #polygon.Add(newPoint)
    return polygon

To run Clone Digger on this file all you have to is issue the below command. Make sure that your shell finds the file clonedigger.py by adding to your path variables or by navigating to its folder.

python clonedigger.py -o D:\output.html D:\CodeTest.py

The output first shows some summary values from the code analysis. Then it shows the different code blocks where duplicate or similar code where found. This is how the output looks like for my short Python code.

Source files: 1

Clones detected: 2

9 of 17 lines are duplicates (52.94%)

Parameters
clustering_threshold = 10
distance_threshold = 5
size_threshold = 5
hashing_depth = 1
clusterize_using_hash = False
clusterize_using_dcup = False

Time elapsed
Construction of AST : 0.00 seconds
Building statement hash : 0.00 seconds
Building patterns : 0.00 seconds
Marking similar statements : 0.02 seconds
Finding similar sequences of statements : 0.00 seconds
Refining candidates : 0.02 seconds
Total time: 0.03
Started at: Wed Jun 24 21:05:15 2009
Finished at: Wed Jun 24 21:05:15 2009

Clone # 1
Distance between two fragments = 4
Clone size = 7

Source file "D:\CodeToTest.py" The first line is 16		Source file "D:\CodeToTest.py" The first line is 13
newPoint = createPoint(xMin, yMax)		newPoint = createPoint(xMin, yMin)
polygon.Add(newPoint)		polygon.Add(newPoint)
newPoint = createPoint(xMax, yMax)		newPoint = createPoint(xMin, yMax)
polygon.Add(newPoint)		polygon.Add(newPoint)
newPoint = createPoint(xMax, yMin)		newPoint = createPoint(xMax, yMax)
polygon.Add(newPoint)		polygon.Add(newPoint)
newPoint = createPoint(xMin, yMin)		newPoint = createPoint(xMax, yMin)

Clone # 2
Distance between two fragments = 4
Clone size = 5

Source file "D:\CodeToTest.py" The first line is 19		Source file "D:\CodeToTest.py" The first line is 13
newPoint = createPoint(xMax, yMax)		newPoint = createPoint(xMin, yMin)
polygon.Add(newPoint)		polygon.Add(newPoint)
newPoint = createPoint(xMax, yMin)		newPoint = createPoint(xMin, yMax)
polygon.Add(newPoint)		polygon.Add(newPoint)
newPoint = createPoint(xMin, yMin)		newPoint = createPoint(xMax, yMax)

Clone Digger is aimed to find software clones in Python and Java programs. It is provided under the GPL license and can be downloaded from the site http://clonedigger.sourceforge.net

I would not try remove all the duplicates or similarities found by Clone Digger. But I think it's a great tool to find code that can be improved. The degree to which you refactor depends greatly on the goal of your code and your time restrictions. Clone Digger is also be useful when working with multiple people on a project or when having to improve some legacy code.
Have any comments ? Any tools you can't live without ? Any suggestions for me ? Feel free to let me know.

PostGIS : Loading and querying data

2009-06-21T22:11:00.004+02:00

Today I'm going to load some POIs (points of interest) in a PostGIS database and query them back. To do this you'll first need an installed version of PostgreSQL with PostGIS. I suggest to first download and install PostgreSQL and then install PostGIS with the latest version. Then you should create a spatially enabled database. I called mine pois_db.

Now we are ready to import the POIs. I first extracted 500.000 POIs from a file geodatabase into a shapefile with the ArcGIS Select tool. Then I used the PostGIS shp2pgsql tool which is located in your bin folder of your PostgreSQL installation to create a text file. This file can then be used to import the data in your database with the PostgreSQL psql tool. My commandline looked like below.

rem -s : srid (spatial reference id)
rem -c : creates a new table and populates it
rem -D : use PostgreSQL dump format
rem -i create a GiST index on the geometry column
"C:\Program Files\PostgreSQL\8.3\bin\shp2pgsql.exe" -s 4326 -c -D -i poi_500000.shp pois_db > pois.sql

More information about the usage of shp2pgsql can be found here. Note that I set the srid to 4326 which stands for WGS84. On spatialreference.org you can find a list of spatial reference ids or you can take a look at the spatial_ref_sys table in your database.
Next thing we are going to is load the pois.sql in the pois_db with the following command line.

"C:\Program Files\PostgreSQL\8.3\bin\psql" -d pois_db -f pois.sql

To query the POIs we issue the following sql command. This query checks whether the bounding box of the point intersects with the bounding box of a created rectangle polygon.

SELECT *
FROM pois
WHERE the_geom && GeomFromText('POLYGON((4.5 50.5, 5.0 50.5, 5.0 51.0, 4.5 51.0, 4.5 50.5))',4326);

This query takes in pgAdmin III with a gist index on the geometry column only 515 ms for more then 4000 found POIs. If I only ask the ids the query is down to 67 ms which is slightly slower then with Rtree but faster then with MongoDb.

I've come to the end of this post. I hope you liked it and leave a comment if you want to.

Related Posts
Pythonnet
MongoDB
Rtree

Python Toolbox 3 : Pythonnet

2009-06-14T21:50:00.009+02:00

After logging and Pygments I want to uncover another useful Python tool called pythonnet. Pythonnet provides a way to call .NET assemblies from within Python. Its usage is very similar to that of IronPython and also comparable to that of ctypes which I used in my post about using the ArcGIS projection engine dll (pe.dll) to project coordinates.

To get working binaries for this project you need to make your hands a little bit dirty but with some small effort you're good to go. And although the last release is rather old I've read about people using it on Python 2.6 with .NET 3.5. Before anything I suggest you to read the readme because this will explain and clarify a lot of things. What I did is download the project from sourceforge and open the visual studio solution. Then I applied the patch described here. For more information about this patch and sample code to test whether the patch was applied successfully I suggest you to read this. Finally you have to build the solution and copiy the clr.pyd and Python.Runtime.dll to the root of your Python installation. In my case this is C:\Python24.
To test if everything went fine you can take a look at the demo directory in the pythonnet folder. I had to add import clr before the other import statements to make the samples work. This is not a bug but merely due to the fact that I used the standard Python executable.

One of the possible uses of pythonnet in a GIS context is when you're creating some arcgisscripting code but want to do things that can't be done with the Python bindings of ArcGIS. One of those things is a spatial filter although you can create a workaround with an extra query layer and the SelectLayerByLocation tool. I first created a static method called Hello to test if pythonnet worked with my dll. Then I created the CreateLayerFromSpatialIntersect method. The goal of this method is to do a bounding box query and export the queried features. This method takes input the path to a shapefile, a destination folder and name and the minimum and maximum values of x and y. The found features are exported to a new shapefile. I also added a helper method to get an IFeatureClass object from the shapefile. I also used the LicenseInitializer class. This class is part of the ArcGIS Console application template and has the sole purpose of getting and releasing an ArcGIS license.

using System;
using ESRI.ArcGIS.Geodatabase;
using ESRI.ArcGIS.Geometry;
using ESRI.ArcGIS.GeoDatabaseUI;
using ESRI.ArcGIS.esriSystem;
using ESRI.ArcGIS.DataSourcesFile;
namespace GisSolved.SpatialQueryTool
{
    public class SpatialQueryTool
    {

        public static string Hello()
        {
            return "Hello !!!";
        }
        public static void CreateLayerFromSpatialIntersect(string shapeSourcePath, string destinationFolder, string destinationName, double xMin, double xMax, double yMin,double yMax)
        {
            LicenseInitializer licenseInitializer = new LicenseInitializer();
            licenseInitializer.InitializeApplication(
    new esriLicenseProductCode[] { esriLicenseProductCode.esriLicenseProductCodeArcView },
    new esriLicenseExtensionCode[] { });

            IFeatureClass sourceFeatureClass = GetShapeFile(shapeSourcePath);

            // create  the spatial filter
            ISpatialFilter spatialFilter = new SpatialFilterClass();
            spatialFilter.SpatialRel = esriSpatialRelEnum.esriSpatialRelIntersects;

            IEnvelope boundingBox = new EnvelopeClass();
            boundingBox.XMax = xMax;
            boundingBox.XMin = xMin;
            boundingBox.YMax = yMax;
            boundingBox.YMin = yMin;

            spatialFilter.Geometry = (IGeometry)boundingBox;

            // create other parameters for the export operation
            IDataset sourceDataset = (IDataset)sourceFeatureClass;
            IDatasetName sourceDatasetName = (IDatasetName)sourceDataset.FullName;

            IFeatureClassName destinationClassName = new FeatureClassNameClass();
            IDatasetName destinationDatasetName = (IDatasetName)destinationClassName;

            destinationDatasetName.Name = destinationName;

            IWorkspaceName destinationWorkspaceName = new WorkspaceNameClass();
            destinationWorkspaceName.PathName = destinationFolder;
            destinationWorkspaceName.WorkspaceFactoryProgID = "esriDataSourcesFile.ShapefileWorkspaceFactory";
            destinationDatasetName.WorkspaceName = destinationWorkspaceName;

            // do the actual export
            IExportOperation exportOp = new ExportOperationClass();
            exportOp.ExportFeatureClass(sourceDatasetName, (IQueryFilter)spatialFilter, null, null, destinationClassName, 0);

            licenseInitializer.ShutdownApplication();
        }

        private static IFeatureClass GetShapeFile(string path)
        {
            string workspacePath = System.IO.Path.GetDirectoryName(path);
            string name = System.IO.Path.GetFileName(path);
            IWorkspaceFactory shapefileWorkspaceFactory = new ShapefileWorkspaceFactoryClass();
            IFeatureWorkspace workspace = (IFeatureWorkspace)shapefileWorkspaceFactory.OpenFromFile(workspacePath, 0);
            return workspace.OpenFeatureClass(name);
        }
    }
}

It's easy to use this code from pythhonnet. The first step is to build the .NET project and copy the resulting dll to the same directory as your Python code. Then you need to import clr and add a reference to the dll. Next thing to do is import the desired class from your assembly. Now you can use your class like you would do in IronPython. To test my code I exported 500.000 POIs to a shapefile which I then queried with the same bounding box as in my two previous posts about MongoDb and Rtree.

import clr
clr.AddReference('SpatialQueryTool')
from GisSolved.SpatialQueryTool import SpatialQueryTool

print SpatialQueryTool.Hello()

source= r"D:\source\poi_500000.shp";
destinationFolder = r"D:\destination";
destinationName = "poi_45_50_505_510_2";
SpatialQueryTool.CreateLayerFromSpatialIntersect(source, destinationFolder, destinationName, 4.5, 5.0, 50.5, 51.0)

I hope that this post was useful to you and as always I welcome any comments.

Spatial indexing a MongoDb with Rtree

2009-06-06T18:29:00.009+02:00

Update : recent versions of MongoDB support native geospatial indexing.

In my previous post about populating a MongoDb I mentioned the existence of the Rtree package as a solution to the lack of a spatial index mechanism in MongoDb. In this post 'm going to try this out. Rtree is a Python package with the sole purpose of creating spatial indexes. To install Rtree I downloaded the Windows installer and ran the installer.

To test Rtree I decided to create a spatial index for my previously created POI database in MongoDb. To do this I first connect to my running MongoDb instance. Then I create an Rtree disk index in the same directory as my script called poi_index. This creates two files (poi_index.dat and poi_index.idx) who will contain the index. Finally I loop over all the POIs in the database and add the x and y values to the index with as reference the id of the POI.

from pymongo.connection import Connection
from rtree import Rtree

mongo_conn = Connection('localhost', 27017)
poidb = mongo_conn.poidb
pois = poidb.pois

poi_index = Rtree("poi_index", overwrite=1)

for poi in pois.find():
    poi_index.add(poi["ID"], (poi["x"], poi["y"]))

Now we have created an index for the 500.002 POIs in my database. The file size of this index is 350 mb which is more than I expected it to be when compared my Mongo POI database which has a total file size of 460 mb. Querying this index is very easy and very fast. To time my query I added a handy wrapper function. Essentially to query the index the only thing you need to do is create a bounding box and call the intersection method of the index. This returns the ids of the intersecting points.

import time

def print_timing(func):
    def wrapper(*arg):
        t1 = time.time()
        res = func(*arg)
        t2 = time.time()
        print '%s took %0.0f ms' % (func.func_name, (t2-t1)*1000.0)
        return res
    return wrapper

minx, maxx = 4.5, 5.0
miny, maxy = 50.5, 51.0
bbox = ( minx, miny, maxx, maxy )

@print_timing
def intersect(index, bbox):
    print len(index.intersection(bbox)), "hits"
intersect(poi_index, bbox)

On my machine this bounding box query took 265 ms on the MongoDb and only 15 ms with Rtree. This is a great improvement for my query. But the disadvantage is that I still need to query my POIs out of the MongoDb with the ids I got from Rtree. Sadly enough this doesn't go very fast so the total improvement is very small and with lots of points even negative. Strange thing is that when I do an IN query on MongoDb the result is always slower than querying the pois one by one by ID. My temporary conclusion is that Rtree is really fast but you only get the ids. The best solution would be that MongoDb could natively use the Rtree index or another spatial index. A small downside of Rtree is that the id can only be an integer value. So you can't store strings as ids. We have come to the end of my small exploration of Rtree as a mean to add a spatial index to MongoDb. I hope you liked it and I still welcome any comments/suggestions.

Populating a MongoDb with POIs

2009-05-31T11:37:00.014+02:00

I have been looking at the alternatives for relational databases for a while now and this week I came across MongoDB. MongoDB is a documented-oriented database which tries to bridge the gap between key/value stores and relational databases. In this post I'm going to install MongoDb and populate it with 500.000 POIs from Tele Atlas.

To install MongoDb I only needed to download the Windows binaries and that's it. Then I followed the guidelines in the Quickstart to start the database. As the default datapath for MongoDb on Windows is c:\data\db I specified a custom datapath.

> mongod.exe --dbpath "D:\mongodb\data" run

Once I got MongoDb up and running I downloaded and installed the Python bindings named PyMongo. This also went very smoothly. On the MongoDb website there is a great introductory tutorial for PyMongo.

To populate the MongoDb I wrote a Python program that loops over the POIs in a layer, creates a dictionary with its properties and coordinates. Then the POI dictionary is inserted in the pois collection in the poi database.
First I connect to my running MongoDb and create a reference to a new database called poidb. Then I create a new collection called pois.

import arcgisscripting
from pymongo.connection import Connection

mongo_conn = Connection('localhost', 27017)
poidb = mongo_conn.poidb
pois = poidb.pois

Then comes some code that uses ESRI functions to create a list with the fields in the POI layer and a featurecursor for the POIs. I removed the SHAPE field from this list because this will be treated separately.

gp = arcgisscripting.create()
poi_path = r'D:\ta.gdb\eureur_____pi'

fields = gp.listfields(poi_path)
names = [f.name.upper() for f in iter(fields.next, None)]
names.remove('SHAPE')
fc = gp.searchcursor(poi_path)

In the next part we loop over the POIs in the cursor. For each poi we extract the x and y coordinate and save them in the poi dictionary. For all the fields we fetch the value and put the none empty values in the dictionary. The strings are stored as unicode. Finally we insert the poi dictionary in the poi collection of the poidb. The script stops when more then 500.000 POIs have been inserted (500.002 to be exact).

for index, feat in enumerate(iter(fc.next, None)):
    poi = {}
    point = feat.shape.getpart(0)
    poi['x'] = point.x
    poi['y'] = point.y
    for field in names:            
        value = feat.getvalue(field)
        if value and len(str(value).strip()) > 0:
            if isinstance(value, str):
                value = unicode(value, 'latin-1')
            poi[field] = value

    pois.insert(poi)
    if index > 500000:
        break

To count the pois you only need to call the count method of the pois collection object.

print pois.count()

Before querying this data I'am going to create ascending indexes for the x and y fields. It would be better if spatial indexes where supported but with reasonably sized datasets a double index is sufficient for doing bounding box queries. I also could try the RTree package but that will be for another time. At the end the information about the indexes is printed.

from pymongo import ASCENDING
poi_db.pois.ensure_index('x', ASCENDING)
poi_db.pois.ensure_index('y', ASCENDING)

print pois.index_information()

Query parameters are passed to MongoDb in a nested directory structure. So defining a bounding box query is an straightforward task. The $gt character means bigger then and $lt means smaller then. In this query I search for all the POIs with a longitude larger then 4.5 and smaller then 5.0 and a latitude larger then 50.5 and smaller then 51.0. Then I print the result count, the number of milliseconds that the query took and the first resulting document.

query_pois = pois.find({"x":{"$gt":4.5, "$lt":5.0}, "y":{"$gt":50.5, "$lt":51.0}})

print query_pois.count()
print query_pois.explain()["millis"]
print query_pois.next()

The result looks like this :

4585

265

SON([(u'HSNUM', u'173'), (u'MUNCD', u'23062'), (u'MUNID', 10560032003493.0), (u'x', 4.5000037996786091), (u'NAME', u'Shell Overijse Frans Verbeekstraat'), (u'OBJECTID', 134471), (u'TELNUM', u'+(32)-(2)-6573607'), (u'STNAME', u'Frans Verbeekstraat'), (u'FEATTYP', 7311), (u'ADDRPID', 10560000020526.0), (u'RELPOS', 75), (u'POSTCODE', u'3090'), (u'STNAMELC', u'DUT'), (u'ARNAMELC', u'DUT'), (u'y', 50.769815399855986), (u'IMPORT', 2), (u'MUNNAME', u'Overijse'), (u'_id', ObjectId('YB\xd6l\x99\xe0$U\x8f\t>\xc3')), (u'CLTRPELID', 10560001815786.0), (u'ID', 10560300013584.0), (u'BUANAME', u'Hoeilaart')])

So it took 265 milliseconds to query 4585 POIs from a collection 500.002 POIs with two range query parameters. If you need to optimize a query you can take a look at MongoDb Site for some tips and tricks.

We have come to the end of my first exploration of MongoDb. What I especially liked is the flexibility you get from this kind of databases and the ease of installation and use. The downside for geographic applications is that at the moment there is no built-in support for geometries. I hope you enjoyed this post and I welcome any comments and suggestions.

Python Toolbox 2 : Pygments

2009-05-19T22:29:00.004+02:00

This is the second post in my series about Python tools and I'm going to talk about Pygments. Pygments is a generic syntax highlighter implemented in Python. I started using it to able to show the code samples on this blog in a decent way.

Pygments can highlight a very large list of programming languages, template languages and other sources like config files, shell scripts, css and many more. The markup of your sources can then be exported to multiple output formats like HTML, GIF or BMP image, Latex, rtf, svg, Bulletin bord code and terminal. More info about the formatters can be found here.

Once installed you can call this library from within Python or you can call the provided executable from the commandline. The documentation for Pygmentize is very clear and extensive so I won't reproduce it here but I'll show you some small samples.

Below is a sample commandline usage. The first command outputs a file with your code wrapped in a div and span elements with css class attributes. The second command creates the stylesheet for your code.

pygmentize -o test.html test.py
pygmentize -S default -f html > style.css

If you want to create a full html file from your code with the styles included then all you need is the following command :

pygmentize -O full=true -o test.html test.py

For the code in my blog I use the following simple Python script that converts the scripts in a folder to a html files. I also added the styles to the css of my blog.
First I import the relevant namespaces from pygments and the os namespace. And then create a basic HtmlFormatter

from pygments import highlight
from pygments.lexers import PythonLexer
from pygments.formatters import HtmlFormatter
import os

formatter = HtmlFormatter()

If you want a printout of all the different classes the only thing you need is the following line. I ran this once and copied to my style definitions on my blog.

print formatter.get_style_defs('.highlight')

The following code converts all the python files in a directory to html files with the same name. We loop over the elements in the folder, read the python files, create a highlight of them and then write this to a new file.

indir= r''
for x in os.listdir(indir):
    if x.endswith('.py'):
        infile = os.path.join(indir, x)
        outfile = infile.replace('.py', '.html')
        code = ''.join(list(open(infile)))
        html = highlight(code, PythonLexer(), formatter)
        f = open(outfile, 'w')
        f.write(html)
        f.close()
        print 'finished :',x

Projections and transformation (2)

2009-05-12T21:12:00.012+02:00

In a comment on my previous post about projections and transformations Gene pointed out that I could use Ogr with it its Python bindings. I already used the C# bindings from Ogr in the past but I thought it would be a nice exercise to use Ogr from Python.

If you are on windows and need to install gdal/ogr and its Python bindings I suggest you to use the following installation guideline.

Below is the short code I came up with but there are numerous other samples on the web sample 1 sample 2 and off course the documentation on the gdal site. Attention I got an error when trying to use the TransformPoints method.

from osgeo import osr

def get_spatialref(epsg_code):
    spatialref = osr.SpatialReference()
    spatialref.ImportFromEPSG(epsg_code)
    return spatialref

lambert72_spatialref = get_spatialref(31300)
etrs89_spatialref = get_spatialref(4258)

coordinate_transformation = osr.CoordinateTransformation(lambert72_spatialref, etrs89_spatialref)

## start location : random point in Belgium
x = 155335.20410
y = 166096.08270
print coordinate_transformation.TransformPoint(x, y)

Another solution is to install and use the Python bindings for the Proj4 library.
The following sample code comes from its documentation.

import pyproj

# projection 1: UTM zone 15, grs80 ellipse, NAD83 datum
# (defined by epsg code 26915)
p1 = pyproj.Proj(init='epsg:26915')
p1_2 = get_spatialref(26915)
# projection 2: UTM zone 15, clrk66 ellipse, NAD27 datum
p2 = pyproj.Proj(init='epsg:26715')
p2_2 = get_spatialref(26715)
# find x,y of Jefferson City, MO.
x1, y1 = p1(-92.199881,38.56694)
# transform this point to projection 2 coordinates.
x2, y2 = pyproj.transform(p1,p2,x1,y1)
print x2, y2

If you need to look up a spatial reference you might wanna take a look at spatialreference.org.

Succes !

Python Toolbox 1 : Logging

2009-05-04T10:07:00.032+02:00

This is the first post in a series about internal and external Python modules. Today I am going to give a brief overview of the built in logging module.

The logging module is a core python module that was introduced in Python 2.3. There are 6 different logging levels (debug, info, warning, error, critical and fatal). You can also log an exception wich outputs a nice traceback. There are numerous output possibilities defined as handlers. The log information can for example be send to a file, an http server, the windows event log, an e-mail address, a stream and more. The logger can also be configured to log to a different output depending on the logging level. You could for example log to a file when the level is debug but you might want to receive an e-mail when a critical or fatal error occurs. In this post we will use a basic file logger and then create a rotating file logger. The rotating file logger will make a back-up of the logfile when it reaches a predefined maximum size.

As the following code shows, basic logging is easy to set up. In short it first imports the logging module. Then it intialize this module by calling the basicConfig method. According to the help the input parameters for this method are :

filename  Specifies that a FileHandler be created, using the specified filename,
          rather than a StreamHandler.
filemode  Specifies the mode to open the file, if filename is specified
          (if filemode is unspecified, it defaults to "a").
format    Use the specified format string for the handler.
datefmt   Use the specified date/time format.
level     Set the root logger level to the specified level.
stream    Use the specified stream to initialize the StreamHandler. Note that
          this argument is incompatible with 'filename' - if both are present,
          'stream' is ignored.

Then we test the different logging level and also deliberately throw an error to test the exception logging. At the end I added some clean up code because the logging module doesn't get unloaded between to test runs in some IDE's.

import logging

logging.basicConfig(filename='gissolved.log',
                    level=logging.DEBUG,
                    format='%(asctime)s %(levelname)s : %(message)s')

try:
    logging.debug('debug')
    logging.info('info')
    logging.warning('warning')
    logging.error('error')
    logging.critical('critical')
    logging.fatal('fatal')
    try:
        raise TypeError, 'Demo error'
    except TypeError:
        logging.exception('exception')
finally: ## close logger properly
    logging.shutdown()
    logging.getLogger(None).handlers = []

Ok now let's create a rotating file logger. To do this we not only import the logging module but also the sub module handlers and the class Formatter. Then we define the log path and create a logger object. Now its time to create the RotatingFileHandler that is defined in the handlers sub module. The parameters are the filename, the mode (defaults to 'a' which means append), maxBytes (max number of bytes before rollover occurs), backupCount (number of backup files that will be created as needed). Then we define a format for the handler and add it to the logger object. At last we set the logging level of our logger.

Now its time to test the logger. I've put the code in a loop in order to log more bytes then the predefined number of maximum bytes per logfile. Internally the rotating file logger first checks whether the max size will be reached by writing the current message to the logfile. If this is the case it closes the current file and adds a suffix to it, a new file with the message is then created. So in this case the first back-up file will be renamed to gissolved_rotate.log.1. The next back-up will then be named gissolved_rotate.log.2 and so on. Until the maximum number of back-up files is reached and the first back-up files gets overwritten. The file being written to is always the original gissolved_rotate.log. In the finally block the current handler is removed from our logger and the handler is closed. If you've defined multiple handlers for your logger you could loop over the logger.handlers and call their close method and then set the logger.handlers equal to an empty list like I did in the first example.

import logging
from logging import handlers, Formatter

try:
    log_path = 'gissolved_rotate.log'
    logger = logging.getLogger('gissolved')
    ## create a rotating file handler and add it to the logger
    handler = handlers.RotatingFileHandler(log_path,
                                           maxBytes=5000000,
                                           backupCount=4)
    log_format = Formatter("%(name)s %(asctime)s %(levelname)s : %(message)s")
    handler.setFormatter(log_format)
    logger.addHandler(handler)
    ## set the root logging level
    logger.setLevel(logging.DEBUG)
    
    for i in range(11000):
        ## write test logging
        logger.debug('debug')
        logger.info('info')
        logger.warning('warning')
        logger.critical('critical')
        logger.error('error')
        logger.fatal('fatal')
        
        try:
            raise TypeError, 'Demo error'
        except TypeError:
            logger.exception('exception')

        logger.info('next round ' + str(i))

finally: ## close logger properly
    logger.removeHandler(handler)
    handler.close()

If you've still have some questions, shoot and I'll try to answer them.

Table to csv

2009-04-30T07:38:00.006+02:00

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

2009-04-20T21:58:00.035+02:00

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

2009-04-16T20:36:00.020+02:00

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