Interesuje mnie zbadanie średniej szerokości wielokąta reprezentującego powierzchnię drogi. Mam również linię środkową drogi jako wektor (która czasami nie jest dokładnie w środku). W tym przykładzie linia środkowa drogi jest czerwona, a wielokąt jest niebieski:

Jednym z podejść brutalnej siły, o którym myślałem, jest buforowanie linii małymi przyrostami, przecięcie bufora siatką kabaretki, przecięcie wielokąta drogowego siatką kabaretki, obliczenie przecinanego obszaru dla obu środków skrzyżowania i kontynuowanie tego aż do błąd jest niewielki. Jest to jednak prymitywne podejście i zastanawiam się, czy istnieje bardziej eleganckie rozwiązanie. Ponadto ukrywałoby to szerokość dużej drogi i małej drogi.

Interesuje mnie rozwiązanie wykorzystujące oprogramowanie ArcGIS 10, PostGIS 2.0 lub QGIS. Widziałem to pytanie i pobrałem narzędzie Dana Pattersona do ArcGIS10, ale nie byłem w stanie obliczyć, czego chcę.

Właśnie odkryłem narzędzie Minimalnej ograniczającej geometrii w ArcGIS 10, które pozwala mi produkować następujące zielone wielokąty:

To wydaje się być dobrym rozwiązaniem dla dróg, które biegną wzdłuż siatki, ale inaczej by nie działały, więc nadal jestem zainteresowany innymi sugestiami.

Częścią problemu jest znalezienie odpowiedniej definicji „średniej szerokości”. Niektóre są naturalne, ale będą się różnić, przynajmniej nieznacznie. Dla uproszczenia rozważ definicje oparte na właściwościach łatwych do obliczenia (co wykluczy te oparte na przykład na transformacji osi środkowej lub sekwencjach buforów).

Jako przykład weźmy pod uwagę, że archetypowa intuicja wielokąta o określonej „szerokości” jest małym buforem (powiedzmy o promieniu rz kwadratowymi końcami) wokół długiej, dość prostej polilinii (powiedzmy o długości L ). Myślimy o 2r = w jako jego szerokości. A zatem:

• Jego obwód P jest w przybliżeniu równy 2L + 2w;

• Jego powierzchnia A jest w przybliżeniu równa w L.

Szerokość w i długość L można następnie odzyskać jako pierwiastki kwadratowe x ^ 2 - (P / 2) x + A; w szczególności możemy oszacować

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

Gdy masz pewność, że wielokąt jest naprawdę długi i chudy, dla dalszego przybliżenia możesz wziąć 2L + 2 w, aby równać się 2L, skąd

• w (z grubsza) = 2A / P.

Błąd względny w tym przybliżeniu jest proporcjonalny do w / L: im chudszy wielokąt, tym bliżej w / L jest do zera, i tym lepsze jest przybliżenie.

Podejście to jest nie tylko niezwykle proste (wystarczy podzielić obszar przez obwód i pomnożyć przez 2), przy dowolnej formule nie ma znaczenia, w jaki sposób wielokąt jest zorientowany lub gdzie się znajduje (ponieważ takie ruchy euklidesowe nie zmieniają ani obszaru, ani obwód).

Możesz rozważyć użycie jednej z tych formuł w celu oszacowania średniej szerokości dowolnych wielokątów reprezentujących segmenty ulic. Błąd popełniany w pierwotnym oszacowaniu w (z formułą kwadratową) pojawia się, ponieważ obszar A zawiera również małe kliny na każdym zagięciu oryginalnej polilinii. Jeśli suma kątów zgięcia wynosi t radianów (jest to całkowita bezwzględna krzywizna polilinii), wtedy naprawdę

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

• A = L w + Pi tw ^ 2.

Podłącz je do poprzedniego rozwiązania (formuła kwadratowa) i uprość. Kiedy dym zniknie, zniknął wkład z terminu t krzywizny ! To, co pierwotnie wyglądało na przybliżenie, jest idealnie dokładne dla nie przecinających się buforów polilinii (z kwadratowymi końcami). W przypadku wielokątów o zmiennej szerokości jest to rozsądna definicja średniej szerokości.

Tutaj pokazuję małą optymalizację dotyczącą rozwiązania @whuber i mówię o „szerokości bufora”, ponieważ jest to przydatne do zintegrowania rozwiązania bardziej ogólnego problemu: czy istnieje funkcja odwrotna st_buffer, która zwraca oszacowanie szerokości?

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;

W tym problem, @celenius pytanie o szerokości ulicy , swrozwiązanie jest

 sw = buffer_width(ST_Length(g1), g2)

gdzie swjest „średnia szerokość”, g1środkowa linia g2, a ulica g2to POLYGON . Użyłem tylko standardowej biblioteki OGC, przetestowałem z PostGIS i rozwiązałem inne poważne praktyczne aplikacje z tą samą funkcją bufor_width.

DEMONSTRACJA

A2Jest to obszar g2, L1długość osi środkowej ( g1I) g2.

Przypuśćmy, że możemy wygenerować g2przez g2=ST_Buffer(g1,w), i że g1jest prosty, więc g2jest prostokątem o długości L1i szerokości 2*w, a

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

Nie jest to ta sama formuła @whuber, ponieważ tutaj wjest połowa g2szerokości prostokąta ( ). Jest to dobry estymator, ale jak widzimy w testach (poniżej), nie jest dokładny, a funkcja wykorzystuje go jako wskazówkę, aby zmniejszyć g2obszar i jako ostateczny estymator.

W tym przypadku nie oceniamy buforów za pomocą „endcap = kwadrat” lub „endcap = okrągły”, które potrzebują sumy do A2 obszaru bufora punktowego o tym samym w.

ODNIESIENIA: na podobnym forum w 2005 r. W. Huber wyjaśnia podobne i inne rozwiązania.

TESTY I POWODY

W przypadku linii prostych wyniki, zgodnie z oczekiwaniami, są dokładne. Ale w przypadku innych geometrii wyniki mogą być rozczarowujące. Być może głównym powodem jest to, że cały model dotyczy dokładnych prostokątów lub geometrii, które można przybliżyć do „prostokąta paska”. Oto „zestaw testowy” do sprawdzenia granic tego przybliżenia (patrz wfactorwyniki powyżej).

 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;

WYNIKI:

Z PROSTOKĄTAMI (linia środkowa to PROSTA LINIA):

         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

Z INNYMI GEOMETRIAMI (linia środkowa złożona):

         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.

O tym btypepatrz przewodnik ST_Buffer z dobrymi ilustracjami i zastosowanymi tutaj LINESTRINGAMI.

WNIOSKI :

• estymator w_estimjest zawsze lepszy niż w_near;
• dla g2geometrii „zbliżonych do prostokątnych” , jest w porządku, dowolnawfactor
• dla innych geometrii (w pobliżu „prostokątnych pasków”), użyj limitu wfactor=~0.011% błędu na w_estim. Do tego współczynnika użyj innego estymatora.

Ostrożność i zapobieganie

Dlaczego występuje błąd oszacowania? Kiedy używasz ST_Buffer(g,w), spodziewasz się, według „modelu paska prostokątnego”, że nowy obszar dodany przez bufor szerokości wwynosi około w*ST_Length(g)lub w*ST_Perimeter(g)… Kiedy nie, zwykle przez nakładki (patrz linie zagięte) lub „stylizację”, kiedy oszacowanie średniej wusterki . To jest główne przesłanie testów.

Aby wykryć ten problem w dowolnym królu bufora , sprawdź zachowanie generowania bufora:

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;

WYNIKI:

         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%

Jeśli możesz połączyć dane wielokąta z danymi linii środkowej (przestrzennie lub tabelarycznie), po prostu zsumuj obszary wielokątów dla każdego wyrównania linii środkowej i podziel przez długość linii środkowej.

Opracowałem wzór na średnią szerokość wielokąta i umieściłem go w funkcji Python / ArcPy. Moja formuła wywodzi się z (ale w znacznym stopniu rozszerza) najprostszego pojęcia średniej szerokości, które widziałem omówione gdzie indziej; to znaczy średnica koła mającego taki sam obszar jak twój wielokąt. Jednak w powyższym pytaniu i moim projekcie bardziej interesowała mnie szerokość najwęższej osi. Ponadto interesowała mnie średnia szerokość dla potencjalnie złożonych, niewypukłych kształtów.

Moim rozwiązaniem było:

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

To jest:

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

Funkcja to:

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

Oto wyeksportowana mapa ze średnią szerokością (i kilkoma innymi atrybutami geometrii dla odniesienia) w różnych kształtach przy użyciu funkcji z góry:

Inne rozwiązanie z przybliżoną osią środkową:

1. Oblicz przybliżoną oś środkową wielokąta;
2. Uzyskaj długość przybliżonej osi środkowej;
3. Uzyskaj odległość od obu końców osi do granicy wielokąta;
4. Suma długości osi i odległości od kroku 3 - jest to przybliżona długość wielokąta;
5. Teraz możesz podzielić powierzchnię wielokąta przez tę długość i uzyskać średnią szerokość wielokąta.

Wynik z pewnością będzie zły dla tych wielokątów, w których przybliżona oś środkowa nie jest pojedynczą linią ciągłą, więc możesz to sprawdzić przed krokiem 1 i wrócić NULLlub coś w tym rodzaju.

Oto przykład funkcji PostgreSQL (uwaga: musisz zainstalować rozszerzenia Postgis i Postgis_sfcgal ):

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;

Niekorzyść:

To rozwiązanie nie będzie działać w przypadkach, gdy wielokąt jest prawie prostokątny, a człowiek może intuicyjnie zdefiniować jego długość, ale przybliżona oś środkowa ma małe rozgałęzienia w pobliżu krawędzi, a zatem algorytm zwraca Brak.

Przykład: