[postgis-users] Great circle as a linestring

[Dane Springmeyer]

Colin,

Okay, I've got a bit closer. The right way to do this is likely a  
function that computes the route between two points borrowing from  
some of the formulas here:
http://williams.best.vwh.net/avform.htm#Crs The postgis function that  
calculates the great circle distance is likely engineered this way,  
but we'll have to hear from the experts on that.

In the meantime, I'm close to a hack that might suffice for you  
without a custom function.

It is:

reproject your points > make a line > interpolate points along that  
line to get the curvature you are after for visual effect.

The projection you choose here will impact the route and how much it  
appears like a great circle.

Then plot the resulting line in openlayers with mapserver in WGS 84  
proj.

This seems to work aside from the problem that the route from hong  
kong to chicago cross the date line. I think there may be some  
workarounds for dateline stuff in openlayers, but I'm not sure.

Here is the code:

drop table if exists airports;
create table airports (id serial,the_geom geometry);
INSERT into airports (the_geom) values (geomfromtext('Point(-87.904841  
41.978603 0)',4326));
INSERT into airports (the_geom) values (geomfromtext('Point(113.914722  
22.308889 0)',4326));

drop table if exists eq_dist_airports;
SELECT id, ST_Transform(the_geom, 2163) as the_geom into  
eq_dist_airports from airports;

drop table if exists route;
select ST_Makeline(the_geom) as the_geom into route from  
eq_dist_airports;

drop table if exists points;
select line_interpolate_point(the_geom,generate_series(1,100)/100.0)  
as the_geom into points from route;

drop table if exists great_circle_wgs;
Select ST_Makeline(ST_Transform(the_geom, 4326)) as the_geom into  
great_circle_wgs from points;

And here is what I get viewing that in uDig:



On Feb 21, 2008, at 8:17 PM, Colin Wetherbee wrote:

> Dane Springmeyer wrote:
>> cool. I don't know much about how mapserver plays with postgis  
>> behind the scenes but I bet that reproject works because postgis is  
>> doing it (being ask by mapserver) and then mapserver can handle  
>> 4326 WGS 84. Am I understanding right that this is your SQL in your  
>> mapserver mapfile?
>
> I tried with SRIDs 4047 and 2163, which seem to be spherical-ish,  
> and those produce straight lines.  Your idea about MapServer not  
> interpreting 53027 as valid was probably correct.
>
> And, yes, that's the SQL from my map file.
>
> Actually, it's wrapped in stuff like...
>
> DATA "line FROM (query goes here) AS foo USING UNIQUE uid USING  
> SRID=xxx"
>
> Colin
>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20080222/5c9adc05/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: pastedGraphic.png
Type: image/png
Size: 83952 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20080222/5c9adc05/attachment.png>