# Finding minimum edge to edge distance of polygons using ArcGIS Desktop?

I have a a map of about 3000 polygons in ArcGIS 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?

## 4 Answers

That is a nice piece of code, but not nearly as nice as (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;
``````

Have I mentioned that spatial databases rock? They do. Oh, they do.

• This will find the distance between the closest vertices, but not the edges themselves -- it doesn't seem that GEOS exposes this more precise answer. Still, quite handy! – scw Nov 5 '10 at 18:17
• Sorry scw, you are wrong in many ways. PostGIS has native distance calculations. GOES is not involved in that. Second, it absolutely gives the closest distance between edges, not only vertices both in geometry distance and in geography type spheroid distance calculation. Paul wrote it. – Nicklas Avén Nov 5 '10 at 19:20
• To see it visually for geometry you can use st_shortestline that returns the line that gives the distance. – Nicklas Avén Nov 5 '10 at 19:21
• Nik is right, in both geometry and geography the distance function returns the distance between the edges. For example, select st_distance('LINESTRING(0 0, 0 100)','LINESTRING(50 1, 51 1)') – Paul Ramsey Nov 5 '10 at 21:44
• wow, spatial databases do rock! i'm calculating the distance between a set of ~8200 polygons and the nearest neighbor in another set of ~8400 polygons. in arcgis 10, the 'generate near table' tool with a search radius of 10000 m took 1 hour and 15 minutes (on 3.4 GHz quad-core i7 desktop). the same query in PostGIS only took 3.5 minutes, and that was on a slower computer (a 2.7 GHz dual-core i7 macbook pro). – pistachionut Mar 9 '12 at 18:25

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;
}
}
``````
• this is a nice piece of code. I did not know about IproximityOperator, and ended up coding something somewhat like this myself (obviously it's slower) – George Silva Nov 5 '10 at 10:39

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.

• This is the solution I would try first but it needs an ArcInfo license level to unlock the Generate Near Table (Analysis) tool. – PolyGeo Dec 23 '10 at 1:48