Hot answers tagged

17

This is terrible code for general-purpose use because it can give erroneous results or even fail altogether for short distances. Use the Haversine Formula instead. (The formula on which your code is based converts two points on the sphere (not an ellipse) into their 3D Cartesian coordinates (xa,ya,za) and (xb,yb,zb) on the unit sphere and forms their dot ...


15

I would recommend checking out: Spherical: http://www.movable-type.co.uk/scripts/latlong.html Great-Circle: http://www.movable-type.co.uk/scripts/gis-faq-5.1.html


10

According to Wikipedia, Vincenty's formula is slower but more accurate: Vincenty's formulae are two related iterative methods used in geodesy to calculate the distance between two points on the surface of a spheroid, developed by Thaddeus Vincenty (1975a) They are based on the assumption that the figure of the Earth is an oblate spheroid, and ...


10

Although geodesics do look a little like sine waves in some projections, the formula is incorrect. Here is one geodesic in an Equirectangular projection. Clearly it is not a sine wave: (The background image is taken from http://upload.wikimedia.org/wikipedia/commons/thumb/e/ea/Equirectangular-projection.jpg/800px-Equirectangular-projection.jpg.) ...


7

The answers provided by others are a little more elegant, but here's an ultrasimple, somewhat unpythonic, bit of Python that provides the basics. The function takes two coordinate pairs and a user-specified number of segments. It yields a set of intermediate points along a great circle path. Output: text ready to write as KML. Caveats: The code does not ...


7

To create a great circle, you need either a custom azimuthal equidistant (aeqd) projection or a gnomonic (gnom) projection on one of your points. In those projections, the great circle is a straight line, which you can densify to have a curved line in other projections. If the great circle crosses the 180° E/W line or the poles, it might be useful to cut ...


6

It seems to be a straight line in whatever projection system pertains when it is created. After that, it is recalculated in each new projection, and the software trys to make it 'stratght'. this is quite noticeable near the poles: a square drawn round the pole in a polar azimuthal projection will invariably turn into a circle (that is, the formerly stratight ...


6

Sadly there isn't a geographic version of ST_ClosestPoint, so you will have to write your own function. There are two ways of calculating the nearest point of a great circle: spherical trigonometry, or 3D vector algebra. Luckily for you I have just written such a function for the latter method; I've not attempted the former because my spherical trig is ...


6

Take a look at this site, http://www.movable-type.co.uk/scripts/latlong.html If you implement it and you get the wrong answers you probably have the wrong units. I think for that site most operations are done in radians instead of decimal degrees.


6

So you know your two latitudes and longitudes, lets say You can calculate the cartesian co-ordinates for each: xa = (Cos(thisLat)) * (Cos(thisLong)); ya = (Cos(thisLat)) * (Sin(thisLong)); za = (Sin(thisLat)); xb = (Cos(otherLat)) * (Cos(otherLong)); yb = (Cos(otherLat)) * (Sin(otherLong)); zb = (Sin(otherLat)); And then calculate the great circle ...


6

If you're using geopy, then the great_circle and vincenty distances are equally convenient to obtain. In this case, you should almost always use the one that gives you the more accurate result, i.e., vincenty. The two considerations (as you point out) are speed and accuracy. Vincenty is two times slower. But probably in a real application the increased ...


6

You could calculate the Geodesics. Saying you want to show the geodesic from A to B, you could first calculate the distance and azimuth from A to B (inverse Geodesic problem) and then calculate points from A to several points between A and B (direct Geodesic problem). I have added a simple script in Python using GeographicLib just outputting the stuff in ...


5

Yes, use ST_Distance_Spheroid, or use a cast to the geography type SELECT ST_Distance('POINT(34.806480 32.083659)'::geography, 'POINT(34.8065365 32.0843373)'::geography); Hard to tell from the coordinates you are using, but make sure the longitude comes first, the latitude comes second in the POINT syntax.


5

GeographicLib has a python interface: This can computer geodesics on an ellipsoid (set flattening to zero to get great circles) and can generate intermediate points on a geodesic (see the "Line" commands in the sample). Here's how to print out intermediate points on the geodesic line from JFK to Changi Airport (Singapore): import sys; sys.path.append("/...


5

PostGIS computes its intersections (and distances (and areas)) differently depending on your column type. if the type is "geometry" the assumption is you want to work on a plane, so you get "straight line" distances / intersection tests if the type is "geography" the assumption is you want to work on a sphere, so you get great circle distances / tests ...


4

While an algebraic solution is possible if one assumes the earth is a sphere, we can still handle an ellipsoidal earth using Newton's method and a Esri's projection engine. The projection engine is a c style dll (pe.dll) and is bundled with the freely downloadable ArcGIS Explorer. I think the question could be rephrased as ... A plane is flying from ...


4

pyproj has the Geod.npts function that will return an array of points along the path. Note that it doesn't include the terminal points in the array, so you need to take them into account: import pyproj # calculate distance between points g = pyproj.Geod(ellps='WGS84') (az12, az21, dist) = g.inv(startlong, startlat, endlong, endlat) # calculate line string ...


4

geopy A Geocoding Toolbox for Python http://code.google.com/p/geopy/wiki/GettingStarted#Calculating_distances


4

You have to reproject every line into an azimutal equidistant projection based on one of its points, then densify the line. Since all points of your grid are connected to three lines, you can densify those three with the same projection. The Densify Geometry tool allows to densify only selected elements. See my example here on how to create great circles ...


4

In that particular projection - Platte Carre - i.e. just plotting latitude and longitude as if they are cartesian coordinates, your great circle line is a straight line. > gcIntermediate(c(168,-50),c(168,90),n=5,addStartEnd=TRUE) lon lat [1,] 168 -50.000000 [2,] 168 -26.666667 [3,] 168 -3.333333 [4,] 168 20.000000 [5,] 168 43.333333 [6,] ...


4

The principles in this blog post transfer over to general purpose PostGIS. http://blog.cartodb.com/jets-and-datelines/ Basically, use ST_Segmentize on geography, and a bit of magic to slice date-line crossing lines.


3

This is the PL/pgSQL version of MerseyViking's code. It also uses PostGIS geography Point and LineString types rather than a custom type to represent coordinates. CREATE OR REPLACE FUNCTION _point_to_cartesian(point geometry(Point), radius float, OUT x float, OUT y float, OUT z float) RETURNS RECORD AS $BODY$ DECLARE lon float; lat float; BEGIN ...


3

One-liner, assuming the POI points are stored in a "geography" column, you supply the ids of the two points and the search radius in meters: WITH line AS ( SELECT ST_MakeLine(p.geog::geometry, q.geog::geometry)::geography AS geog FROM pois p, pois q WHERE p.id = :id1 and q.id = :id2 ) SELECT p.name, p.id FROM pois p JOIN line ON ST_DWithin(p.geog, ...


3

Here are some links that might help: http://www.koders.com/python/fid0A930D7924AE856342437CA1F5A9A3EC0CAEACE2.aspx http://code.activestate.com/recipes/393241-calculating-the-distance-between-zip-codes/


3

The more general problem, posed for an ellipsoid of revolution, is considered in Section 8 of http://arxiv.org/abs/1109.4448v2 This gives solutions of the interception problem (the problem at hand) and the intersection problem using the ellipsoidal gnomonic projection. The same technique will apply to a sphere, of course.


3

See this link for code (using GeographicLib) to compute the cross track distance to a geodesic. This does the ellipsoidal calculation and is accurate to 20 nanometers or so. Also look at "Calculation 2. Distance from a point to a line" on this site. This calculates the cross track path with Javascript and displays the result.


3

If the wheels are pointing straight ahead, you are taking the shortest route toward some point directly ahead and thus following a great circle (or geodesic). If you wish to cross each meridian at exactly the same azimuth each time (and thus follow a rhumb line), you would have to gradually steer slightly more towards the meridian as you progress. Generally,...


3

Using a spherical model of the earth may give adequate accuracy and leads to simple fast calculations. Convert all coordinates into earth-centered (3D) cartesian coordinates. For example, the formula (cos(lon)*cos(lat), sin(lon)*cos(lat), sin(lat)) will do. (It uses a distance measure in which the earth's radius is one unit, which is convenient.) ...


2

You could use a similar hack to the ST_Intersection for geography. Would look like this: CREATE OR REPLACE FUNCTION st_closestpoint(geography, geography) RETURNS geography AS $$SELECT geography(ST_Transform(ST_ClosestPoint(ST_Transform(geometry($1), _ST_BestSRID($1,$2)),ST_Transform(geometry($2), _ST_BestSRID($1,$2)) ),4326)) $$ ...


2

How about transforming your starting coordinate to one on an azimuthal projection centered at the starting point. Then move at the appropriate azimuth for the desired distance, and transform back to geographic coordinates? A few resources: Wikipedia on Azimuthal Equidistant Projections Proj Syntax for the same



Only top voted, non community-wiki answers of a minimum length are eligible