Encontrando distância mínima de ponta a ponta dos polígonos usando o ArcGIS Desktop?

9

Eu tenho um mapa de cerca de 3000 polígonos no ArcGIS 10. Estou procurando encontrar a distância entre cada um deles. Eu sei como fazê-lo usando as coordenadas lat e long do centróide, mas estou procurando a menor distância em linha reta da borda mais próxima de um polígono à borda mais próxima do outro polígono. Alguma ideia?

Kevin
fonte

Respostas:

11

Esse é um bom código, mas não tão bom quanto (supondo que sua tabela esteja em coordenadas geográficas, se não apenas remova as transmissões para a geografia)

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;

Eu mencionei que os bancos de dados espaciais balançam? Eles fazem. Oh eles fazem.

Paul Ramsey
fonte
Isso encontrará a distância entre os vértices mais próximos, mas não as arestas - não parece que o GEOS exponha essa resposta mais precisa. Ainda assim, bastante útil!
Scw #
11
Desculpe scw, você está errado de várias maneiras. O PostGIS possui cálculos de distâncias nativas. GOES não está envolvido nisso. Segundo, fornece absolutamente a distância mais próxima entre as arestas, não apenas os vértices, tanto na distância da geometria quanto no cálculo da distância do esferóide do tipo geografia. Paulo escreveu.
Nicklas Avén
Para visualizá-lo visualmente para geometria, você pode usar st_shortestline que retorna a linha que indica a distância.
Nicklas Avén
11
Nik está certo, tanto na geometria quanto na geografia, a função de distância retorna a distância entre as arestas. Por exemplo, selecione st_distance ('LINESTRING (0 0, 0 100)', 'LINESTRING (50 1, 51 1)')
Paul Ramsey
2
uau, bancos de dados espaciais são demais! Estou calculando a distância entre um conjunto de ~ 8200 polígonos e o vizinho mais próximo em outro conjunto de ~ 8400 polígonos. no arcgis 10, a ferramenta 'generate near table' com um raio de pesquisa de 10.000 m levou 1 hora e 15 minutos (na área de trabalho quad-core i7 de 3,4 GHz). a mesma consulta no PostGIS levou apenas 3,5 minutos, e foi em um computador mais lento (um macbook i7 dual-core de 2,7 GHz i7).
pistachionut
8

A distância de A a B é igual a B a A e a distância de A a A é zero, portanto, uma meia matriz poupará algum trabalho.

IProximityOperator retorna a distância da borda. O código abaixo usa uma projeção azimutal centralizada no centróide de cada polígono (também deve funcionar com linhas). Se os polígonos não forem muito complexos (ou se você tiver muita memória), carregar todas as geometrias na memória e projetá-las seria mais rápida. (Isso não foi totalmente testado).

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
fonte
este é um belo pedaço de código. Eu não sabia sobre IproximityOperator, e acabou de codificação algo mais ou menos assim me (obviamente é mais lento)
George Silva
2

Eu acho que a ferramenta Near Table funcionaria para o que você deseja:

Determina as distâncias de cada recurso nos recursos de entrada para um ou mais recursos próximos nos recursos próximos, dentro do raio de pesquisa. Os resultados são registrados na tabela de saída.

Apenas deixe o raio de pesquisa em branco.

Jason Scheirer
fonte
Esta é a solução que eu tentaria primeiro, mas precisa de um nível de licença do ArcInfo para desbloquear a ferramenta Gerar tabela próxima (análise).
PolyGeo