7

I've seen lots of discussion on this, especially using IIndexQuery.NearestFeature() on a FeatureIndex instance. Something like this:

private IdDistancePair FindNearest(IGeometry g, IFeatureClass fc, string where = "")
{
IFeatureIndex fi = new FeatureIndexClass();
fi.FeatureClass = fc;

//you can optionally set a FeatureCursor on the IFeatureIndex 
//to reduce the number of features indexed.
//also, limiting subfields may speed up the indexing operation
IQueryFilter qfilter = null;
if (!string.IsNullOrEmpty(where))
{
  qfilter = new QueryFilterClass();
  qfilter.WhereClause = where;
  IFeatureCursor fcursor = fc.Search(qfilter, true);
  fi.FeatureCursor = fcursor;
}

//create the index
fi.Index(null, ((IGeoDataset)fc).Extent);

int fid;
double distance;
((IIndexQuery)fi).NearestFeature(g, out fid, out distance);

return new IdDistancePair(fid, distance);
}

However, it appears that the FeatureIndex is an in-memory thing unrelated to any existing spatial index on the query feature class. It takes a long time (dozens of seconds) to build. (and the above doesn't work for me anyway, always returning a FID of -1 after a 45s wait).

I only need to do one query per layer, repeated for ~5 layers. The docs on the FeatureIndex class suggest that the above approach is expensive in such a circumstance.

I need this to complete in say 10 seconds or less. I CAN limit the search to a radius, say a few thousand meters. But there are still likely to be thousands of query features inside that radius, so I'd prefer to avoid any brute-force loop/measure the distance solution in favor of one that uses a spatial index, if possible.

Any suggestions?

4

Here I go answering my own question. It turns out that a brute-force "spatial query a search radius, loop over every feature in it, measure distance to each" approach is faster than I expected. This implementation finishes basically instantaneously with a search radius that contains hundreds of target features.

This has the advantage of using the existing spatial index for the spatial query, if not the loop, and it can also honor any definition query and selection set in the query feature layer.

private function DoWork(IPoint inputpoint, IFeatureLayer queryflayer, double searchradius)
{
    IdDistancePair result = FindNearest((IGeometry)inputpoint, queryflayer, searchradius);

   //Do stuff with the result...
}

/// <summary>
/// This does a spatial query within a search radius, then loops
/// over features in that area, checking distance between the source geometry and each feature. 
/// At least the spatial query can use the spatial index. It is surprisingly fast.
/// This will optionally honor any pre-existing selection in the query FeatureLayer - 
/// that is, it can do a 'select from' the existing selection, if any.
/// </summary>
private IdDistancePair FindNearest(IGeometry g, IFeatureLayer flayer, double searchradius, string where = "", bool honorselection = true)
{
    IFeatureCursor fcursor = DoSpatialQuery(flayer, g, searchradius, where, honorselection);

    IdDistancePair returnval = FindNearest(g, fcursor);
    Marshal.FinalReleaseComObject(fcursor);
    return returnval;
}

///<summary>
/// Creates a spatial query which performs a spatial search for 
/// features in the supplied feature class and has the option to also apply an attribute query via a where clause.
/// See http://edndoc.esri.com/arcobjects/9.2/NET/7b4b8987-a3f0-4954-980f-720e61965449.htm
/// By accepting an IFeatureLayer (instead of IFeatureClass) any definition query on the layer will be honored, and we have the option to honor any preexisting selection set.
/// </summary>
private IFeatureCursor DoSpatialQuery(IFeatureLayer flayer, IGeometry searchGeometry, double bufferdistance = 0.0, string where = "", bool honorselection = true, esriSpatialRelEnum spatialRelation = esriSpatialRelEnum.esriSpatialRelIntersects, bool recycling = true)
{
    ISpatialFilter spatialFilter = CreateSpatialFilter(flayer, searchGeometry, spatialRelation, bufferdistance, where);

    // perform the query and use a cursor to hold the results
    IFeatureCursor featureCursor = null;
    if (!honorselection)
        featureCursor = flayer.Search((IQueryFilter)spatialFilter, recycling);
    else
    {
        if (((IFeatureSelection)flayer).SelectionSet.Count > 0)
        {
            ICursor cursor = null;
            ((IFeatureSelection)flayer).SelectionSet.Search((IQueryFilter)spatialFilter, recycling, out cursor);
            featureCursor = (IFeatureCursor)cursor;
        }
        else
            featureCursor = flayer.Search((IQueryFilter)spatialFilter, recycling);
    }

    Marshal.FinalReleaseComObject(spatialFilter);

    return featureCursor;
}

/// <summary>
/// Loops over rows in a FeatureCursor, measuring the distance from 
/// the target Geometry to each feature in the FeatureCursor, to
/// find the closest one.
/// </summary>
private IdDistancePair FindNearest(IGeometry g, IFeatureCursor fcursor)
{
    int closestoid = -1;
    double closestdistance = 0.0;
    IFeature feature = null;
    while ((feature = fcursor.NextFeature()) != null)
    {
        double distance = ((IProximityOperator)g).ReturnDistance(feature.Shape);
        if (closestoid == -1 || distance < closestdistance)
        {
            closestoid = feature.OID;
            closestdistance = distance;
        }
    }
    return new IdDistancePair(closestoid, closestdistance);
}

/// <summary>
/// Utility function to create a spatial filter.
/// </summary>
private ISpatialFilter CreateSpatialFilter(IFeatureLayer flayer, IGeometry searchGeometry, esriSpatialRelEnum spatialRelation, double bufferdistance, string where)
{
    // create a spatial query filter
    ISpatialFilter spatialFilter = new SpatialFilterClass();
    spatialFilter.GeometryField = flayer.FeatureClass.ShapeFieldName;
    spatialFilter.SpatialRel = spatialRelation;

    // specify the geometry to query with. apply a buffer if desired
    if (bufferdistance > 0.0)
    {
        // Use the ITopologicalOperator interface to create a buffer.
        ITopologicalOperator topoOperator = (ITopologicalOperator)searchGeometry;
        IGeometry buffer = topoOperator.Buffer(bufferdistance);
        spatialFilter.Geometry = buffer;
    }
    else
        spatialFilter.Geometry = searchGeometry;

    //limit the fields included on the cursor to those needed 
    //(you don't need to include columns used in the where clause)
    spatialFilter.SubFields = string.Format("{0},{1}", flayer.FeatureClass.OIDFieldName, flayer.FeatureClass.ShapeFieldName);

    // apply the where clause if desired
    if (!string.IsNullOrEmpty(where))
        spatialFilter.WhereClause = where;

    return spatialFilter;
}
  • Looks like you've found a solution but I was wondering if you could attribute your points with some sort of region code so that you only even test against points in a region which would knock out a vast majority of initial points. For example you are doing nearest distance on points along a river within a specific catchment, you would not test across catchment boundaries so be pre-labelling the points with catchment ID's and filter for those before you do all the intensive distance calculations should improve processing time. Just an idea? – Hornbydd Oct 30 '13 at 23:21
  • sure! my thinking too - the above allows for use of a WHERE to limit rows. – MC5 Oct 31 '13 at 13:56
3

Try this approach: Stick all of your 5 layers that you want to compare to your point class in one IGeometryBag (a geometry bag can hold points, lines or polygons). Then stick all of your points in another geometry bag. For each point, get distance to every layer in geometry bag using IProximityOperator.ReturnDistance(). Take that and the FID and store them in a list. Once you have looped through them all, get the minimum distance from the list, and that FID is your nearest feature. This method is much faster than IIndexQuery because you don't have to build an index (~1-2 seconds is usually how long it takes to run over several thousand points/lines/polygons).

Here's my working code that I grabbed that does something similar using this method:

    ESRI.ArcGIS.Geometry.IGeometryBag pGeometryBag = new GeometryBagClass();
    ESRI.ArcGIS.Geometry.IGeometryCollection pGeometryCollection = (IGeometryCollection)pGeometryBag;

    ESRI.ArcGIS.Carto.IEnumLayer pEnumLayer; //get all layers in the map
    pEnumLayer = pDoc.FocusMap.Layers;


    ESRI.ArcGIS.Carto.ILayer pLayer; //create a layer to store looped values of point combobox layers
    pLayer = pEnumLayer.Next();

    while (pLayer != null)//Get all the features and stick them into the geometry collection
    {

        ESRI.ArcGIS.Carto.IFeatureLayer pFeatureLayer = (IFeatureLayer)pLayer;
        if (pLayer == statusforselectedlayer) 
        {
             ITable pTable = (ITable)pFeatureLayer.FeatureClass;
             List<int> OIDList = new List<int>();
             IQueryFilter queryFilter = new QueryFilterClass();
             queryFilter.SubFields = pFeatureLayer.FeatureClass.OIDFieldName;
             ICursor pCursor = pTable.Search(queryFilter, false);
             IRow pRow = null;
             while((pRow = pCursor.NextRow()) != null)
             {
                 OIDList.Add(pRow.OID);
             }

             IFeatureCursor blocksCursor = pFeatureLayer.FeatureClass.GetFeatures(OIDList.ToArray(), false);
             IFeature blockFeature = null;
             while ((blockFeature = blocksCursor.NextFeature()) != null)
             {
                 pGeometryCollection.AddGeometry(blockFeature.Shape);
             }

        }

    ISpatialIndex spatialIndex = (ISpatialIndex)pGeometryBag;
    spatialIndex.AllowIndexing = true;
    spatialIndex.Invalidate();

    ISpatialFilter spatialFilter = new SpatialFilterClass();
    spatialFilter.Geometry = pGeometryBag;

    ESRI.ArcGIS.Geometry.IEnumGeometry pEnumGeometry = (IEnumGeometry)pGeometryBag;
    pEnumGeometry.Reset();
    ESRI.ArcGIS.Geometry.IGeometry pGeometry = pEnumGeometry.Next();

    List<double> distances = new List<double>();

    while (pGeometry != null)
    {
        IProximityOperator pProximity = (IProximityOperator)pGeometry;         

        ESRI.ArcGIS.Geometry.IEnumGeometry pEnumGeomBag = (IEnumGeometry)pGeometryBag;
        ESRI.ArcGIS.Geometry.IGeometry pGeomBagIterable = pEnumGeomBag.Next();

        while (pGeomBagIterable != null)
        {

            distances.Add(pProximity.ReturnDistance(pGeomBagIterable));
            pGeomBagIterable = pEnumGeomBag.Next();
        }
        pGeometry = pEnumGeometry.Next();


    }

Here's the full explanation of what the code does, and I'll come back and edit it later for your situation when I have more time, but that should be enough information for you to derive the solution.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.