前言
搞到了6万多点数据,点数据中存储了点所在的县名和坐标值,要转成矢量。涉及四百多个县,且经度上跨度较大,不在同一坐标系下,坐标值横坐标即有6位也有8位(有无带号都有),还存在无效的坐标。
文中数据点此下载。
处理思路
一、数据清洗(初步去除无效坐标)
去除无效坐标,X非6或者8位,Y非7位
此处根据X无效删除141条记录,X无效删除19条记录,剩余记录60892条,得到 2原始坐标删除无效坐标.xlsx 文件。
二、初步推理坐标系
已确定坐标使用 XIAN80 坐标系,推出中央经线或者带号即可。
8位可直接取前两位作为带号,6位就要根据县名去推测了。这里使用获取到的全国县界矢量去对应带号及中央经线,这里不区分3度带和6度带,因3度带中央经线包含了所有6度带的中央经线,更精细些,这里假设6位坐标均为3度带。
带号与中央经线换算
带号与中央经线换算:
3°带:x=3n
6°带:x=6n-3
经度与最近分带带号换算(x指要素中心点坐标):
3°带:n=int((x+1.5)/3)
6°带:n=int(x/6)+1
注意:这里根据全国县界矢量确定的中央经线为理论值,而坐标最初使用的中央经线有可能并非此值;匹配县名时注意国内是否存在一名多县的情况,如长安区、城关区、鼓楼区。
使用去除无效坐标点数据与县级行政界线进行挂接,取得各个点对应的投影。导出dbf删除多余字段,检查长安区(陕西)、城关区(甘肃)、鼓楼区(河南),分别对其修改,得到 3PointInfo.dbf 文件。
三、创建初步点数据,确定坐标推测错误的县
创建点数据,一般思路是根据不同坐标分别去创建,过于麻烦,这里利用 PointGeometry 类创建点时可指定坐标的特性用Python脚本实现。
脚本涉及到 arcpy 的 PointGeometry、POINT、SpatialReference类和 arcpy.da 游标的相关知识。
脚本(请勿放置在中文路径下运行)
import arcpy, os
arcpy.env.overwriteOutput = True
in_dbf = os.getcwd() + os.sep + "3PointInfo.dbf"
spatialReference = arcpy.SpatialReference(4610)
out_fc = os.getcwd() + os.sep + "TestCreatePoint.shp"
arcpy.CreateFeatureclass_management(os.getcwd(), "TestCreatePoint.shp", "POINT",
spatial_reference=spatialReference)
arcpy.AddField_management(out_fc, "XIAN", "TEXT", field_length=50)
arcpy.AddField_management(out_fc, "ZUOBIAO", "TEXT", field_length=100)
arcpy.AddField_management(out_fc, "X", "LONG")
arcpy.AddField_management(out_fc, "Y", "LONG")
cursor = arcpy.da.InsertCursor(out_fc, ["SHAPE@", "Id", "XIAN", "ZUOBIAO","X","Y"])
try:
with arcpy.da.SearchCursor(in_dbf, ["b1","b3","b4","b5","X_LENGTH","D_3_cen"]) as s_cursor:
i = 0
for row in s_cursor:
i += 1
print i,row[1], row[2], row[3]
if row[4] == 6:
spatialReference_prj = os.getcwd() + os.sep + "XIAN80" + os.sep +"Xian 1980 3 Degree GK CM " +
str(row[5]) + "E.prj"
elif row[4] == 8:
if int(str(row[2])[:2]) > 25:
spatialReference_prj = os.getcwd() + os.sep + "XIAN80" + os.sep + "Xian 1980 3 Degree GK Zone " +
str(row[2])[:2] + ".prj"
else:
spatialReference_prj = os.getcwd() + os.sep + "XIAN80" + os.sep + "Xian 1980 GK Zone " +
str(row[2])[:2] + ".prj"
spatialReference_temp = arcpy.SpatialReference(spatialReference_prj)
point = arcpy.Point(row[2], row[3])
point_geo = arcpy.PointGeometry(point, spatialReference_temp)
cursor.insertRow([point_geo, row[0], row[1],os.path.basename(spatialReference_prj).split('.')[0],row[2],row[3]])
del cursor
except arcpy.ExecuteError:
arcpy.GetMessages()
脚本运行过程中会发现新乡县点无法创建成功,查看坐标形式带号为24,而 XIAN80 3度带带号范围为25-45,6度带带号范围为13-23,并无24的带号,说明此X坐标也是无效,将其前两位改为38(b4+14000000),得到4PointInfo.dbf。
再次运行,没出啥问题,得到 TestCreatePoint.shp!
4、使用理论坐标点与县界相交,对于某县坐标完全偏离的人工推测投影带
观察矢量部分位于所有县级行政界线以外,垦利县,对其纠正,得到5PointInfo.dbf。
使用生成的点数据与县界相交,这里假定每个县所有坐标要么全部错误,要么全部正确。选出点县名称不等于县行政区域名称,且接近100%的点均不在理论县界范围内的县,确定这些县(取),对其假定中央经线进行校正。
如何调整中央经线?
要知道的是1度距离大概是110cosαkm(α指纬度值,这里只认定为110km)。
坐标为避免负值均会做500km东偏。
3度带范围则为:335000-665000,其中央为500000
用其坐标减去500000,若为正值则带号-1,中央经线-3;若为负值则带号+1,中央经线+3
注:为什么可以采用上面方法呢,因为几乎所有点均不在县级行政界线内,假定中央经线一定是错误的,故这种纠正即为最合理的,稍后用数据进行验证。
1、生成的点数据(TestCreatePoint.shp)与县界相交,得到 点与县级行政区域相交.shp ,分别对 所有XIAN 字段和 XIAN<>NAME 选择条件下 XIAN 字段进行汇总(6Sum_OutputAll.dbf 、7Sum_Output.dbf),挂接确定不在范围占比县共98个(8Sum_Output_JOIN_Sum_OutputAll.dbf)。
2、用5PointInfo.dbf挂接8Sum_Output_JOIN_Sum_OutputAll.dbf后,对中央经线进行纠正,
"8Sum_Output_JOIN_Sum_OutputAll.XIAN" IS NOT NULL AND "5PointInfo.b4" >500000筛选条件:
D_3_cen=D_3_cen-3
"8Sum_Output_JOIN_Sum_OutputAll.XIAN" IS NOT NULL AND "5PointInfo.b4" <500000筛选条件:
D_3_cen=D_3_cen+3
3、导出"8Sum_Output_JOIN_Sum_OutputAll.XIAN" IS NOT NULL筛选条件下5PointInfo.dbf中的点(9PointInfo.dbf),再次运行脚本,得到 TestCreatePoint1.shp。
检验TestCreatePoint1.shp正确性,10212个点中不在县级行政界线内的点有986个,占比为100%为二七区和新乡县(老城区、西工区点数较少不做分析),经观察此两县坐标格式有效值无效,不做校正。
4、最后利用5PointInfo.dbf中的点(5PointInfo.dbf),再次运行脚本,得到 TestCreatePoint2.shp。
最终检验TestCreatePoint2.shp存在6458个点未在县级行政界线内,此部分多数是两者所依据的行政界线有差异,部分点位于相邻县内,或者有些县用到了两套坐标系,还有一个原因是所使用的测试点为随机创建,坐标原本就位于相邻县内。
5、到此就完成了坐标到矢量的转换,此项数据只作为密度分析,6458个点也均匀分布,对作业精度影响不大。后续只需对明显的错误离散点进行删除操作即可。
技术总结
操作难点:
1、对于误差的把握及情况的假定
2、对于投影坐标系的理解
3、Python脚本的使用