Tag Info

Hot answers tagged

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


14

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 ...


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

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 ...


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 ...


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

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 ...


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

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.


4

GeographicLib has a python interface http://geographiclib.sourceforge.net/html/other.html#python 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 ...


4

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


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

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 ...


4

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 ...


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

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

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/


2

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 ...


2

Here is a full algebraic solution for two arbitrary small circles, using Cartesian coordinates to simplify the maths. It's assumed the earth is a spherical, but that's fine for my purposes.


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


2

This is a great question and one that I have a part answer too. The problem is a combination of the database and the rendering. Maybe there are smarter simpler answers, but here goes. The best solution is to create a MultiLineString that contains one linestring to the east of the dateline and one to the west. It also depends a bit on projection. I've gone ...


2

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.


1

FWIW here is some R code to explore the minimum distance to use between generated points. This uses great circle distance and intermediate points functions in geosphere, and reprojects to Mercator using rgdal. There's a hack to cut Antarctica a bit to stop it spreading to Infinity. ## background data library(maptools) data(wrld_simpl) dodge <- ...


1

It works for me this way: Create a text file with the following content: Nr;WKT 1;Linestring (0 0, 0 19113000) Create a custom CRS around the point of your interest: +proj=aeqd +lat_0=51 +lon_0=7 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs Note that I took the projection on a sphere, and my Linestring takes 3 times the radius. Set ...


1

great circles would be straight line from a single point if your projection is azimuthal centered on this point. so you can use the points2one plugin to create your line (basic trigonometry to have the coodinates), then reproject to your original coordinate system. Just make sure that you have dense lines (many vertices) so that you have a smooth curve after ...


1

CartoDB allows you to submit custom SQL queries so you can use the same commands as suggested in this post to create a great circle directly in CartoDB. For this demo I created two tables, with one point each (source and destination). Then I created a third empty table to hold the line string. How exactly you do this depends on your particular situation. ...


1

What you're looking for is the initial bearing (or forward azimuth), which if followed in a straight line along a great-circle arc will take you from the start point to the end point. Here is some simple JavaScript from this link: var y = Math.sin(dLon) * Math.cos(lat2); var x = Math.cos(lat1)*Math.sin(lat2) - ...



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