[postgis-users] Great circle as a linestring

Dane Springmeyer blake at hailmail.net
Fri Feb 22 09:56:36 PST 2008

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