arcgis desktop - ArcObjects - find nearest feature to a point - has to be FAST

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?

Answer

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...
}

/// 

/// 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.

/// 
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;
}


///

/// 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.
/// 
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;
}


/// 

/// Loops over rows in a FeatureCursor, measuring the distance from 
/// the target Geometry to each feature in the FeatureCursor, to
/// find the closest one.
/// 
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);
}

/// 

/// Utility function to create a spatial filter.
/// 
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;
}