高德地图部分兴趣点是有矢量范围，本文分享一个简单方法，不使用高德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, math
errfile = 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 math
x_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.0
ret += (20.0 * math.sin(lat * pi) + 40.0 * math.sin(lat / 3.0 * pi)) * 2.0 / 3.0
ret += (160.0 * math.sin(lat / 12.0 * pi) + 320 * math.sin(lat * pi / 30.0)) * 2.0 / 3.0
return 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.0
ret += (20.0 * math.sin(lng * pi) + 40.0 * math.sin(lng / 3.0 * pi)) * 2.0 / 3.0
ret += (150.0 * math.sin(lng / 12.0 * pi) + 300.0 * math.sin(lng / 30.0 * pi)) * 2.0 / 3.0
return 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 * pi
magic = math.sin(radlat)
magic = 1 - ee * magic * magic
sqrtmagic = 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 + dlat
mglng = lng + dlng
return [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 * pi
magic = math.sin(radlat)
magic = 1 - ee * magic * magic
sqrtmagic = 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 + dlat
mglng = lng + dlng
return [mglng, mglat]
if __name__ == '__main__':
"""脚本单独使用时运行以下内容"""
"""----------------------------------------------"""
"""---------------------PARA---------------------"""
try:
logprint("start of program")
# TODO
import geopandas as gpd
from shapely.geometry import Polygon
import matplotlib.pyplot as plt
from pyproj import Transformer
transformer_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)
# 创建 GeoDataFrame
gdf = gpd.GeoDataFrame(geometry=[polygon])
gdf["BZ"] = name
gdf['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 message
import traceback
import sys
tb = 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
扫一扫
1193
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
