Menghitung rata-rata lebar poligon?

41

Saya tertarik memeriksa lebar rata-rata poligon yang mewakili permukaan jalan. Saya juga memiliki garis tengah jalan sebagai vektor (yang kadang-kadang tidak persis di tengah). Dalam contoh ini, garis tengah jalan berwarna merah, dan poligon berwarna biru:

masukkan deskripsi gambar di sini

Salah satu pendekatan brute force yang saya pikirkan, adalah untuk menyangga garis dengan sedikit demi sedikit, memotong penyangga dengan jaring jaring, memotong jalan-poligon dengan jaring jaring, menghitung daerah berpotongan untuk kedua langkah persimpangan dan terus melakukan ini sampai kesalahannya kecil. Ini adalah pendekatan yang kasar, dan saya bertanya-tanya apakah ada solusi yang lebih elegan. Selain itu, ini akan menutupi lebar jalan besar dan jalan kecil.

Saya tertarik dengan solusi yang menggunakan perangkat lunak ArcGIS 10, PostGIS 2.0 atau QGIS. Saya telah melihat pertanyaan ini dan mengunduh alat Dan Patterson untuk ArcGIS10, tetapi saya tidak dapat menghitung apa yang saya inginkan dengannya.

Saya baru saja menemukan alat Geometri Batas Minimum di ArcGIS 10 yang memungkinkan saya menghasilkan poligon hijau berikut:

masukkan deskripsi gambar di sini

Ini sepertinya solusi yang baik untuk jalan yang mengikuti kisi-kisi, tetapi tidak akan bekerja sebaliknya, jadi saya masih tertarik dengan saran lain.

djq
sumber
Apakah Anda mengesampingkan kemungkinan solusi di bilah sisi? yaitu gis.stackexchange.com/questions/2880/… tampaknya ditandai sebagai jawaban yang sepele terhadap posting yang berpotensi rangkap
@ DanPatterson Saya tidak melihat pertanyaan seperti ini (banyak yang terkait, tentu saja). Apakah maksud Anda pertanyaan saya ditandai? Saya tidak mengerti baris kedua Anda.
djq
Pertanyaan terkait itu, @Dan, menyangkut interpretasi "lebar" yang berbeda (well, sebenarnya, interpretasinya tidak begitu jelas). Balasan di sana tampaknya berfokus pada menemukan lebar pada titik terlebar, daripada lebar rata-rata.
whuber
Ketika @whuber ingin memusatkan diskusi di sini, menutup pertanyaan lain, saya menyarankan agar pertanyaan itu diedit untuk orang-orang mengerti " perkiraan lebar rata-rata strip persegi panjang "
Peter Krauss
@ Peter: Karena strip persegi panjang adalah fortiori sebuah poligon, judul yang lebih umum harus berdiri.
whuber

Jawaban:

41

Bagian dari masalahnya adalah menemukan definisi yang sesuai dari "lebar rata-rata." Beberapa alami tetapi akan berbeda, setidaknya sedikit. Untuk kesederhanaan, pertimbangkan definisi berdasarkan properti yang mudah dihitung (yang akan mengesampingkan definisi yang didasarkan pada transformasi sumbu medial atau urutan buffer, misalnya).

Sebagai contoh, perhatikan bahwa intuisi pola dasar poligon dengan "lebar" yang pasti adalah penyangga kecil (katakanlah jari-jari r dengan ujung-ujung kuadrat) di sekitar polyline yang panjang dan cukup lurus (katakanlah panjang L ). Kami menganggap 2r = w sebagai lebarnya. Demikian:

  • Perimeter P- nya kira-kira sama dengan 2L + 2w;

  • Luasnya A kira-kira sama dengan w L.

Lebar w dan panjang L kemudian dapat diperoleh kembali sebagai akar kuadrat x ^ 2 - (P / 2) x + A; khususnya, kita dapat memperkirakan

  • w = (P - Sqrt (P ^ 2 - 16A)) / 4 .

Ketika Anda yakin poligon itu benar-benar panjang dan kurus, sebagai perkiraan lebih lanjut Anda dapat mengambil 2L + 2w untuk sama dengan 2L, dari mana

  • w (kasar) = 2A / P.

Kesalahan relatif dalam pendekatan ini sebanding dengan w / L: semakin kurus poligonnya, semakin dekat w / L ke nol, dan semakin baik perkiraannya.

Tidak hanya pendekatan ini sangat sederhana (hanya membagi area dengan perimeter dan kalikan dengan 2), dengan rumus apa pun itu tidak masalah bagaimana poligon berorientasi atau di mana ia berada (karena gerakan Euclidean seperti itu tidak mengubah area maupun perimeter).

Anda dapat mempertimbangkan menggunakan salah satu dari rumus ini untuk memperkirakan lebar rata-rata untuk setiap poligon yang mewakili segmen jalan. Kesalahan yang Anda buat dalam perkiraan awal w (dengan rumus kuadrat) terjadi karena area A juga mencakup irisan kecil di setiap belokan polyline asli. Jika jumlah sudut tikungan adalah t radian (ini adalah kelengkungan mutlak total dari polyline), maka benar-benar

  • P = 2L + 2w + 2 Pi tw dan

  • A = L w + Pi tw ^ 2.

Hubungkan ini ke solusi sebelumnya (rumus kuadrat) dan sederhanakan. Ketika asap hilang, kontribusi dari istilah kelengkungan t telah menghilang! Apa yang awalnya tampak seperti perkiraan sangat akurat untuk buffer polyline yang tidak memotong sendiri (dengan ujung kuadrat). Untuk poligon lebar variabel, oleh karena itu ini merupakan definisi yang wajar dari lebar rata-rata.

whuber
sumber
Terima kasih @whuber itu adalah jawaban yang bagus dan telah membantu saya untuk memikirkannya lebih jelas.
djq
@whuber: Saya sedang menulis makalah dan saya perlu memberikan referensi ('akademik') yang tepat untuk metode yang Anda jelaskan di sini. Apakah Anda punya referensi seperti itu? Apakah ukuran ini memiliki nama? Jika tidak, saya dapat memberi nama setelah Anda! Bagaimana dengan "ukuran lebar Huber"?
julien
@ Julien Saya tidak punya referensi. Format ini mungkin berfungsi: MISC {20279, TITLE = {Menghitung lebar rata-rata poligon?}, AUTHOR = {whuber ( gis.stackexchange.com/users/664/whuber )}, CARA DAPATKAN = {GIS}, CATATAN = {URL: gis.stackexchange.com/q/20279/664 (versi: 2013-08-13)}, EPRINT = { gis.stackexchange.com/q/20279 }, URL = { gis.stackexchange.com/q/20279 }}
whuber
1
Ini adalah jawaban yang luar biasa dan sangat sederhana dan dijelaskan dengan baik. Terima kasih!
Amball
18

Di sini saya menunjukkan sedikit optimasi tentang solusi @whuber, dan saya memasukkan istilah "lebar buffer", karena ini berguna untuk mengintegrasikan solusi dari masalah yang lebih umum: Apakah ada fungsi invers st_buffer, yang mengembalikan estimasi lebar?

CREATE FUNCTION buffer_width(
        -- rectangular strip mean width estimator
    p_len float,   -- len of the central line of g
    p_geom geometry, -- g
    p_btype varchar DEFAULT 'endcap=flat' -- st_buffer() parameter
) RETURNS float AS $f$
  DECLARE
    w_half float;
    w float;    
  BEGIN
         w_half := 0.25*ST_Area(p_geom)/p_len;
         w      := 0.50*ST_Area( ST_Buffer(p_geom,-w_half,p_btype) )/(p_len-2.0*w_half);
     RETURN w_half+w;
  END
$f$ LANGUAGE plpgsql IMMUTABLE;

Untuk masalah ini, pertanyaan @celenius tentang lebar jalan , sw, solusinya adalah

 sw = buffer_width(ST_Length(g1), g2)

di mana sw"lebar rata-rata", g1garis tengah g2, dan jalan g2adalah POLYGON . Saya hanya menggunakan pustaka standar OGC, diuji dengan PostGIS , dan menyelesaikan aplikasi praktis serius lainnya dengan fungsi buffer_width yang sama.

DEMONSTRASI

A2adalah luas g2, L1panjang garis tengah ( g1) dari g2.

Misalkan kita dapat menghasilkan g2dengan g2=ST_Buffer(g1,w), dan itu g1adalah lurus, jadi g2adalah persegi panjang dengan panjang L1dan lebar 2*w, dan

    A2 = L1*(2*w)   -->  w = 0.5*A2/L1

Ini bukan formula yang sama dengan @whuber, karena wini adalah setengah dari g2lebar persegi panjang ( ). Ini adalah penduga yang baik, tetapi seperti yang dapat kita lihat dengan tes (di bawah), tidak tepat, dan fungsinya menggunakannya sebagai petunjuk, untuk mengurangi g2area, dan sebagai penduga akhir.

Di sini kita tidak mengevaluasi buffer dengan "endcap = square" atau "endcap = round", yang membutuhkan penjumlahan A2 dari area buffer titik dengan yang sama w.

DAFTAR PUSTAKA: dalam forum serupa tahun 2005 , W. Huber menjelaskan solusi sejenis dan lainnya.

UJI DAN ALASAN

Untuk garis lurus hasilnya, seperti yang diharapkan, tepat. Tetapi untuk geometri lainnya hasilnya bisa mengecewakan. Alasan utamanya adalah, mungkin, semua model adalah untuk persegi panjang yang tepat, atau untuk geometri yang dapat diperkirakan menjadi "strip persegi panjang". Di sini "test kit" untuk memeriksa batas perkiraan ini (lihat wfactorhasil di atas).

 SELECT *, round(100.0*(w_estim-w)/w,1) as estim_perc_error
    FROM (
        SELECT btype, round(len,1) AS len, w, round(w/len,3) AS wfactor,
               round(  buffer_width(len, gbase, btype)  ,2) as w_estim ,
               round(  0.5*ST_Area(gbase)/len       ,2) as w_near
        FROM (
         SELECT
            *, st_length(g) AS len, ST_Buffer(g, w, btype) AS gbase
         FROM (
               -- SELECT ST_GeomFromText('LINESTRING(50 50,150 150)') AS g, -- straight
               SELECT ST_GeomFromText('LINESTRING(50 50,150 150,150 50,250 250)') AS g,
            unnest(array[1.0,10.0,20.0,50.0]) AS w
              ) AS t, 
             (SELECT unnest(array['endcap=flat','endcap=flat join=bevel']) AS btype
             ) AS t2
        ) as t3
    ) as t4;

HASIL:

DENGAN RECTANGLES (jalur sentral adalah LURUS LURUS):

         btype          |  len  |  w   | wfactor | w_estim | w_near | estim_perc_error 
------------------------+-------+------+---------+---------+--------+------------------
 endcap=flat            | 141.4 |  1.0 |   0.007 |       1 |      1 |                0
 endcap=flat join=bevel | 141.4 |  1.0 |   0.007 |       1 |      1 |                0
 endcap=flat            | 141.4 | 10.0 |   0.071 |      10 |     10 |                0
 endcap=flat join=bevel | 141.4 | 10.0 |   0.071 |      10 |     10 |                0
 endcap=flat            | 141.4 | 20.0 |   0.141 |      20 |     20 |                0
 endcap=flat join=bevel | 141.4 | 20.0 |   0.141 |      20 |     20 |                0
 endcap=flat            | 141.4 | 50.0 |   0.354 |      50 |     50 |                0
 endcap=flat join=bevel | 141.4 | 50.0 |   0.354 |      50 |     50 |                0

DENGAN GEOMETRI LAINNYA (garis tengah dilipat):

         btype          | len |  w   | wfactor | w_estim | w_near | estim_perc_error 
 -----------------------+-----+------+---------+---------+--------+------------------
 endcap=flat            | 465 |  1.0 |   0.002 |       1 |      1 |                0
 endcap=flat join=bevel | 465 |  1.0 |   0.002 |       1 |   0.99 |                0
 endcap=flat            | 465 | 10.0 |   0.022 |    9.98 |   9.55 |             -0.2
 endcap=flat join=bevel | 465 | 10.0 |   0.022 |    9.88 |   9.35 |             -1.2
 endcap=flat            | 465 | 20.0 |   0.043 |   19.83 |  18.22 |             -0.9
 endcap=flat join=bevel | 465 | 20.0 |   0.043 |   19.33 |  17.39 |             -3.4
 endcap=flat            | 465 | 50.0 |   0.108 |   46.29 |  40.47 |             -7.4
 endcap=flat join=bevel | 465 | 50.0 |   0.108 |   41.76 |  36.65 |            -16.5

 wfactor= w/len
 w_near = 0.5*area/len
 w_estim is the proposed estimator, the buffer_width function.

Tentang btypelihat panduan ST_Buffer , dengan ilustratin yang baik dan LINESTRING yang digunakan di sini.

KESIMPULAN :

  • penaksir w_estimselalu lebih baik daripada w_near;
  • untuk g2geometri "dekat ke persegi panjang" , ok, apa sajawfactor
  • untuk geometri lain (dekat dengan "strip segi empat"), gunakan batas wfactor=~0.01untuk kesalahan 1% w_estim. Sampai di wfactor ini, gunakan estimator lain.

Perhatian dan pencegahan

Mengapa kesalahan estimasi terjadi? Ketika Anda menggunakan ST_Buffer(g,w), Anda harapkan, dengan "model strip persegi panjang", bahwa area baru yang ditambahkan oleh buffer lebar wsekitar w*ST_Length(g)atau w*ST_Perimeter(g)... Ketika tidak, biasanya dengan overlay (lihat garis terlipat) atau dengan "gaya", adalah ketika estimasi wkesalahan rata-rata . Ini adalah pesan utama dari tes ini.

Untuk mendeteksi masalah ini di sembarang raja buffer , periksa perilaku generasi buffer:

SELECT btype, w, round(100.0*(a1-len1*2.0*w)/a1)::varchar||'%' AS straight_error,  
                 round(100.0*(a2-len2*2.0*w)/a2)::varchar||'%' AS curve2_error,
                 round(100.0*(a3-len3*2.0*w)/a3)::varchar||'%' AS curve3_error
FROM (
 SELECT
    *, st_length(g1) AS len1, ST_Area(ST_Buffer(g1, w, btype)) AS a1,
    st_length(g2) AS len2, ST_Area(ST_Buffer(g2, w, btype)) AS a2,
    st_length(g3) AS len3, ST_Area(ST_Buffer(g3, w, btype)) AS a3
 FROM (
       SELECT ST_GeomFromText('LINESTRING(50 50,150 150)') AS g1, -- straight
              ST_GeomFromText('LINESTRING(50 50,150 150,150 50)') AS g2,
              ST_GeomFromText('LINESTRING(50 50,150 150,150 50,250 250)') AS g3,
              unnest(array[1.0,20.0,50.0]) AS w
      ) AS t, 
     (SELECT unnest(array['endcap=flat','endcap=flat join=bevel']) AS btype
     ) AS t2
) as t3;

HASIL:

         btype          |  w   | straight_error | curve2_error | curve3_error 
------------------------+------+----------------+--------------+--------------
 endcap=flat            |  1.0 | 0%             | -0%          | -0%
 endcap=flat join=bevel |  1.0 | 0%             | -0%          | -1%
 endcap=flat            | 20.0 | 0%             | -5%          | -10%
 endcap=flat join=bevel | 20.0 | 0%             | -9%          | -15%
 endcap=flat            | 50.0 | 0%             | -14%         | -24%
 endcap=flat join=bevel | 50.0 | 0%             | -26%         | -36%

        waspada

Peter Krauss
sumber
13

Jika Anda dapat menggabungkan data poligon ke data garis tengah Anda (baik dengan cara spasial atau tabular), maka jumlahkan area poligon untuk setiap penyelarasan garis tengah dan bagi dengan panjang garis tengah.

Scro
sumber
itu benar! Dalam hal ini, garis tengah saya tidak sama panjangnya, tetapi saya selalu bisa bergabung sebagai satu, dan membaginya per poligon.
djq
Jika data Anda dalam postgreSQL / postGIS dan Anda memiliki bidang id jalan untuk centerlines dan poligon, maka tidak perlu digabung / dipisah, dan menggunakan fungsi agregat, jawaban Anda hanya tinggal pertanyaan saja. Saya lambat di SQL, atau saya akan mengirim contoh. Beri tahu saya jika ini yang akan Anda selesaikan, dan saya akan membantu mengatasinya (jika perlu.)
Scro
Terima kasih Scro, saat ini tidak dalam PostGIS, tetapi cukup cepat untuk memuat. Saya pikir saya akan mencoba pendekatan @ whuber pertama tetapi saya akan membandingkannya dengan hasil dari PostGIS (dan terima kasih atas tawaran bantuan SQL, tapi saya harus dapat dikelola). Terutama mencoba untuk mendapatkan pendekatan yang jelas di kepala saya dulu.
djq
+1 Ini adalah solusi sederhana yang bagus untuk keadaan di mana tersedia.
whuber
9

Saya telah mengembangkan formula untuk lebar rata-rata poligon dan memasukkannya ke dalam fungsi Python / ArcPy. Rumus saya berasal dari (tetapi secara substansial meluas) gagasan paling mudah tentang rata-rata lebar yang saya lihat dibahas di tempat lain; yaitu, diameter lingkaran yang memiliki area yang sama dengan poligon Anda. Namun, dalam pertanyaan di atas dan di proyek saya, saya lebih tertarik pada lebar sumbu tersempit. Selain itu, saya tertarik pada lebar rata-rata untuk bentuk yang berpotensi rumit dan tidak cembung.

Solusi saya adalah:

(perimeter / pi) * area / (perimeter**2 / (4*pi))
= 4 * area / perimeter

Itu adalah:

(Diameter of a circle with the same perimeter as the polygon) * Area / (Area of a circle with the same perimeter as the polygon)

Fungsinya adalah:

def add_average_width(featureClass, averageWidthField='Width'):
    '''
    (str, [str]) -> str

    Calculate the average width of each feature in the feature class. The width
        is reported in units of the feature class' projected coordinate systems'
        linear unit.

    Returns the name of the field that is populated with the feature widths.
    '''
    import arcpy
    from math import pi

    # Add the width field, if necessary
    fns = [i.name.lower() for i in arcpy.ListFields(featureClass)]
    if averageWidthField.lower() not in fns:
        arcpy.AddField_management(featureClass, averageWidthField, 'DOUBLE')

    fnsCur = ['SHAPE@LENGTH', 'SHAPE@AREA', averageWidthField]
    with arcpy.da.UpdateCursor(featureClass, fnsCur) as cur:
        for row in cur:
            perim, area, width = row
            row[-1] = ((perim/pi) * area) / (perim**2 / (4 * pi))
            cur.updateRow(row)

    return averageWidthField

Berikut adalah peta yang diekspor dengan lebar rata-rata (dan beberapa atribut geometri lain untuk referensi) di berbagai bentuk menggunakan fungsi dari atas:

masukkan deskripsi gambar di sini

Tom
sumber
4
Jika Anda menyederhanakan ekspresi, itu akan adil area / perimeter * 4.
culebrón
Terima kasih, @ culebrón. Saya mencari kejelasan konsep lebih dari kesederhanaan formula, dan saya bahkan tidak pernah berpikir untuk menyederhanakan persamaan. Ini akan menghemat waktu pemrosesan.
Tom
0

Solusi lain dengan perkiraan medial axis:

  1. Hitung perkiraan medial axis poligon;
  2. Dapatkan panjang perkiraan medial axis;
  3. Dapatkan jarak dari kedua ujung sumbu ke batas poligon;
  4. Jumlah panjang sumbu dan jarak dari langkah 3 - itu adalah perkiraan panjang poligon;
  5. Sekarang Anda dapat membagi area poligon dengan panjang ini dan mendapatkan lebar rata-rata poligon.

Hasil pasti akan salah untuk poligon-poligon di mana sumbu medial perkiraan bukan garis kontinu tunggal, sehingga Anda dapat memeriksanya sebelum langkah 1 dan kembali NULLatau sesuatu.

contoh

Berikut ini adalah contoh dari fungsi PostgreSQL (catatan: Anda perlu menginstal PostGIS dan postgis_sfcgal ekstensi):

CREATE FUNCTION ApproximatePolygonLength(geom geometry)
RETURNS float AS $$
    SELECT
        CASE
            /* in case when approximate medial axis is empty or simple line
             * return axis length
             */
            WHEN (ST_GeometryType(axis.axis) = 'ST_LineString' OR ST_IsEmpty(axis.axis))
                THEN axis_length.axis_length
                    + start_point_distance.start_point_distance
                    + end_point_distance.end_point_distance
            /* else geometry is too complex to define length */
            ELSE NULL
        END AS length
    FROM
        LATERAL (
            SELECT
                ST_MakeValid(geom) AS valid_geom
        ) AS valid_geom,
        LATERAL (
            SELECT
                /* `ST_LineMerge` returns:
                 *  - `GEOMETRYCOLLECTION EMPTY`, if `ST_ApproximateMedialAxis` is an empty line (i.e. for square);
                 *  - `LINESTRING ...`, if `ST_ApproximateMedialAxis` is a simple line;
                 *  - `MULTILINESTRING ...`, if `ST_ApproximateMedialAxis` is a complex line
                 *     that can not be merged to simple line. In this case we should return `NULL`.
                 */
                ST_LineMerge(
                    ST_ApproximateMedialAxis(
                        valid_geom.valid_geom
                    )
                ) AS axis
        ) AS axis,
        LATERAL (
            SELECT
                ST_Boundary(valid_geom.valid_geom) AS border
        ) AS border,
        LATERAL (
            SELECT
                ST_Length(axis.axis) AS axis_length
        ) AS axis_length,
        LATERAL (
            SELECT
                ST_IsClosed(axis.axis) AS axis_is_closed
        ) AS axis_is_closed,
        LATERAL (
            SELECT
                CASE WHEN axis_is_closed.axis_is_closed THEN 0
                ELSE
                    ST_Distance(
                        border.border,
                        CASE ST_GeometryType(axis.axis)
                            WHEN 'ST_LineString' THEN ST_StartPoint(axis.axis)
                            /* if approximate medial axis is empty (i.e. for square),
                             * get centroid of geometry
                             */
                            ELSE ST_Centroid(valid_geom.valid_geom)
                        END
                    )
                END AS start_point_distance
        ) AS start_point_distance,
        LATERAL (
            SELECT
                CASE WHEN axis_is_closed.axis_is_closed THEN 0
                ELSE
                    ST_Distance(
                        border.border,
                        CASE ST_GeometryType(axis.axis)
                            WHEN 'ST_LineString' THEN ST_EndPoint(axis.axis)
                            /* if approximate medial axis is empty (i.e. for square),
                             * get centroid of geometry
                             */
                            ELSE ST_Centroid(valid_geom.valid_geom)
                        END
                    )
                END AS end_point_distance
        ) AS end_point_distance;
$$ LANGUAGE SQL;

CREATE FUNCTION ApproximatePolygonWidth(geom geometry)
RETURNS float AS $$
    SELECT
        CASE
            WHEN length IS NULL THEN NULL
            ELSE area.area / length.length
        END AS width
    FROM
        (
            SELECT ApproximatePolygonLength(geom) AS length
        ) AS length,
        (
            SELECT
                ST_Area(
                    ST_MakeValid(geom)
                ) AS area
        ) AS area;
$$ LANGUAGE SQL;

Kerugian:

Solusi ini tidak akan bekerja dengan kasus ketika poligon hampir empat persegi panjang dan manusia secara intuite dapat menentukan panjangnya tetapi perkiraan medial axis memiliki cabang kecil di dekat tepi dan dengan demikian algoritma mengembalikan Tidak ada.

Contoh:

contoh rusak

Andrey Semakin
sumber