背景及思路
第一次用python操作 ArcGIS, 写下这篇博文记录一下。板块是房地产行业专有术语,很多维度特征都在板块层级讨论的。这不,要把上海市所有所有板块画出来,作为其他数据,如产业,人口数据的一个划分,作为最小研究地理空间单位。以前咨询过淘宝,把上海市板块用arcgis画出来大约需要150元人民币,恰好,手里头有板块边界的位置数据,即板块边界的经纬度,如下图所示
琢磨着,自己用python写个脚本,每个板块抽象成一个由多个折点构成的多边形,实现思路如下
- 读取板块边界经纬度数据
- 采用split()函数将经纬度数据分割,构成折点
- 将折点连起来形成封闭多边形,即板块
完整代码
# -*- coding: utf-8 -*-
"""
Project_name:draw plate for shanghai
@author: 帅帅de三叔
Created on Wed Oct 30 15:16:44 2019
"""
import sys
arcpy_path = [r'C:\Python27\ArcGIS10.6\Lib\site-packages',
r'C:\Program Files (x86)\ArcGIS\Desktop10.6\arcpy',
r'C:\Program Files (x86)\ArcGIS\Desktop10.6\bin',
r'C:\Program Files (x86)\ArcGIS\Desktop10.6\ArcToolbox\Scripts']
sys.path.extend(arcpy_path)
stdi,stdo,stde=sys.stdin,sys.stdout,sys.stderr
reload(sys) #通过import引用进来时,setdefaultencoding函数在被系统调用后被删除了,所以必须reload一次
sys.stdin,sys.stdout,sys.stderr=stdi,stdo,stde
sys.setdefaultencoding('utf-8')
import arcpy #导入地理处理模块
from arcpy import env #导入环境类
env.workspace=r"D:\python for ArcGIS\绘制上海板块"
env.overwriteOutput=True #是否开启复写
import pandas as pd #导入数据分析模块
plate_data=pd.read_excel("lat_lng.xlsx",sheetname='sheet0') #读取板块数据
rows,cols=plate_data.shape #数据框尺寸
lng_lat=plate_data[u'板块经纬度'] #经纬度数据
plate_name=plate_data[u'板块名称'] #板块名称
#factoryCode = arcpy.GetParameterAsText(4490) #WGS_1984_World_Mercator投影坐标系工厂代码4490,3395
#spatial_ref = arcpy.SpatialReference(factoryCode) #设置空间参考参数
spatial_ref = arcpy.SpatialReference('China Geodetic Coordinate System 2000') #China Geodetic Coordinate System 2000 or WGS 1984 World Mercator
polygonPoints=arcpy.Array() #用来存放构成多边形的折点
polygonGeometryList=[] #用来存放多边形几何对象组
for row in range(0,rows): #按行循环
points=lng_lat[row].split(";") #折点
for spot in points:
xy=spot.split(",") #折点经纬度
if len(xy)==2:
point=arcpy.Point() #几何对象,用来存放折点对象
point.id=row;point.X=float(xy[0]);point.Y=float(xy[1]) #转为点对象
polygonPoints.add(point) #构成一串折点
#print(point.id)
polygon=arcpy.Polygon(polygonPoints,spatial_ref,"","") #利用折点构造多边形带空间参考
polygonGeometryList.append(polygon) #把多边形追加到数组
polygonPoints.removeAll() #移除折点
arcpy.CopyFeatures_management(polygonGeometryList,r"D:\python for ArcGIS\绘制上海板块\plate_sh.shp","POLYGON")
结果截图
图中淡黄绿色就是一个一个板块。
后话
思路简单代码清晰, 所有的板块做成一个列表 polygonGeometryList ,然后一起画出来,也可以通过游标insert的办法插入。画出来的多边形没有其他属性,甚至没有板块名称和id,且两个相邻板块的公共分割线还存在不重合的情行,放大你会发现公共边会有空隙。