5

I have a problem which draws my back in some project for some time now.

I'm basically looking to trap a polygon using x, y points drawn by some script I've written. lat, lon are the center GPS cords of the polygon and I'm looking for its surrounding polygon.

here is a part of my code in python:

def getcords(lat, lon, dr, bearing):
    lat2=asin(sin(lat)*cos(dr)+cos(lat)*sin(dr)*cos(bearing))
    lon2=lon+atan2(sin(bearing)*sin(dr)*cos(lat),cos(dr)-sin(lat)*sin(lat2))
    return [lat2,lon2]

my input goes like this:

  • lat, lon - are given in decimal degrees.
  • dr - is the angular computed by dividing the distance in miles by the earth's -radius(=3958.82)
  • bearing - between 0-360 degrees.

However for the input:

getcorsds1(42.189275, -76.85823, 0.5/3958.82, 30)

I get output: [-1.3485899508698462, -76.8576637627568], however [42.2516666666667, -76.8097222222222] is the right answer.

as for the angular distance, I calculate it simply by dividing the distance in miles by the earth's radius(=3958.82).

anybody?

Neuron
  • 5,141
  • 5
  • 38
  • 59
242Eld
  • 225
  • 2
  • 5
  • 13

3 Answers3

5

With geopy v2.0.0 (+ kilometers instead miles)

from geopy import Point                                                                                                                                                                       
from geopy.distance import geodesic                                                                                                                                                           
                                                                                                                                                                                              
distKm = 1                                                                                                                                                                                    
lat1 = 35.68096477080332                                                                                                                                                                      
lon1 = 139.76720809936523                                                                                                                                                                     
                                                                                                                                                                                              
print('center', lat1, lon1)                                                                                                                                                                   
print('north', geodesic(kilometers=distKm).destination(Point(lat1, lon1), 0).format_decimal())                                                                                                
print('east', geodesic(kilometers=distKm).destination(Point(lat1, lon1), 90).format_decimal())                                                                                                
print('south', geodesic(kilometers=distKm).destination(Point(lat1, lon1), 180).format_decimal())                                                                                              
print('west', geodesic(kilometers=distKm).destination(Point(lat1, lon1), 270).format_decimal())

result is

center 35.6809647708 139.767208099
north 35.6899775841, 139.767208099
east 35.680964264, 139.778254714
south 35.6719519439, 139.767208099
west 35.680964264, 139.756161485
Eric Duminil
  • 52,989
  • 9
  • 71
  • 124
tacosan
  • 174
  • 1
  • 4
  • 1
    **Don't use `geopy`!** It's no longer maintained and [have become problematic](https://github.com/geopy/geopy/pull/144). – not2qubit Jan 29 '18 at 06:10
  • 1
    @not2qubit That issues was resolved. Is there further evidence that we should stay clear of it? – Chris Jan 04 '19 at 20:11
  • For a better understanding, see [this answer](https://stackoverflow.com/a/48485675/1147688). – not2qubit Jan 05 '19 at 19:42
  • geopy had a [release](https://pypi.org/project/geopy/#history) six weeks ago, and eight last year. It would seem that sometime after @not2qubit 's comment, someone picked it up. – hlongmore Jan 30 '19 at 01:14
  • I'm glad that someone tries to maintain their own code base, but unfortunately some of the most basic info/documentation is still missing from the repo. With the result that the common user have no idea what they are calculating and how. – not2qubit Jan 30 '19 at 10:21
  • @not2qubit I updated the example above, with geodesics. I only did a quick test (starting at 0°N and °E) and moving 111km, and it seems to deliver plausible results. – Eric Duminil Jul 22 '20 at 19:49
4

Why don't you use nice libraries?

from geopy import Point
from geopy.distance import distance, VincentyDistance

# given: lat1, lon1, bearing, distMiles
lat2, lon2 = VincentyDistance(miles=distMiles).destination(Point(lat1, lon1), bearing)

For lat1, lon1, distMiles, bearing = 42.189275,-76.85823, 0.5, 30 it returns 42.1955489, -76.853359.

eumiro
  • 207,213
  • 34
  • 299
  • 261
  • I have tried to apply this solution in QGIS for points in EPSG:4326 but the result is not the proper one, it seems there is problem with the bearing, for example when I try bearing = 0 the new point is located in ~90º Any clue about the potential problem? I am using PyQGIS to create the layers –  Apr 17 '17 at 17:03
  • **Don't use `geopy`!** It's no longer maintained and [have become problematic](https://github.com/geopy/geopy/pull/144). – not2qubit Jan 29 '18 at 06:11
  • 3
    Is it not maintained? Many new versions since this answer. This syntax does not work with the current version (1.20.0). Point is no longer needed, Vincenty is deprectated, and a tuple with lat, long, and altitude is returned This works: `lat2,lon2,_ = distance(miles=distMiles).destination((lat1, lon1), bearing)` – erk Sep 19 '19 at 02:21
  • Is there a default we can use for 'bearing'? – Brendan Metcalfe Sep 25 '19 at 00:03
  • 1
    @Brendan - it doesn't make sense to have a default bearing, as you would then get lat/lon coordinates offset in an arbitrary direction; that's like saying to stand at one point on earth, walk a certain distance, then tell me your new position, without indicating what direction you should have walked. This answer has bearings along each N/S/E/W axis though: https://stackoverflow.com/a/40645383/181509 – EmpireJones Sep 28 '19 at 03:21
3

The sin and cos functions expect their arguments in radians, not in degrees. The asin and atan2 functions produce a result in radians, not in degrees. In general, one needs to convert input angles (lat1, lon1 and bearing) from degrees to radians using math.radians() and convert output angles (lat2 and lon2) from radians to degrees using math.degrees().

Note that your code has two other problems:

(1) It doesn't allow for travel across the 180-degrees meridian of longitude; you need to constrain your answer such that -180 <= longitude_degrees <= +180.

(2) If you are going to use this function extensively, you might like to remove the redundant calculations: sin(lat1), cos(dr), cos(lat1), and sin(dr) are each calculated twice.

John Machin
  • 81,303
  • 11
  • 141
  • 189