坐标系转换
"""
坐标转换工具类
xll--->2021-05-19 developer
"""
import math
import pandas as pd
import numpy as np
from pyproj import Proj, transform, Transformer
from xxx.settings import BASE_DIR
import os
from warnings import simplefilter
simplefilter(action='ignore', category=FutureWarning)
simplefilter(action='ignore', category=DeprecationWarning)
# 坐标转换
class TransCoordinates:
"""
WGS84，BD09、GCJ02(火星坐标系)和UTM坐标系相互转换
"""
def __init__(self, longitude, latitude):
"""
初始化
:param longitude: 经度：float
:param latitude: 纬度：float
"""
self.longitude = longitude
self.latitude = latitude
self.pi = 3.1415926535897932384626 # π
self.x_pi = 3.14159265358979324 * 3000.0 / 180.0
self.a = 6378245.0 # 长半轴
self.ee = 0.00669342162296594323 # 偏心率平方
def out_of_china(self, longitude, latitude):
"""
判断经纬度是否在国内
:return: 在：True,不在：False
"""
return not (73.66 < longitude < 135.05 and 3.86 < latitude < 53.55)
def _transformlat(self, lng, lat):
"""
纬度转换
"""
ret = -100.0 + 2.0 * lng + 3.0 * lat + 0.2 * lat * lat + \
0.1 * lng * lat + 0.2 * math.sqrt(math.fabs(lng))
ret += (20.0 * math.sin(6.0 * lng * self.pi) + 20.0 *
math.sin(2.0 * lng * self.pi)) * 2.0 / 3.0
ret += (20.0 * math.sin(lat * self.pi) + 40.0 *
math.sin(lat / 3.0 * self.pi)) * 2.0 / 3.0
ret += (160.0 * math.sin(lat / 12.0 * self.pi) + 320 *
math.sin(lat * self.pi / 30.0)) * 2.0 / 3.0
return ret
def _transformlng(self, lng, lat):
"""
经度转换
"""
ret = 300.0 + lng + 2.0 * lat + 0.1 * lng * lng + \
0.1 * lng * lat + 0.1 * math.sqrt(math.fabs(lng))
ret += (20.0 * math.sin(6.0 * lng * self.pi) + 20.0 *
math.sin(2.0 * lng * self.pi)) * 2.0 / 3.0
ret += (20.0 * math.sin(lng * self.pi) + 40.0 *
math.sin(lng / 3.0 * self.pi)) * 2.0 / 3.0
ret += (150.0 * math.sin(lng / 12.0 * self.pi) + 300.0 *
math.sin(lng / 30.0 * self.pi)) * 2.0 / 3.0
return ret
def wgs84_to_gcj02(self):
"""
WGS84转GCJ02(火星坐标系)
"""
dlat = self._transformlat(self.longitude - 105.0, self.latitude - 35.0)
dlng = self._transformlng(self.longitude - 105.0, self.latitude - 35.0)
radlat = self.latitude / 180.0 * self.pi
magic = math.sin(radlat)
magic = 1 - self.ee * magic * magic
sqrtmagic = math.sqrt(magic)
dlat = (dlat * 180.0) / ((self.a * (1 - self.ee)) / (magic * sqrtmagic) * self.pi)
dlng = (dlng * 180.0) / (self.a / sqrtmagic * math.cos(radlat) * self.pi)
mglat = self.latitude + dlat
mglng = self.longitude + dlng
return [mglng, mglat]
def bd09_to_gcj02(self):
"""
百度坐标系转GCJ02(火星坐标系)
"""
if self.out_of_china(self.longitude, self.latitude): # 判断是否在国内
return [self.longitude, self.latitude]
x = self.longitude - 0.0065
y = self.latitude - 0.006
z = math.sqrt(x * x + y * y) - 0.00002 * math.sin(y * self.x_pi)
theta = math.atan2(y, x) - 0.000003 * math.cos(x * self.x_pi)
gg_lng = z * math.cos(theta)
gg_lat = z * math.sin(theta)
return [gg_lng, gg_lat]
def gcj02_to_bd09(self):
"""
GCJ02(火星坐标系)转百度坐标系
"""
z = math.sqrt(self.longitude * self.longitude + self.latitude * self.latitude) + 0.00002 * math.sin(
self.latitude * self.x_pi)
theta = math.atan2(self.latitude, self.longitude) + 0.000003 * math.cos(self.longitude * self.x_pi)
bd_lng = z * math.cos(theta) + 0.0065
bd_lat = z * math.sin(theta) + 0.006
return [bd_lng, bd_lat]
def gcj02_to_wgs84(self):
"""
火星坐标系(GCJ-02)转WGS84
"""
dlat = self._transformlat(self.longitude - 105.0, self.latitude - 35.0)
dlng = self._transformlng(self.longitude - 105.0, self.latitude - 35.0)
radlat = self.latitude / 180.0 * self.pi
magic = math.sin(radlat)
magic = 1 - self.ee * magic * magic
sqrtmagic = math.sqrt(magic)
dlat = (dlat * 180.0) / ((self.a * (1 - self.ee)) / (magic * sqrtmagic) * self.pi)
dlng = (dlng * 180.0) / (self.a / sqrtmagic * math.cos(radlat) * self.pi)
mglat = self.latitude + dlat
mglng = self.longitude + dlng
return [self.longitude * 2 - mglng, self.latitude * 2 - mglat]
def wgs84_to_bd09(self):
"""
WGS84转百度坐标系
"""
self.longitude, self.latitude = self.wgs84_to_gcj02()
return self.gcj02_to_bd09()
def bd09_to_wgs84(self):
"""
百度坐标系转WGS84
"""
self.longitude, self.latitude = self.bd09_to_gcj02()
return self.gcj02_to_wgs84()
def wgs84_to_utm(self):
"""
WGS84转UTM
"""
zone = int(self.longitude / 6) + 31
inProj = Proj(init="epsg:4326")
outProj = Proj(init="EPSG:326{}".format(zone))
UTM_Easting, UTM_Northing = transform(inProj, outProj, self.longitude, self.latitude)
return [UTM_Easting, UTM_Northing, zone]
def utm_to_wgs84(self, zone, hemisphere='N'):
"""
UTM转wgs84
"""
inProj = Proj(init="EPSG:326{}".format(zone))
outProj = Proj(init='EPSG:4326')
lon, lat = transform(inProj, outProj, self.longitude, self.latitude)
return [lon, lat]
# 读取excel,批量转换生成csv
class BatchTrans:
@staticmethod
def batch_bd09_to_wgs84(datafile):
data = pd.read_excel(datafile)
def process(x):
bd_lon = float(x['longitude'])
bd_lat = float(x['latitude'])
res = TransCoordinates(bd_lon, bd_lat).bd09_to_wgs84()
return res
data["wgs84"] = data[["longitude", "latitude"]].apply(lambda x: process(x), axis=1)
data.head()
file_path = os.path.join(BASE_DIR, 'static', 'trans_file', 'data_bd09_wgs84.csv')
data.to_csv(file_path, index=False)
return file_path
@staticmethod
def batch_bd09_to_utm(datafile):
data = pd.read_excel(datafile)
def process(x):
bd_lon = float(x['longitude'])
bd_lat = float(x['latitude'])
wgs = TransCoordinates(bd_lon, bd_lat).bd09_to_wgs84()
res = TransCoordinates(wgs[0], wgs[1]).wgs84_to_utm()
return res
data["UTM Easting"] = data[["longitude", "latitude"]].apply(lambda x: process(x)[0], axis=1)
data["UTM Northing"] = data[["longitude", "latitude"]].apply(lambda x: process(x)[1], axis=1)
data["Zone"] = data[["longitude", "latitude"]].apply(lambda x: process(x)[2], axis=1)
data.head()
file_path = os.path.join(BASE_DIR, 'static', 'trans_file', 'data_bd09_utm.csv')
data.to_csv(file_path, index=False)
return file_path
@staticmethod
def batch_wgs84_to_bd09(datafile):
data = pd.read_excel(datafile)
def process(x):
lon = float(x['longitude'])
lat = float(x['latitude'])
res = TransCoordinates(lon, lat).wgs84_to_bd09()
return res
data["bd09"] = data[["longitude", "latitude"]].apply(lambda x: process(x), axis=1)
data.head()
file_path = os.path.join(BASE_DIR, 'static', 'trans_file', 'data_wgs84_bd09.csv')
data.to_csv(file_path, index=False)
return file_path
@staticmethod
def batch_wgs84_to_utm(datafile):
data = pd.read_excel(datafile)
def process(x):
lon = float(x['longitude'])
lat = float(x['latitude'])
res = TransCoordinates(lon, lat).wgs84_to_utm()
return res
data["UTM Easting"] = data[["longitude", "latitude"]].apply(lambda x: process(x)[0], axis=1)
data["UTM Northing"] = data[["longitude", "latitude"]].apply(lambda x: process(x)[1], axis=1)
data["Zone"] = data[["longitude", "latitude"]].apply(lambda x: process(x)[2], axis=1)
data.head()
file_path = os.path.join(BASE_DIR, 'static', 'trans_file', 'data_wgs84_utm.csv')
data.to_csv(file_path, index=False)
return file_path
@staticmethod
def batch_utm_to_bd09(datafile):
data = pd.read_excel(datafile)
def process(x):
utm_x = float(x['utm_x'])
utm_y = float(x['utm_y'])
utm_zone = int(x['utm_zone'])
utm_hemisphere = str(x['utm_hemisphere'])
wgs = TransCoordinates(utm_x, utm_y).utm_to_wgs84(utm_zone, hemisphere=utm_hemisphere)
res = TransCoordinates(wgs[0], wgs[1]).wgs84_to_bd09()
return res
data["bd09"] = data[["utm_x", "utm_y", "utm_zone", "utm_hemisphere"]].apply(lambda x: process(x), axis=1)
data.head()
file_path = os.path.join(BASE_DIR, 'static', 'trans_file', 'data_utm_bd09.csv')
data.to_csv(file_path, index=False)
return file_path
@staticmethod
def batch_utm_to_wgs84(datafile):
data = pd.read_excel(datafile)
def process(x):
utm_x = float(x['utm_x'])
utm_y = float(x['utm_y'])
utm_zone = int(x['utm_zone'])
utm_hemisphere = str(x['utm_hemisphere'])
res = TransCoordinates(utm_x, utm_y).utm_to_wgs84(utm_zone, hemisphere=utm_hemisphere)
return res
data["wgs84"] = data[["utm_x", "utm_y", "utm_zone", "utm_hemisphere"]].apply(lambda x: process(x), axis=1)
data.head()
file_path = os.path.join(BASE_DIR, 'static', 'trans_file', 'data_utm_wgs84.csv')
data.to_csv(file_path, index=False)
return file_path
if __name__ == '__main__':
# wgs84转bd09
bd09 = TransCoordinates(122.8650033003745, 29.655506581424735).wgs84_to_bd09()
print(bd09)
# wgs84转utm
utm = TransCoordinates(121.802651, 31.150972).wgs84_to_utm()
print(utm)
# bd09转wgs84
wgs84 = TransCoordinates(122.87551044369413, 29.65919236613756).bd09_to_wgs84()
print(wgs84)
# utm转wgs84
wgs84_2 = TransCoordinates(646846.0009044164, 3865654.0003936305).utm_to_wgs84(50, hemisphere='N')
print(wgs84_2)
# 批量处理wgs84转bd09
aa = BatchTrans().batch_wgs84_to_bd09(r'D:\Project\yyr_tool\InternalTool\static\template_file\template.xls')
print(aa)