Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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?

-Kevin

share|improve this question
add comment

4 Answers 4

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.

share|improve this answer
add comment

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;
    }
}
share|improve this answer
    
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 Nov 5 '10 at 10:39
add comment

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.

share|improve this answer
    
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
    
I clearly overstated my GIS knowledge. I'm going to need to learn to use POSTGIS. From a quick review, it seems very similar to SQL which I have a simplistic understanding of. –  Kevin Nov 5 '10 at 22:08
show 4 more comments

I think the near table tool would work for what you want, just leave the search radius blank.

share|improve this answer
    
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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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