1.梳理、总结经纬度处理在Maxcompute平台上的实战应用，如bd09、gcj02、wgs84经纬度坐标系转换UDF函数注册与使用。
2.欢迎批评指正，跪谢一键三连！
1.参考代码
- 坐标系转换
bd09坐标系（百度坐标系）转gcj02坐标系（中国国家测绘局（g表示guojia国家，c表示cehui测绘，j表示ju局）制订的地理信息系统的坐标系统）gcj02坐标系（中国国家测绘局（g表示guojia国家，c表示cehui测绘，j表示ju局）制订的地理信息系统的坐标系统）转bd09坐标系（百度坐标系）bd09坐标系（百度坐标系）转地球坐标系（world geodetic system）地球坐标系（world geodetic system）转gcj02坐标系（中国国家测绘局（g表示guojia国家，c表示cehui测绘，j表示ju局）制订的地理信息系统的坐标系统）gcj02坐标系（中国国家测绘局（g表示guojia国家，c表示cehui测绘，j表示ju局）制订的地理信息系统的坐标系统）转地球坐标系（world geodetic system）地球坐标系（world geodetic system）转bd09坐标系（百度坐标系）
# coding:utf-8 from odps.udf import annotate import math @annotate("double,double,string,string->string") class CoordTransform(object): ''' function: 坐标系转换 input: lng,lat,src,dest parms: lng:经度 double lat:纬度 double src:输入坐标系 ['bd','gcj','wgs'] dest:输出坐标系 ['bd','gcj','wgs'] output：string: 'lng,lat' ''' def evaluate(self, lng, lat, src, dest): if src == 'bd' and dest == 'gcj': return ','.join([str(x) for x in GeoUtils.bd09togcj02(float(lng), lat)]) elif src == 'bd' and dest == 'wgs': # return GeoUtils.bd09towgs84(lng,lat) return ','.join([str(x) for x in GeoUtils.bd09towgs84(lng, lat)]) elif src == 'gcj' and dest == 'wgs': return ','.join([str(x) for x in GeoUtils.gcj02towgs84(lng, lat)]) elif src == 'gcj' and dest == 'bd': return ','.join([str(x) for x in GeoUtils.gcj02tobd09(lng, lat)]) elif src == 'wgs' and dest == 'bd': return ','.join([str(x) for x in GeoUtils.wgs84tobd09(lng, lat)]) if src == 'wgs' and dest == 'gcj': return ','.join([str(x) for x in GeoUtils.wgs84togcj02(lng, lat)]) class GeoUtils(): PI = 3.1415926535897932384626 a = 6378245.0 ee = 0.00669342162296594323 @staticmethod def bd09togcj02(bd_lng, bd_lat): """bd09坐标系（百度坐标系）转gcj02坐标系（中国国家测绘局（g表示guojia国家，c表示cehui测绘，j表示ju局）制订的地理信息系统的坐标系统）""" x = bd_lng - 0.0065 y = bd_lat - 0.006 z = math.sqrt(x * x + y * y) - 0.00002 * math.sin(y * GeoUtils.PI) theta = math.atan2(y, x) - 0.000003 * math.cos(x * GeoUtils.PI) gcj_lng = z * math.cos(theta) gcj_lat = z * math.sin(theta) return [gcj_lng, gcj_lat] @staticmethod def gcj02tobd09(gcj_lng, gcj_lat): """gcj02坐标系（中国国家测绘局（g表示guojia国家，c表示cehui测绘，j表示ju局）制订的地理信息系统的坐标系统）转bd09坐标系（百度坐标系）""" z = math.sqrt(gcj_lng * gcj_lng + gcj_lat * gcj_lat) + 0.00002 * math.sin(gcj_lat * GeoUtils.PI) theta = math.atan2(gcj_lat, gcj_lng) + 0.000003 * math.cos(gcj_lng * GeoUtils.PI) bd_lng = z * math.cos(theta) + 0.0065 bd_lat = z * math.sin(theta) + 0.006 return [bd_lng, bd_lat] @staticmethod def bd09towgs84(lng, lat): """bd09坐标系（百度坐标系）转地球坐标系（world geodetic system）""" return GeoUtils.gcj02towgs84(*GeoUtils.bd09togcj02(lng, lat)) @staticmethod def wgs84togcj02(lng, lat): """地球坐标系（world geodetic system）转gcj02坐标系（中国国家测绘局（g表示guojia国家，c表示cehui测绘，j表示ju局）制订的地理信息系统的坐标系统）""" dlat = GeoUtils.transformlat(lng - 105.0, lat - 35.0) dlng = GeoUtils.transformlng(lng - 105.0, lat - 35.0) radlat = lat / 180.0 * GeoUtils.PI magic = math.sin(radlat) magic = 1 - GeoUtils.ee * magic * magic sqrtmagic = math.sqrt(magic) dlat = (dlat * 180.0) / ((GeoUtils.a * (1 - GeoUtils.ee)) / (magic * sqrtmagic) * GeoUtils.PI) dlng = (dlng * 180.0) / (GeoUtils.a / sqrtmagic * math.cos(radlat) * GeoUtils.PI) mglat = lat + dlat mglng = lng + dlng return [mglng, mglat] @staticmethod def gcj02towgs84(lng, lat): """gcj02坐标系（中国国家测绘局（g表示guojia国家，c表示cehui测绘，j表示ju局）制订的地理信息系统的坐标系统）转地球坐标系（world geodetic system）""" dlat = GeoUtils.transformlat(lng - 105.0, lat - 35.0) dlng = GeoUtils.transformlng(lng - 105.0, lat - 35.0) radlat = lat / 180.0 * GeoUtils.PI magic = math.sin(radlat) magic = 1 - GeoUtils.ee * magic * magic sqrtmagic = math.sqrt(magic) dlat = (dlat * 180.0) / ((GeoUtils.a * (1 - GeoUtils.ee)) / (magic * sqrtmagic) * GeoUtils.PI) dlng = (dlng * 180.0) / (GeoUtils.a / sqrtmagic * math.cos(radlat) * GeoUtils.PI) mglat = lat + dlat mglng = lng + dlng return [lng * 2 - mglng, lat * 2 - mglat] @staticmethod def wgs84tobd09(lng, lat): """地球坐标系（world geodetic system）转bd09坐标系（百度坐标系）""" dlat = GeoUtils.transformlat(lng - 105.0, lat - 35.0) dlng = GeoUtils.transformlng(lng - 105.0, lat - 35.0) radlat = lat / 180.0 * GeoUtils.PI magic = math.sin(radlat) magic = 1 - GeoUtils.ee * magic * magic sqrtmagic = math.sqrt(magic) dlat = (dlat * 180.0) / ((GeoUtils.a * (1 - GeoUtils.ee)) / (magic * sqrtmagic) * GeoUtils.PI) dlng = (dlng * 180.0) / (GeoUtils.a / sqrtmagic * math.cos(radlat) * GeoUtils.PI) mglat = lat + dlat mglng = lng + dlng z = math.sqrt(mglng * mglng + mglat * mglat) + 0.00002 * math.sin(mglat * GeoUtils.PI) theta = math.atan2(mglat, mglng) + 0.000003 * math.cos(mglng * GeoUtils.PI) bd_lng = z * math.cos(theta) + 0.0065 bd_lat = z * math.sin(theta) + 0.006 return [bd_lng, bd_lat] @staticmethod def transformlat(lng, lat): 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 * GeoUtils.PI) + 20.0 * math.sin(2.0 * lng * GeoUtils.PI)) * 2.0 / 3.0 ret += (20.0 * math.sin(lat * GeoUtils.PI) + 40.0 * math.sin(lat / 3.0 * GeoUtils.PI)) * 2.0 / 3.0 ret += (160.0 * math.sin(lat / 12.0 * GeoUtils.PI) + 320 * math.sin(lat * GeoUtils.PI / 30.0)) * 2.0 / 3.0 return ret @staticmethod def transformlng(lng, lat): 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 * GeoUtils.PI) + 20.0 * math.sin(2.0 * lng * GeoUtils.PI)) * 2.0 / 3.0 ret += (20.0 * math.sin(lng * GeoUtils.PI) + 40.0 * math.sin(lng / 3.0 * GeoUtils.PI)) * 2.0 / 3.0 ret += (150.0 * math.sin(lng / 12.0 * GeoUtils.PI) + 300.0 * math.sin(lng / 30.0 * GeoUtils.PI)) * 2.0 / 3.0 return ret @staticmethod def outofchina(lng, lat): """ 判断是否在国内，不在国内不做处理，异常点位 :param lng: :param lat: :return: """ return not (lng > 73.66 and lng < 135.05 and lat > 3.86 and lat < 53.55)
