Я заинтересован в изучении средней ширины многоугольника, который представляет дорожное покрытие. У меня также есть осевая линия дороги как вектор (который иногда не точно в центре). В этом примере линия дороги - красная, а многоугольник - синий:

Один из подходов, основанных на грубой силе, который я придумал, заключается в том, чтобы буферизовать линию небольшими приращениями, пересечь буфер с сеткой рыболовной сети, пересечь полигон дороги с сеткой рыболовной сети, вычислить площадь пересечения для обеих мер пересечения и продолжать делать это до тех пор, пока ошибка маленькая. Это грубый подход, и мне интересно, есть ли более элегантное решение. Кроме того, это скрыло бы ширину большой дороги и небольшой дороги.

Меня интересует решение, использующее программное обеспечение ArcGIS 10, PostGIS 2.0 или QGIS. Я видел этот вопрос и загрузил инструмент Дэна Паттерсона для ArcGIS10, но я не смог рассчитать, что я хочу с ним.

Я только что обнаружил инструмент Minimum Bounding Geometry в ArcGIS 10, который позволяет мне создавать следующие зеленые полигоны:

Это похоже на хорошее решение для дорог, которые следуют за сеткой, но не будут работать иначе, поэтому я все еще заинтересован в любых других предложениях.

Часть проблемы состоит в том, чтобы найти подходящее определение «средней ширины». Некоторые из них являются естественными, но будут отличаться, по крайней мере, немного. Для простоты рассмотрим определения, основанные на свойствах, которые легко вычислить (что, например, исключит те, которые основаны на преобразовании средней оси или последовательностях буферов).

В качестве примера рассмотрим, что архетипическая интуиция многоугольника с определенной «шириной» представляет собой небольшой буфер (скажем, радиуса r с квадратными концами) вокруг длинной, довольно прямой ломаной (скажем, длины L ). Мы думаем о 2r = w как о его ширине. Таким образом:

• Его периметр P приблизительно равен 2L + 2w;

• Его площадь A приблизительно равна w L.

Ширина w и длина L могут быть восстановлены как корни квадратичного x ^ 2 - (P / 2) x + A; в частности, мы можем оценить

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

Когда вы уверены, что многоугольник действительно длинный и тощий, в качестве дальнейшего приближения вы можете взять 2L + 2w, чтобы получить 2L, откуда

• w (грубо) = 2A / P.

Относительная ошибка в этом приближении пропорциональна w / L: чем тоньше многоугольник, тем ближе w / L к нулю и тем лучше становится приближение.

Мало того, что этот подход чрезвычайно прост (просто разделите площадь по периметру и умножьте на 2), с любой формулой не имеет значения, как ориентирован многоугольник или где он расположен (потому что такие евклидовы движения не изменяют ни площадь, ни периметр).

Вы можете рассмотреть возможность использования любой из этих формул для оценки средней ширины для любых полигонов, которые представляют сегменты улиц. Ошибка, которую вы делаете в первоначальной оценке w (с помощью квадратичной формулы), возникает из-за того, что область A также содержит крошечные клинья на каждом изгибе исходной полилинии. Если сумма углов изгиба равна t радиан (это полная абсолютная кривизна ломаной), то на самом деле

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

• A = L w + Pi tw ^ 2.

Включите их в предыдущее (квадратная формула) решение и упростите. Когда дым рассеивается, вклад от термина кривизны t исчез! То, что первоначально выглядело как приближение, совершенно точно для несамопересекающихся буферов ломаной линии (с квадратами на концах). Поэтому для полигонов переменной ширины это разумное определение средней ширины.

Здесь я показываю небольшую оптимизацию в отношении решения @whuber и говорю о «ширине буфера», потому что это полезно для интеграции решения более общей задачи: есть ли обратная функция st_buffer, которая возвращает оценку ширины?

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;

Для этой проблемы, @celenius вопрос о ширине улицы , sw, решение

 sw = buffer_width(ST_Length(g1), g2)

где sw«средняя ширина», g1центральная линия g2, а улица g2- это полигон . Я использовал только стандартную библиотеку OGC, протестировал PostGIS и решил другие серьезные практические приложения с той же функцией buffer_width.

ДЕМОНСТРАЦИЯ

A2это площадь g2, L1длина центральной линии ( g1) of g2.

Предположим , что мы можем генерировать с g2помощью g2=ST_Buffer(g1,w), и это g1является прямой, так g2это прямоугольник с длиной L1и шириной 2*w, и

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

Это не та же формула @whuber, потому что здесь wполовина g2ширины прямоугольника ( ). Это хорошая оценка, но, как мы видим из тестов (ниже), она не является точной, и функция использует ее как подсказку, чтобы уменьшить g2площадь, и как окончательную оценку.

Здесь мы не оцениваем буферы с «endcap = square» или «endcap = round», которые нуждаются в сумме с A2 областью точечного буфера с тем же w.

СПИСОК ЛИТЕРАТУРЫ: на аналогичном форуме 2005 года В. Хубер объясняет подобные и другие решения.

ИСПЫТАНИЯ И ПРИЧИНЫ

Для прямых результаты, как и ожидалось, являются точными. Но для других геометрий результаты могут быть неутешительными. Основная причина, пожалуй, в том, что вся модель предназначена для точных прямоугольников или для геометрий, которые могут быть аппроксимированы «прямоугольником полосы». Вот «тестовый набор» для проверки пределов этого приближения (см. wfactorРезультаты выше).

 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;

ПОЛУЧЕННЫЕ РЕЗУЛЬТАТЫ:

С прямоугольниками (центральная линия - прямая):

         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

С ДРУГИМИ ГЕОМЕТРИЯМИ (сложенная осевая линия):

         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.

Об этом btypeсмотрите руководство ST_Buffer , с хорошими иллюстрациями и LINESTRING, используемыми здесь.

ВЫВОДЫ :

• оценка w_estimвсегда лучше чем w_near;
• для "почти прямоугольной" g2геометрии, все в порядке, любойwfactor
• для другой геометрии (рядом с «прямоугольными полосами») используйте ограничение wfactor=~0.01на 1% ошибки w_estim. До этого фактора используйте другой оценщик.

Осторожность и профилактика

Почему возникает ошибка оценки? Когда вы используете ST_Buffer(g,w)«модель прямоугольной полосы», вы ожидаете, что новая область, добавленная буфером ширины, wравна w*ST_Length(g)или w*ST_Perimeter(g)... Если нет, обычно с помощью наложений (см. Согнутые линии) или «стилизации», это когда оценка средней wнеисправности . Это главное сообщение тестов.

Чтобы обнаружить эту проблему на любом типе буфера , проверьте поведение генерации буфера:

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;

ПОЛУЧЕННЫЕ РЕЗУЛЬТАТЫ:

         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%

Если вы можете объединить данные многоугольника с данными центральной линии (пространственными или табличными средствами), просто сложите площади многоугольника для каждого выравнивания центральной линии и разделите на длину центральной линии.

Я разработал формулу для средней ширины многоугольника и поместил ее в функцию Python / ArcPy. Моя формула получена из (но существенно расширяет) самого простого понятия средней ширины, которое я видел в других местах; то есть диаметр круга, имеющего ту же площадь, что и ваш многоугольник. Однако в приведенном выше вопросе и в моем проекте меня больше интересовала ширина самой узкой оси. Кроме того, меня интересовала средняя ширина для потенциально сложных, невыпуклых форм.

Мое решение было:

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

То есть:

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

Функция:

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

Вот экспортированная карта со средней шириной (и некоторыми другими атрибутами геометрии для справки) по множеству фигур с использованием функции сверху:

Другое решение с приблизительной медиальной осью:

1. Рассчитать примерную среднюю ось многоугольника;
2. Получить длину приблизительной медиальной оси;
3. Получить расстояние от обоих концов оси до границы многоугольника;
4. Суммируйте длину оси и расстояния от шага 3 - это приблизительная длина многоугольника;
5. Теперь вы можете разделить площадь многоугольника на эту длину и получить среднюю ширину многоугольника.

Результат, безусловно, будет неправильным для тех многоугольников, где приблизительная медиальная ось не является единой непрерывной линией, поэтому вы можете проверить это перед шагом 1 и вернуться NULLили что-то еще.

Вот пример функции PostgreSQL (примечание: вам нужно установить расширения postgis и 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;

Недостаток:

Это решение не будет работать со случаями, когда многоугольник является почти прямоугольным, и человек может интуитивно определить его длину, но приблизительная средняя ось имеет небольшие ветви вблизи края, и, таким образом, алгоритм возвращает None.

Пример: