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
- 使用ArcScan提取道路中心线
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解释器的位置。
如果使用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)
4. 提取结果质量评价
整体效果不错,但肉眼检查可以发现还是有一些不尽人意的地方,例如交叉路口拓扑不对、靠得比较近的道路被连在一起、道路线条不平整(这个问题用arcscan会稍微好一些)等。如果有更好的方法欢迎交流。