
# -*- coding: utf-8 -*-"""从图片提取坐标并生成点 Shapefile 文件坐标格式:度 分 秒(DMS)转十进制度(DD)"""import arcpyimport os# 设置覆盖输出arcpy.env.overwriteOutput = True# 输出路径output_dir = r"C:\Users\yl\Desktop\GIS"output_shp = os.path.join(output_dir, "p.shp")# 确保输出目录存在if not os.path.exists(output_dir): os.makedirs(output_dir) print("创建目录:" + output_dir)# 从图片中提取的坐标数据(度 分 秒格式)# 格式:经度度 经度分 经度秒,纬度度 纬度分 纬度秒,属性值coord_data = """108 7 46.333549, 30 42 24.997819, 2110001108 7 46.330161, 30 41 1.729522, 2110001108 7 19.072221, 30 41 1.733426, 2110001108 7 19.072120, 30 40 43.169068, 2110001108 6 41.499120, 30 40 43.167721, 2110001108 6 41.495237, 30 39 0.224086, 2110001107 57 4.279536, 30 39 0.228607, 2110001107 57 4.280175, 30 36 52.300680, 2110001107 53 30.199039, 30 36 52.296873, 2110001107 53 30.198787, 30 38 50.827047, 2110001107 54 1.600147, 30 38 50.825773, 2110001107 54 1.598039, 30 39 37.799964, 2110001108 1 23.846435, 30 39 37.795230, 2110001108 1 23.846020, 30 40 17.245728, 2110001108 2 53.997860, 30 40 17.252483, 2110001108 2 54.001286, 30 41 39.199217, 2110001108 4 10.801306, 30 41 39.201409, 2110001108 4 10.796124, 30 42 24.996779, 2110001108 7 46.333549, 30 42 24.997819, 2110001"""def dms_to_dd(degrees, minutes, seconds):"""度分秒转十进制度"""dd = abs(float(degrees)) + float(minutes)/60.0 + float(seconds)/3600.0if float(degrees) < 0: dd = -ddreturn ddprint("输出路径:" + output_shp)print()# 解析坐标数据points = []lines = coord_data.strip().split("\n")for line in lines: parts = line.strip().split(",")if len(parts) >= 3:# 经度lon_parts = parts[0].strip().split() lon_dd = dms_to_dd(lon_parts[0], lon_parts[1], lon_parts[2])# 纬度lat_parts = parts[1].strip().split() lat_dd = dms_to_dd(lat_parts[0], lat_parts[1], lat_parts[2])# 属性attr = parts[2].strip() points.append({"lon": lon_dd,"lat": lat_dd,"attr": attr,"lon_dms": parts[0].strip(),"lat_dms": parts[1].strip() })print("解析坐标点数:" + str(len(points)))print()# 创建 Shapefilespatial_ref = arcpy.SpatialReference(4326) # WGS84print("创建点 Shapefile(WGS84 坐标系)...")# 创建点要素类arcpy.management.CreateFeatureclass( out_path=output_dir, out_name="p.shp", geometry_type="POINT", spatial_reference=spatial_ref)# 添加字段arcpy.management.AddField(output_shp, "ID_Code", "TEXT", field_length=50)arcpy.management.AddField(output_shp, "Lon_DMS", "TEXT", field_length=50)arcpy.management.AddField(output_shp, "Lat_DMS", "TEXT", field_length=50)# 插入数据print("插入点数据...")with arcpy.da.InsertCursor(output_shp, ["SHAPE@", "ID_Code", "Lon_DMS", "Lat_DMS"]) as cursor:for pt in points:# 创建点几何point = arcpy.Point(pt["lon"], pt["lat"]) point_geometry = arcpy.PointGeometry(point, spatial_ref)# 插入记录cursor.insertRow([point_geometry, pt["attr"], pt["lon_dms"], pt["lat_dms"]])# 创建.cpg 文件解决中文编码cpg_path = output_shp.replace(".shp", ".cpg")with open(cpg_path, "w", encoding="utf-8") as f: f.write("CP936")print()print("完成!")print("输出文件:" + output_shp)print("点数量:" + str(len(points)))print()# 显示前 5 个点的信息print("前 5 个点坐标(十进制度):")for i, pt in enumerate(points[:5]): print(" " + str(i+1) + ": 经度=" + str(round(pt["lon"], 6)) + ", 纬度=" + str(round(pt["lat"], 6)) + ", 属性=" + pt["attr"])
夜雨聆风