|
|
本文记录下笔者基于高德地图开放的API使用地理解码和逆编码以及从高德地图里获取的要素从火星坐标系转WGS84坐标（粗糙的转换，采用从网上找到的转换算法）
1.地理解码及逆编码
地理解码
笔者定义了地理解码的函数，具体如下：
# 地理解码
def geocoding(address_list,key,city):
'''
address_list : 用于地理解码的地名列表,
key : 高德地图地理解码应用对应的密钥
city : 检索地名所在的地级市行政区划名称
return : 对应的地理解码结果
'''
print('此函数返回由地名，火星坐标系下地名对应的经度和纬度组成的pandas格式数据')
#导入必要库
import requests
import json
import pandas as pd
#定义请求地理解码工作的高德服务网址
url='https://restapi.amap.com/v3/geocode/geo'
#保存输出结果的变量
lon,lat,spatial = [],[],[]
for a in address_list:
try:
params = { 'key': key,
'address': a,
'city': city }
res = requests.get(url, params)
jd = json.loads(res.text)
coords = jd['geocodes'][0]['location']
lon.append(coords[:10])
lat.append(coords[-9:])
except:
lon.append('-1')
lat.append('-1')
for l1,l2 in zip(lon,lat):
spatial.append([float(l1),float(l2)])
spatial = pd.DataFrame(spatial,columns=['lon','lat'])
spatial['name']=address_list
spatial=spatial[['name','lon','lat']]
return spatial 函数中设置了三个参数，address_list为你要进行地理解码的地名列表，字符串格式存到列表，元组里；第二个参数key为你在高德地图开发者应用里设置用于web应用的密钥（读者可以自己去高德地图开放平台注册成为开发者，申请对应的Key，这个每天有限额)，如下所示：
图1： 高德地图开放平台中用于地理解码的
第三个为地名所属的地级市行政区划，用以匹配搜索。最后输出结果为pandas库的DataFrame格式，示例如下：
图2： 地理解码示例
逆地理解码
笔者封装的函数：
# 逆地理解码
def inverse_geocoding(lon,lat,key):
'''
lon : 要查询地名对应的经度
lat : 要查询地名对应的纬度
key : 高德服务对应的key
'''
import requests
import json
location_name = list()
def getname(lon,lat):
url = 'https://restapi.amap.com/v3/geocode/regeo?key={}&location={},{}'.format(key,lon,lat)
res = requests.get(url)
json_data = json.loads(res.text)
json_data=json_data['regeocode']['formatted_address']
return json_data
for l1,l2 in zip(lon,lat):
location_name.append(getname(l1,l2))
return location_name函数参数具体代码里有介绍，这里不再赘述，使用示例：
图3： 逆地理解码示例
火星坐标系转WGS84
笔者书写的函数如下：
# 火星坐标系转WGS84
def gcj02_to_wgs84(lon,lat):
'''
lon : 要转换的点的经度
lat : 要转换的点的纬度
return : 粗略转换为WGS84坐标的经纬度列表
'''
import numpy as np
lon_wgs84,lat_wgs84 = [],[]
def gcj02_to_wgs84_one(lng, lat):
"""
GCJ02(火星坐标系)转GPS84
:param lng:火星坐标系的经度
:param lat:火星坐标系纬度
:return:
"""
x_pi = 3.14159265358979324 * 3000.0 / 180.0
pi = 3.1415926535897932384626 # π
a = 6378245.0 # 长半轴
ee = 0.00669342162296594323 # 偏心率平方
def _transformlng(lng, lat):
lng, lat = np.array(lng), np.array(lat)
ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + \
0.1 * lng * lat + 0.1 * np.sqrt(np.fabs(lng))
ret += (20.0 * np.sin(6.0 * lng * pi) + 20.0 *
np.sin(2.0 * lng * pi)) * 2.0 / 3.0
ret += (20.0 * np.sin(lng * pi) + 40.0 *
np.sin(lng / 3.0 * pi)) * 2.0 / 3.0
ret += (150.0 * np.sin(lng / 12.0 * pi) + 300.0 *
np.sin(lng / 30.0 * pi)) * 2.0 / 3.0
return ret
def _transformlat(lng, lat):
lng, lat = np.array(lng), np.array(lat)
ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + \
0.1 * lng * lat + 0.2 * np.sqrt(np.fabs(lng))
ret += (20.0 * np.sin(6.0 * lng * pi) + 20.0 *
np.sin(2.0 * lng * pi)) * 2.0 / 3.0
ret += (20.0 * np.sin(lat * pi) + 40.0 *
np.sin(lat / 3.0 * pi)) * 2.0 / 3.0
ret += (160.0 * np.sin(lat / 12.0 * pi) + 320 *
np.sin(lat * pi / 30.0)) * 2.0 / 3.0
return ret
lng, lat = np.array(lng), np.array(lat)
dlat = _transformlat(lng - 105.0, lat - 35.0)
dlng = _transformlng(lng - 105.0, lat - 35.0)
radlat = lat / 180.0 * pi
magic = np.sin(radlat)
magic = 1 - ee * magic * magic
sqrtmagic = np.sqrt(magic)
dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * pi)
dlng = (dlng * 180.0) / (a / sqrtmagic * np.cos(radlat) * pi)
mglat = lat + dlat
mglng = lng + dlng
return lng * 2 - mglng, lat * 2 - mglat
for l1,l2 in zip(lon,lat):
lon_wgs84.append(gcj02_to_wgs84_one(l1,l2)[0])
lat_wgs84.append(gcj02_to_wgs84_one(l1,l2)[1])
return lon_wgs84,lat_wgs84使用示例：
图4： 坐标转换使用示例
<hr/>所使用的库均可使用pip install 安装，较为简单 |
|
发表于 2022-12-14 19:21:09
