Znaleźć minimalną odległość od krawędzi do krawędzi wielokątów za pomocą ArcGIS Desktop?

9

Mam mapę około 3000 wielokątów w ArcGIS 10. Szukam odległości między nimi. Wiem, jak to zrobić, używając współrzędnych szerokości i długości środka ciężkości, ale szukam najkrótszej linii prostej od najbliższej krawędzi jednego wielokąta do najbliższej krawędzi drugiego wielokąta. Jakieś pomysły?

Kevin
źródło

Odpowiedzi:

11

To niezły kawałek kodu, ale nie tak ładny jak (zakładając, że twoja tabela ma współrzędne geograficzne, jeśli nie po prostu usuń rzutowania do geografii)

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;

Czy wspominałem już, że przestrzenne bazy danych się zmieniają? Oni robią. Och, tak.

Paul Ramsey
źródło
Znajduje to odległość między najbliższymi wierzchołkami, ale nie samych krawędzi - nie wydaje się, że GEOS ujawnia tę dokładniejszą odpowiedź. Mimo to całkiem przydatny!
scw
1
Przepraszam scw, mylisz się na wiele sposobów. PostGIS ma natywne obliczenia odległości. GOES nie bierze w tym udziału. Po drugie, absolutnie podaje najbliższą odległość między krawędziami, nie tylko wierzchołkami zarówno w odległości geometrycznej, jak i przy obliczaniu odległości sferoidalnej typu geograficznego. Paul to napisał.
Nicklas Avén,
Aby zobaczyć to wizualnie dla geometrii, możesz użyć st_shortestline, która zwraca linię, która podaje odległość.
Nicklas Avén
1
Nik ma rację, zarówno w geometrii, jak i geografii funkcja odległości zwraca odległość między krawędziami. Na przykład wybierz st_distance („LINESTRING (0 0, 0 100)”, „LINESTRING (50 1, 51 1)”)
Paul Ramsey,
2
wow, przestrzenne bazy danych rządzą! obliczam odległość między zestawem ~ 8200 wielokątów a najbliższym sąsiadem w innym zestawie ~ 8400 wielokątów. w arcgis 10 narzędzie „generuj przy stole” o promieniu wyszukiwania 10000 m zajęło 1 godzinę i 15 minut (na czterordzeniowym komputerze stacjonarnym i7 3,4 GHz). to samo zapytanie w PostGIS zajęło tylko 3,5 minuty, i to na wolniejszym komputerze (dwurdzeniowy i7 macbook pro 2,7 GHz).
pistachionut
8

Odległość od A do B jest taka sama jak B do A, a odległość od A do A wynosi zero, dlatego pół macierzy pozwoli ci zaoszczędzić trochę pracy.

IProximityOperator zwraca odległość od krawędzi. Poniższy kod wykorzystuje odwzorowanie azymutalne wyśrodkowane na środku ciężkości każdego wielokąta (powinno również działać z liniami). Jeśli wielokąty nie są zbyt złożone (lub jeśli masz dużo pamięci) ładowanie wszystkich geometrii do pamięci, ich rzutowanie byłoby szybsze. (To nie jest dokładnie testowane).

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;
    }
}
Kirk Kuykendall
źródło
to niezły kawałek kodu. Nie wiedziałem o IproximityOperator i sam kodowałem coś w tym stylu (oczywiście wolniej)
George Silva,
2

Myślę, że narzędzie przy bliskiej tabeli działałoby tak, jak chcesz:

Określa odległości od każdej operacji w obiektach wejściowych do jednej lub więcej pobliskich obiektów w obiektach bliskich, w promieniu wyszukiwania. Wyniki są zapisywane w tabeli wyników.

Po prostu pozostaw pole wyszukiwania puste.

Jason Scheirer
źródło
Jest to rozwiązanie, które najpierw wypróbuję, ale wymaga ono poziomu licencji ArcInfo, aby odblokować narzędzie Generuj tabelę (Analiza).
PolyGeo