Samuel Bosch

Showing posts with label ArcGIS. Show all posts

Using yield to create an iterator

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

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

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

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

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

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

Drag Drop from ArcCatalog

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 3 : Pythonnet

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.

Populating a MongoDb with POIs

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.

Table to csv

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

import arcgisscripting, string
from cursoriterator import CursorIterator

gp = arcgisscripting.create()

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

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

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

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


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

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

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

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

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

Finally we write the outputlist to the defined outputpath.

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

The full code can be found on the arcscripts site.

Projections and transformation

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

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

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

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

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

    import ctypes

The pe.dll is accessed in the following way.

    pe = ctypes.windll.pe

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

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

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

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

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

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

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

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

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

import ctypes

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

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

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

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

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

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

I hope you liked this post. Have fun.

Pages