• Arcgis-Tools_03混乱点转矢量


    前言

    搞到了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脚本实现。

    脚本涉及到 arcpyPointGeometryPOINTSpatialReference类和 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

    image-20200315122657576

    再次运行,没出啥问题,得到 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.dbf7Sum_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

    image-20200315133033172

    检验TestCreatePoint1.shp正确性,10212个点中不在县级行政界线内的点有986个,占比为100%为二七区和新乡县(老城区、西工区点数较少不做分析),经观察此两县坐标格式有效值无效,不做校正。

    4、最后利用5PointInfo.dbf中的点(5PointInfo.dbf),再次运行脚本,得到 TestCreatePoint2.shp

    image-20200315140550910

    最终检验TestCreatePoint2.shp存在6458个点未在县级行政界线内,此部分多数是两者所依据的行政界线有差异,部分点位于相邻县内,或者有些县用到了两套坐标系,还有一个原因是所使用的测试点为随机创建,坐标原本就位于相邻县内。

    5、到此就完成了坐标到矢量的转换,此项数据只作为密度分析,6458个点也均匀分布,对作业精度影响不大。后续只需对明显的错误离散点进行删除操作即可。

    技术总结

    操作难点:

    1、对于误差的把握及情况的假定

    2、对于投影坐标系的理解

    3、Python脚本的使用

  • 相关阅读:
    Asp.Net 之 前台绑定常用总结
    http无状态
    整理Py小demo
    RawURL
    整理的笔记
    单例模式
    DevExpress v18.1新版亮点——Analytics Dashboard篇(一)
    DevExpress v18.1新版亮点——Reporting篇(四)
    MyEclipse教程:使用UML创建模块库——第二部分(一)
    DevExpress v18.1新版亮点——Reporting篇(三)
  • 原文地址:https://www.cnblogs.com/bigmonk/p/12498034.html
Copyright © 2020-2023  润新知