wgs84计算2点距离(纬度不同、1度代表程度不同)

"""
LENGTH OF A DEGREE OF LATITUDE AND LONGITUDE BY COORDINATE
Calculates length of a degree of latitude and longitude based on geodetic 
meridian for any latitude and longitude position on an elipsoid without need of 
any external API or data. Constants based on elipsoid values used in WGS84, 
replicating calculation used b National Geospatial Agency (NDA) and CSGnet.
Formula is in a format that minimizes error at high latitudes by not dividing 
by cosines (like haversine calculations).
:param lat_degrees:		Latitude at postion for which you want degree lengths
:param units: 			Unit of measure for distnace returned. Default is 'km'
            			Other values: 'mi', 'nm', 'm', 'ft'
Result is a pair of floating point distances (lat_length, long_length)
	in unit distance 
Recommended for simple distance and geofence calculations below 88 degrees 
latutude and for distances and dimensions up to 200 km. Calculation has 
error rate <0.1% error at equator, under 3% (longitude) at the poles.  
"""
import sys
from math import sin, cos, pi, pow, degrees, radians

## Values for Earth, could swap in Mars, the Moon, etc.
EARTH_EQUATORIAL_RADIUS = 6378.137 # km
EARTH_POLAR_RADIUS = 6356.7523142 # km
EARTH_ECCENTRICITY_SQUARED = 0.00669437999014 

a = EARTH_EQUATORIAL_RADIUS
b = EARTH_POLAR_RADIUS
e_squared = EARTH_ECCENTRICITY_SQUARED

def format_latitude(latitude):
	# Cleanup latitudes that have gone over 90 degrees
	# Needed for correct radian conversion.
	# Then conversion to radians
	latitude = float(latitude or 0.0)
	if abs(latitude) > 90:
		latitude = ((latitude + 90) % 180) - 90
	return radians(latitude)

def format_longitude(longitude):
	# Cleanup longitudes that have gone over 180 degrees
	# Needed for correct radian conversion.
	# Then conversion to radians
	longitude = float(longitude or 0.0)
	if abs(longitude) > 180:
		longitude = ((longitude + 180) % 360) - 180
	return radians(longitude)

def get_degree_len(lat_degrees, units='km'):
	lat_radians = format_latitude(lat_degrees)

	# Meters (and Km) are based on EARTH_POLAR RADIUS
	# As such, meters and km distances are default of merdian calculation
	# Non-metric units will need to be adjusted based on conversion
	if units == 'km':
		conversion = 1.0 # km per km
	elif units == 'mi':
		conversion = 1/1.609 # mi per km
	elif units == 'm':
		conversion = 1000.0 # m per km
	elif units == 'ft':
		conversion = 5280/1.609 # ft per km
	elif units == 'nm':
		conversion = 1/1.852 # nm per km
	else:
		# Bail on invalid unity, can raise exception instead
		sys.exit("Error: %r is an invalid unit." % units)

	lat_length = (pi * a * (1.0 - e_squared)) / \
		(180.0 * pow(1 - (e_squared * pow(sin(lat_radians), 2)), 1.5)) \
		* conversion

	long_length = (pi * a * cos(lat_radians)) / \
		(180.0 * pow(1 - (e_squared * pow(sin(lat_radians), 2)), 0.5)) \
		* conversion
	return lat_length, long_length


### TESTING IT OUT FOR YOURSELF
# Input is latitude in degrees and unit measure
lat = float(input("Latitude (degrees)?: "))
unit = str(input("Units (km, m, mi, nm, ft)?: "))

# Showing how pair values are returned
dlat, dlong = get_degree_len(lat, unit)
print ("At %f degrees Latitude:" % lat)
print ("\t 1 Degree of Latitude is %f %s long" % (dlat, unit))
print ("\t 1 Degree of Longitude is %f %s long\n" % (dlong, unit))

 

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值