我有兴趣检查代表道路表面的多边形的平均宽度。我还将道路中心线作为矢量(有时不完全位于中心)。在此示例中,道路中心线为红色,多边形为蓝色:



我想到的一种蛮力方法是通过以下方法缓冲线以较小的增量递增,将缓冲区与鱼网栅格相交,将道路多边形与鱼网栅格相交,计算两个交叉点度量的相交面积,并继续这样做直到误差很小。但是,这是一种粗略的方法,我想知道是否有更优雅的解决方案。此外,这将掩盖大道路和小道路的宽度。

我对使用ArcGIS 10,PostGIS 2.0或QGIS软件的解决方案感兴趣。我已经看到了这个问题,并下载了ArcGIS 10的Dan Patterson的工具,但无法使用它来计算。

我刚刚在ArcGIS 10中发现了最小边界几何工具,该工具可以启用我可以生成以下绿色多边形:



对于沿着网格的道路,这似乎是一个很好的解决方案,但无法正常工作,因此我仍然对任何其他建议。

评论

您是否在侧边栏排除了可能的解决方案?即gis.stackexchange.com/questions/2880/…显然标记为对可能重复的帖子的平凡答案

@DanPatterson我没有看到任何类似的问题(当然,很多问题都是相关的)。您是说我的问题被举报了吗?我不明白你的第二行。

@Dan这个相关的问题涉及对“宽度”的另一种解释(实际上,这种解释还不是很清楚)。那里的答复似乎集中在寻找最宽处的宽度,而不是平均宽度。

由于@whuber希望在这里集中讨论,结束另一个问题,我建议对该问题进行编辑,以使人们理解“矩形带的平均宽度的估计值”

@Peter:因为矩形是一个富裕的多边形,所以应该有更通用的标题。

#1 楼

问题的一部分是找到“平均宽度”的合适定义。有几种是自然的,但至少会略有不同。为了简单起见,请考虑基于易于计算的属性的定义(例如,将排除基于中间轴变换或缓冲区序列的那些定义)。

作为示例,请考虑具有确定的“宽度”的多边形的原型直觉是围绕一条较长的,相当笔直的折线(长度为L)的小缓冲区(半径为r的平方端)。我们认为2r = w为宽度。因此:


其周长P大约等于2L + 2w;
其面积A大约等于wL。

宽度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),无论使用哪种公式,多边形的方向或位置都无关紧要(因为这样的欧氏运动不会改变面积或周长)。 br />
您可以考虑使用这两个公式之一来估计代表街道段的任何多边形的平均宽度。您对w的原始估计(使用二次公式)产生的误差是由于区域A在原始折线的每个折弯处还包括微小的楔形而引起的。如果折角的总和为t弧度(这是折线的总绝对曲率),则实际上


P = 2L + 2w + 2 Pi tw和
A = L w + Pi tw ^ 2。

将它们插入以前的(二次方程式)解决方案并进行简化。当烟雾清除时,曲率项t的贡献就消失了!对于非自相交折线缓冲区(两端为正方形),最初看起来像是近似值是非常准确的。因此,对于宽度可变的多边形,这是平均宽度的合理定义。

评论


谢谢@whuber,这是一个很好的答案,它帮助我更加清楚地考虑它。

– djq
2012-2-14在0:09

@whuber:我正在写一篇论文,我需要对您在此处描述的方法给出适当的(“学术”)参考。你有这样的参考吗?这个措施有名字吗?如果没有,我可以用你的名字命名!那“ Huber的宽度度量”呢?

–朱利安
18 Mar 15 '18在10:34

@julien我没有任何引用。此格式可能适用:MISC {20279,TITLE = {计算多边形的平均宽度?},AUTHOR = {whuber(gis.stackexchange.com/users/664/whuber)},HOWPUBLISHED = {GIS},NOTE = {URL: gis.stackexchange.com/q/20279/664(版本:2013-08-13)},EPRINT = {gis.stackexchange.com/q/20279}、URL = {gis.stackexchange.com/q/20279}}

– hu
18-3-15在13:43



这是一个了不起的答案,所以很简单并且得到了很好的解释。谢谢!

–足球
19年11月16日在18:29

#2 楼

在这里,我没有显示有关@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是“平均宽度”,g1g2的中心线,而街道g2是POLYGON。我仅使用OGC标准库,并通过PostGIS进行了测试,并使用相同的buffer_width函数解决了其他严重的实际应用。
DEMONSTRATION
A2
假设我们可以通过g2生成L1,并且g1是直的,所以g2是长度为g2且宽度为g2=ST_Buffer(g1,w)的矩形,并且
    A2 = L1*(2*w)   -->  w = 0.5*A2/L1

不一样@whuber的公式,因为这里g1是矩形(g2)宽度的一半。这是一个很好的估计器,但是正如我们在测试中看到的那样(下图),它并不精确,该函数将其用作线索,以减小L1面积,并作为最终估计器。
这里我们不用“ endcap = square”或“ endcap = round”评估缓冲区,需要对具有相同2*w的点缓冲区的面积求和。
参考:在2005年的一个类似论坛中,W。Huber解释说这样的解决方案和其他解决方案。
测试和原因
对于直线,预期的结果是准确的。但是对于其他几何形状,结果可能令人失望。主要原因可能是,所有模型都是针对精确的矩形,或者是近似于“条形矩形”的几何形状。这里是一个“测试套件”,用于检查该近似值的限制(请参见上面结果中的w)。
 g2 

结果:
带有矩形(中心线是一条直线):
         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

其他GEOMETRIES(中心线折叠):
         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.

关于g2,请参见ST_Buffer指南,其中包含良好的ilustratins和LINESTRING。
结论:

A2的估计量始终是比w好;
对于“接近矩形”的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用于1%的w_estim上的错误。在此wfactor之前,请使用另一个估计器。宽度为w_near的缓冲区所增加的新区域大约为g2wfactor...。否则,通常是通过覆盖(请参见折线)或“样式”来估计平均wfactor=~0.01故障。这是测试的主要信息。
要在任何缓冲区中检测到此问题,请检查缓冲区生成的行为:
 w_estim 

结果:
         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%



#3 楼

如果可以将多边形数据连接到中心线数据(通过空间或表格方式),则只需将每个中心线路线的多边形面积相加,然后除以中心线的长度即可。

评论


确实如此!在这种情况下,我的中心线长度不一样,但是我总是可以将它们合并为一个,然后按每个多边形分割它们。

– djq
2012年2月14日,0:10

如果您的数据在postgreSQL / postGIS中,并且您有一个中心线和多边形的街道ID字段,则无需合并/分割,而使用聚合函数,您的答案就可以查询了。我在SQL方面很慢,否则我会举一个例子。让我知道这是否是您要解决的方式,如果需要,我会帮助您解决。

– Scro
2012-2-14在0:47



谢谢Scro,它现在不在PostGIS中,但是可以很快加载。我想我会先尝试@whuber的方法,但是我将其与PostGIS的结果进行比较(并感谢您提供SQL帮助,但是应该能够管理)。主要是想先弄清楚我的想法。

– djq
2012-2-14的2:11

+1这是一个适用的简单解决方案。

– hu
2012年2月14日下午6:01

#4 楼

我已经为多边形的平均宽度开发了一个公式,并将其放入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
 


这里是使用上面的函数导出的具有各种形状的平均宽度(以及其他一些几何属性供参考)的地图:



评论


如果简化表达式,则仅为面积/周长* 4。

–culebrón
17-2-21在22:14



谢谢,@culebrón。我只是想让概念更清晰,而不是公式简单,我甚至从未想过简化方程式。这应该为我节省一些处理时间。

–汤姆
17年2月22日在0:07

对于我遇到的问题,这是一个很好的解决方案,保留了大多数非道路多边形。

–Stonetip
20/12/7在23:20

#5 楼

具有近似中间轴的另一种解决方案:


计算多边形的近似中间轴;
获得近似中间轴的长度;
获得轴的两端到多边形的边界;
求和轴的长度和与第3步的距离-这是近似的多边形的长度;
现在,您可以将多边形的面积除以该长度,得到多边形的平均宽度。 >对于那些近似中间轴不是一条连续线的多边形,结果肯定是错误的,因此您可以在第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;


示例: