ArcGIS / arcpy 提取OSM道路中心线

酸奶太酸

已于 2023-07-09 17:37:28 修改

阅读量2.1k

点赞数 3

分类专栏： GIS 文章标签： arcgis

于 2023-07-09 17:30:48 首次发布

本文链接：https://blog.youkuaiyun.com/Yvelyne/article/details/131624074

版权

GIS 专栏收录该内容

6 篇文章

订阅专栏

1. 概述

目的
在做一些研究时会用到道路中心线，但OSM道路数据都是双线，需要一种提取道路中心线的方法，最好能通过代码实现批量处理。
思路
矢量路网数据转栅格，先膨胀以合并双线道路，然后细化、转矢量。
版本：ArcMap 10.5
参考：
- 使用ArcScan提取道路中心线
  https://gis.stackexchange.com/questions/29863/creating-centrelines-from-road-polygons-casings-using-arcgis-desktop
  https://gisgeography.com/how-to-vectorize-image-arcscan/
- ArcScan不能通过python调用，可以使用thin等ArcToolbox里的工具作为平替
  https://community.esri.com/t5/python-questions/access-arcscan-from-python/td-p/278660

2. ArcGIS操作步骤

数据准备
下载OSM上的主要道路（关键字highway，值trunk、primary、secondary、motorway）
投影（Project）
如果原始数据是WGS84坐标系，需要先投影到平面坐标系。
矢量转栅格（Polyline to Raster）
栅格的图形学运算比较方便，因此转成栅格处理。这里的Value field、Cell assignment type其实不重要，因为后面要转成二值图像。为了写代码方便，这里Value field使用highway字段。Cellsize设为20 m。
重分类（Reclassify）
将对应道路的像元值设为1，背景（NODATA）设为0。
膨胀（Expand）
矢量转为栅格后，双线道路之间还有空隙，因此通过expand让它们合并到一起。Number of cells可根据实际情况选择，这里使用了3。
收缩（Thin）
为了方便后续栅格转矢量，进行收缩。Shape for corners选择了SHARP（据说更适合道路之类的人工造物，但我看提取结果也不是很sharp），Maximum thickness使用默认值。
栅格转矢量（Raster to Polyline）
简化（Simplify Lines）
虽然前一步栅格转矢量的时候勾选了simplify选项，但得到的矢量线条还是不太平整，这里再次进行简化。Simplification Tolerance设为50 m，去掉Keep collapsed points。

3. arcpy代码

3.1 arcpy在哪里

参考这个文档，寻找ArcGIS的python解释器的位置。
pythn interpreter
如果使用PyCharm，可以按下图设置添加这个python解释器。
添加解释器

3.2 如何使用arcpy

参阅官方文档。ArcGIS的文档还是挺完整的，以Polyline to Raster工具为例，点击右下角的Tool help，即可看到该工具的参数说明和代码示例。
在这里插入图片描述

3.3 完整arcpy代码

import arcpy
from arcpy import env


def polylineToRaster(inFeatures, outRaster, cellSize):
    print "polylineToRaster: " + inFeatures
    # Set local variables
    valField = "highway"
    assignmentType = "MAXIMUM_COMBINED_LENGTH"
    priorityField = None

    # Execute PolylineToRaster
    arcpy.PolylineToRaster_conversion(inFeatures, valField, outRaster,
                                      assignmentType, priorityField, cellSize)


def reclassify(inRaster, outRaster):
    print "reclassify: " + inRaster
    # Set local variables
    reclassField = "highway"
    remap = arcpy.sa.RemapValue([["motorway", 1], ["primary", 1], ["secondary", 1], ["trunk", 1], ["NODATA", 0]])

    # Execute Reclassify
    outReclassify = arcpy.sa.Reclassify(inRaster, reclassField, remap)

    # Save the output
    outReclassify.save(outRaster)


def expand(inRaster, outRaster, numberCells):
    print "expand: " + inRaster
    # Set local variables
    zoneValues = [1]

    # Execute Expand
    outExpand = arcpy.sa.Expand(inRaster, numberCells, zoneValues)

    # Save the output
    outExpand.save(outRaster)


def thin(inRaster, outRaster):
    print "thin: " + inRaster
    # Execute Thin
    thinOut = arcpy.sa.Thin(inRaster, "ZERO", "FILTER", "SHARP")

    # Save the output
    thinOut.save(outRaster)
    return outRaster


def rasterToPolyline(inRaster, outLines, dangleTolerance):
    print "rasterToPolyline: " + inRaster
    # Set local variables
    backgrVal = "ZERO"
    field = "VALUE"

    # Execute RasterToPolygon
    arcpy.RasterToPolyline_conversion(inRaster, outLines, backgrVal,
                                      dangleTolerance, "SIMPLIFY", field)
    return outLines


def simplifyLines(inLines, outLines, tolerance):
    print "simplifyLines: " + inLines
    # Simplify all lines.
    arcpy.cartography.SimplifyLine(inLines,
                    outLines,
                    "POINT_REMOVE",
                    tolerance,
                    collapsed_point_option="NO_KEEP")


def pipeline(city):
    inFeatureName = city + "_main"
    rasterName = city + "_raster"
    polylineToRaster(inFeatureName, rasterName, 20)

    binaryRasterName = city + "_raster_bi"
    reclassify(rasterName, binaryRasterName)
    arcpy.Delete_management(rasterName)

    expandRasterName = city + "_expand"
    expand(binaryRasterName, expandRasterName, 3)
    arcpy.Delete_management(binaryRasterName)

    thinRasterName = city + "_thin"
    thin(expandRasterName, thinRasterName)
    arcpy.Delete_management(expandRasterName)

    centerlineName = city + "_centerline"
    rasterToPolyline(thinRasterName, centerlineName, 0)
    arcpy.Delete_management(thinRasterName)

    simpCenterlineName = city + "_centerline_simp"
    simplifyLines(centerlineName, simpCenterlineName, 50)
    arcpy.Delete_management(centerlineName)


if __name__ == "__main__":
    city = "beijing"
    # Set environment settings
    env.workspace = "D:/building/shp/arcgis.gdb"
    env.overwriteOutput = True  # overwrite existing output
    # Check out the ArcGIS Spatial Analyst extension license
    arcpy.CheckOutExtension("Spatial")

    pipeline(city)