当前位置：首页 > news >正文

分享高德地图获取矢量面的方法和python脚本

news 来源：原创 2025/5/1 9:40:06

高德地图部分兴趣点是有矢量范围，本文分享一个简单方法，不使用高德API，直接通过浏览器F12开发工具获取坐标点数据，方式如箭头所示（打开右侧边框后重新点击矢量面，shape字段右键将字符串复制为JSON文本）

EDGE浏览器F12开发工具截图

python代码如下（自行修改121-123行，包括截图的坐标信息，位置名和保存路径）：

# -*- coding:utf-8 -*-
# ---------------------------------------------------------------------------
# Author: LGZ
# Created on: 
# Reference:
# coding:cp936 or coding:utf-8
# ---------------------------------------------------------------------------import logging, os, pprint
import pandas as pd
import numpy as np# import arcpy
# arcpy.env.overwriteOutput = True
# import itertools, random, matherrfile = r"ex.txt"
loggingfile = r"lg.txt"# logging.disable(logging.CRITICAL)  #禁用CRITICAL级别以下的日志记录
# 按此格式显示DEBUG级别以上的日志记录
# logging.basicConfig(filename=loggingfile, level=logging.DEBUG, format="%(asctime)s-%(levelname)s-%(message)s")
logging.basicConfig(level=logging.DEBUG, format="%(asctime)s-%(levelname)s-%(message)s")# pprint.pformat函数只能输入一个参数，logprint(*shape)中*shape为元组，将元组放进pprint.pformat
def logprint(*shape):return logging.debug(pprint.pformat(shape))# file =
# if os.path.isfile(file):  # 判断文件是否已存在，已存在则删除
#     os.unlink(file)
# if os.path.dirname(file) == "" :  # 判断是否没有完整路径
#     file = os.getcwd()+"//"+file  # 如没有则连接工作目录制造完整路径
# file_list=os.listdir(file)  # 根据需要的文件类型获取目录下文件名列表
# file_list=[i for i in file_list if i.endswith("xlsx") and \  
# not i.startswith("~$")]# save_fold = 
# os.makedirs(save_fold, exist_ok=True)  # 根据文件夹是否已存在制造文件夹import mathx_pi = float(3.14159265358979324 * 3000.0 / 180.0)
# //pai
pi = float(3.1415926535897932384626)
# //离心率
ee = float(0.00669342162296594323)
# //长半轴
a = float(6378245.0)
# //经度转换
def transformlat(lng:float, lat:float)-> float:ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + 0.1 * lng * lat + 0.2 * math.sqrt(abs(lng))ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 * math.sin(2.0 * lng * pi)) * 2.0 / 3.0ret += (20.0 * math.sin(lat * pi) + 40.0 * math.sin(lat / 3.0 * pi)) * 2.0 / 3.0ret += (160.0 * math.sin(lat / 12.0 * pi) + 320 * math.sin(lat * pi / 30.0)) * 2.0 / 3.0return ret# //纬度转换
def transformlng(lng:float, lat:float)-> float:ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + 0.1 * lng * lat + 0.1 * math.sqrt(abs(lng))ret += (20.0 * math.sin(6.0 * lng * pi) + 20.0 * math.sin(2.0 * lng * pi)) * 2.0 / 3.0ret += (20.0 * math.sin(lng * pi) + 40.0 * math.sin(lng / 3.0 * pi)) * 2.0 / 3.0ret += (150.0 * math.sin(lng / 12.0 * pi) + 300.0 * math.sin(lng / 30.0 * pi)) * 2.0 / 3.0return ret
# //国测局转84
def gcj02towgs84(lng:float, lat:float)-> list:dlat = transformlat(lng - 105.0, lat - 35.0)dlng = transformlng(lng - 105.0, lat - 35.0)radlat = lat / 180.0 * pimagic = math.sin(radlat)magic = 1 - ee * magic * magicsqrtmagic = math.sqrt(magic)dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi)mglat = lat + dlatmglng = lng + dlngreturn [lng * 2 - mglng, lat * 2 - mglat]def wgs84_to_gcj02(lng, lat):"""WGS84转GCJ02(火星坐标系):param lng:WGS84坐标系的经度:param lat:WGS84坐标系的纬度:return:"""dlat = transformlat(lng - 105.0, lat - 35.0)dlng =  transformlng(lng - 105.0, lat - 35.0)radlat = lat / 180.0 * pimagic = math.sin(radlat)magic = 1 - ee * magic * magicsqrtmagic = math.sqrt(magic)dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)dlng = (dlng * 180.0) / (a / sqrtmagic * math.cos(radlat) * pi)mglat = lat + dlatmglng = lng + dlngreturn [mglng, mglat]if __name__ == '__main__':"""脚本单独使用时运行以下内容""""""----------------------------------------------""""""---------------------PARA---------------------"""try:logprint("start of program")# TODOimport geopandas as gpdfrom shapely.geometry import Polygonimport matplotlib.pyplot as pltfrom pyproj import Transformertransformer_plane = Transformer.from_crs("EPSG:4490", "EPSG: 4526", always_xy=True)# 提供的坐标对coordinates_str = "113.254935,23.218678;113.254873,23.219087;113.255672,23.219361;113.255763,23.218744;113.255897,23.218481;113.256155,23.218113;113.256313,23.21819;113.256455,23.217946;113.256549,23.21798;113.256793,23.217596;113.256874,23.217598;113.256919,23.217588;113.256842,23.217423;113.256506,23.217179;113.256321,23.21729;113.256101,23.217366;113.255855,23.217411;113.25568,23.217502;113.255415,23.217748;113.255216,23.21802;113.255055,23.218308;113.254935,23.218678"name = "夏茅海口科创园"path = r"G:\PyProject\地理编码\其他工业园\\" + name + ".shp"# 分割坐标对coordinate_pairs = [pair.split(',') for pair in coordinates_str.split(';')]# 将坐标对转换为浮点数# x, y = transformer_plane.transform(gcj02towgs84(lng, lat)[0], gcj02towgs84(lng, lat)[1])coordinates = [(float(x), float(y)) for x, y in coordinate_pairs]WGScoordinates = [(gcj02towgs84(lng, lat)[0], gcj02towgs84(lng, lat)[1]) for (lng, lat) in coordinates]planecoordinates = [(transformer_plane.transform(lng, lat)[0], transformer_plane.transform(lng, lat)[1]) for (lng, lat) in WGScoordinates]# 创建多边形polygon = Polygon(planecoordinates)# 创建 GeoDataFramegdf = gpd.GeoDataFrame(geometry=[polygon])gdf["BZ"] = namegdf['MJ'] = gdf.geometry.area# gdf.to_file(r"G:\PyProject\地理编码\新建文件地理数据库.gdb", driver='OpenFileGDB', layer=name)gdf.to_file( path ,encoding='gbk')# assert a != 0, "a不能为0"# raise Exception("XX")except Exception as ex:# If an error occurred, print line number and error messageimport tracebackimport systb = sys.exc_info()[2]print(u"捕获到的异常信息是：{}".format(ex))  # 或者使用print("捕获到的异常信息是：",ex.args[0]或者str(ex)或者直接ex);3.9版本ex.message已不可用)print(u"捕获到的异常代码行号是：Line {0}".format(tb.tb_lineno))print(traceback.format_exc())  # 显示完整错误路径with open(errfile, "a") as err:err.write(traceback.format_exc())print(u"将traceback信息写入文件成功")# except arcpy.ExecuteError:# print(arcpy.GetMessages())# 无错误运行else后代码else:logprint("program success")# arcpy.AddMessage("program success")# 有没有错误均运行finally后代码finally:pass