I have a a map of about 3000 polygons in ArcMap 10. I'm looking to find the distance between each one of them. I know how to do it using the lat and long coordinates of the centroid, but I'm looking for the shortest straight line distance from the closest edge of one polygon to the closest edge of the other polygon.
Any ideas?
Assuming your table is in geographic coordinates, if not, just remove the casts to geography.
CREATE TABLE mytable_distances AS
SELECT a.id, b.id, ST_Distance(a.geom::geography, b.geom::geography) as distance
FROM mytable a, mytable b;
Distance from A to B is same as B to A, and distance from A to A is zero, therefore a half matrix will save you some work.
IProximityOperator returns distance from the edge. The code below uses an azimuthal projection centered on the centroid of each polygon (should work with lines too). If the the polygons aren't too complex (or if you have a lot of memory) loading all the geometries into memory an projecting them would be faster. (This is not thoroughly tested).
public class Pair
{
public int Oid1;
public int Oid2;
public double Dist;
public static void TestGetDistances()
{
IWorkspaceFactory wsf = new ESRI.ArcGIS.DataSourcesGDB.FileGDBWorkspaceFactoryClass();
string path = @"C:\Program Files\ArcGIS\DeveloperKit10.0\Samples\data\Usa\USA.gdb";
var fws = wsf.OpenFromFile(path, 0) as IFeatureWorkspace;
IFeatureClass fc = fws.OpenFeatureClass("states");
var halfMatrix = Pair.GetPairs(fc);
}
/// <summary>
/// key is oid of each feature, value is pairs for features with smaller oids.
/// </summary>
/// <param name="fc"></param>
/// <returns></returns>
public static SortedList<int, List<Pair>> GetPairs(IFeatureClass fc)
{
ISpatialReferenceFactory3 srf = new SpatialReferenceEnvironmentClass();
IProjectedCoordinateSystem pcs =
srf.CreateProjectedCoordinateSystem((int)esriSRProjCSType.esriSRProjCS_WGS1984N_PoleAziEqui);
var outList = new SortedList<int, List<Pair>>();
IFeatureCursor fCur = fc.Search(null, true);
IFeature f;
while ((f = fCur.NextFeature()) != null)
{
var pairs = GetDistances(f, pcs);
Debug.Print("{0} has {1} pairs", f.OID, pairs.Count);
outList.Add(f.OID, pairs);
}
System.Runtime.InteropServices.Marshal.FinalReleaseComObject(fCur);
return outList;
}
private static IPoint GetGCSCentroid(IGeometry geom)
{
if (geom.SpatialReference is IProjectedCoordinateSystem)
{
geom.Project(((IProjectedCoordinateSystem)geom.SpatialReference).GeographicCoordinateSystem);
}
IArea a = geom is IArea ? geom as IArea : geom.Envelope as IArea;
return a.Centroid;
}
/// <summary>
/// return a list of all other features whose OID is lesser than f1
/// </summary>
/// <param name="f1"></param>
/// <param name="pcs"></param>
/// <returns></returns>
private static List<Pair> GetDistances(IFeature f1, IProjectedCoordinateSystem pcs)
{
IPoint centroid = GetGCSCentroid(f1.ShapeCopy);
pcs.set_CentralMeridian(true, centroid.X);
((IProjectedCoordinateSystem2)pcs).LatitudeOfOrigin = centroid.Y;
var g1 = f1.ShapeCopy;
g1.Project(pcs);
var outList = new List<Pair>();
var fc = f1.Class as IFeatureClass;
var proxOp = g1 as IProximityOperator;
IFeatureCursor fCur = fc.Search(null, true);
IFeature f2 = null;
while ((f2 = fCur.NextFeature()) != null)
{
if (f2.OID < f1.OID)
{
var g2 = f2.ShapeCopy;
g2.Project(pcs);
outList.Add(new Pair()
{
Oid1 = f1.OID,
Oid2 = f2.OID,
Dist = proxOp.ReturnDistance(g2)
});
}
}
System.Runtime.InteropServices.Marshal.FinalReleaseComObject(fCur);
return outList;
}
}
The calculation isn't something that ArcGIS natively provides, but I gave a solution to the problem in Shapely which may help, or this thread on the old ArcGIS forums.
I think the near table tool would work for what you want:
Determines the distances from each feature in the input features to one or more nearby features in the near features, within the search radius. The results are recorded in the output table.
Just leave the search radius blank.
