GIS Solved

Creating a node.js GeometryService : part 3

2012-04-04T21:06:00.001+02:00

Its been a while since I posted my first and second post on a GeometryService based on node.js and PostGIS.

Today I want to announce that I open sourced my work in progress. It is in a pre-alpha state and it can be found here. The only thing it can do at the moment is convert ESRI Json geometries from one coordinate system to another. This is done by converting the geometry to WKT, creating a sql statement for the re-projection and then converting the GeoJSON result back to ESRI Json.

A sample input url when you started the geometryservice is : http://127.0.0.1:3000/rest/services/Geometry/GeometryServer/project?inSR=4326&outSR=3819&geometries={%22geometryType%22:%22esriGeometryPoint%22,%22geometries%22:[{%22x%22:-117,%22y%22:34},%20{%22x%22:-115,%22y%22:25}]}. This returns

{"geometries":[{"x":-115.00495449693217,"y":24.99533030379053},{"x":-117.00524224953409,"y":33.99507554185914}]}

as a result. The intermediate sql statement looks like this :

SELECT st_asgeojson(st_transform(st_geomfromtext(geometries.geometry, '4326'), '3819')) g FROM  (SELECT ST_ASTEXT('POINT(-117 34)') geometry UNION SELECT ST_ASTEXT('POINT(-115 25)')) geometries

Most of the logic of the code happens in geometryService.js. The conversion from and to ESRI Json happens in datatransformer.js. As stated in my previous post I also created some tests for the functions in the datatransformer.js. There is still a lot to do, especially everything related to error handling and logging and off course the implementation of the other geometry operations.

The main dependency is a runnning instance of PostgreSQL with PostGIS. As the connection string is not configurable yet you'll have to change it in the executeSQL function in geometryService.js.

The code can be found on following url https://bitbucket.org/gissolved/node_geometryservice and it is licensed under the MIT license.

If you like this then consider buying me a book or a license for Sublime Text 2.

Timing functions in Python

2012-02-12T13:54:00.001+01:00

Recently I saw a blogpost on tracking time in Python. This is inspired me to share how I track time in my python scripts. First I put the following python decorator function in my code or preferably import it from a common module.

import time                                                

def timeit(method):

    def timed(*args, **kw):
        ts = time.time()
        result = method(*args, **kw)
        te = time.time()

        print('%r (%r, %r) %2.2f sec' % (method.__name__, args, kw, te-ts))
        return result

    return timed

You can use the timeit function like this:

@timeit
def longRunningOperation(dummyArg, optionalDummyArg=None):
    time.sleep(3)
    
longRunningOperation("test", "optional test")

The output of the timeit decorator function looks like this:

'longRunningOperation' (('test', 'optional test'), {}) 3.00 sec

As you can see, the timit function reports the function name, the passed arguments and the elapsed time

If you want to learn more about decorator functions I would suggest to read the python.org wiki article on decorators. On the same site there is also an article with useful decorator functions for things like caching, properties, pre- and post conditions, memoization and many more.

Roundup January 2012

2012-02-05T20:01:00.000+01:00

This post is a roundup post for the month january.

Blogging

Just before the end of the year I created a tumblr blog where I shared some inspiring or interesting quotes.

On this blog I wrote 2 posts. The first post of january was my second post on developing a GeometryService in node.js. Some things I talked about where nodemonw for monitoring your files for changes and restarting node.js, node-inspector for debugging and vowsjs for writing tests. The second post of this month elaborated on the Geometry Visualizer that I created this month. The Geometry Visualizer uses the ArcGIS JavaScript API and openlayers to provide a .

GIS/Geography

List of free GIS datasets
GDAL/OGR 1.9 with support for lots of new drivers like DWG, CouchDB, ESRI FileGDB, SVG, Google Fusion Tables and many more new goodies.
Designing a Minimalist OpenStreetMap Baselayer for MapBox
Tilemill now also available for windows

Programming

Clean Project: a simple tool for e-mailing Visual Studio solutions.
SqlDiffFramework: a visual differencing engine for dissimilar data source.
FluentMigrator: a migration framework for .net much like Ruby Migrations.
Use DateTime.UtcNow, it's faster and you never have problems with daylight saving time 1 hour jumps (Scott Hanselmann got bitten by this one).
11 Things every Software Developer should be doing in 2012 : Twitter, StackOverflow, blog, get out there, have a modern phone, embrace mobile, learn a design pattern, set a reachable goal, learn a programming language, boost your confidence and read programming blogs,/books/magazines.
Haskell for the Evil Genius

Teaching/learning

Physicists Seek To Lose The Lecture As Teaching Tool
Learn jQuery in 30 days
Learn JavaScript programming this year thanks to Codecademy (Ruby and Python are in the making).
Learn jQuery in 30 days with tuts+.
Follow a Stanford course online, see at the bottom for more courses.
I'm reading Real World Functional Programming and here is why you should pick up a functional programming language.

Other

Facebook, Google, ... Tax avoidance 2.0
Minimum viable SEO (long tail vs head tail, use bread crumb navigation and use page titles.
Duck programming, the pitfalls of building a system around rules engines, programmable workflow or domain-specific language and how to avoid them.
Your time is $1000/hour, and you need to act accordingly: How should a startup founder value her time?
Small teams are dramatically more efficient than large teams
The Case Against Buying Christmas Presents by Leo Babauta.
Awesome answers on the question "How accurate are the power numbers from a Tacx Flow?" at bicycles.stackexchange.com. Yet another reason to buy a power meter. Although not everyone agrees.
Salary Negotiation: Make More Money, Be More Valued.

Entertainment

Youtube: "Video Game Math is The New Reading" and "Shopping Math is The New Free
Funky music: Goldfish - Get Busy Living
Amazing C book review on Amazon
Fascinating: Smart crows and ravens

New blogs to follow

The Morning Brew: daily .NET software development link blog.
ViperChill: a blog about viral marketing, the best post I read was The Post on Productivity I Wish Someone Else Had Written

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 options = {
  maxExtent: new OpenLayers.Bounds(-1000000000, -1000000000, 1000000000, 1000000000),
}
var map = new OpenLayers.Map('openLayersMap', options);

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.

As I announced in this post, I open sourced the code of this project. It can be found in this bitbucket repository.

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.

As I announced in this post, I open sourced the code of this project. It and can be found in this bitbucket repository.

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