Wednesday 26 April 2017

carto - What is this PostGIS query doing to show great circle connections?


I'm trying to draw a flight map in CartoDB. My starting point is a single lat long point. We'll pretend it is this:


ST_GeomFromText('POINT(-71.060316 48.432044)', 4326)


And I want draw lines between that point and the points in my CartoDB table. I had a hard time even getting started, but someone pointed me to this map:


http://bl.ocks.org/andrewxhill/raw/8406976/


Which is fancier than what I need (I don't need to redraw everytime you click someplace) but includes the basic query:


SELECT ST_Transform(ST_Segmentize(ST_MakeLine(
ST_Transform(the_geom, 953027),
ST_Transform(
(SELECT ST_PointOnSurface(the_geom)
FROM lower_48_zips
WHERE cartodb_id = "+data.cartodb_id+"), 953027)), 100000), 4326)
the_geom FROM us_ski_areas WHERE mindist < 1000


Someone also pointed me towards an example which uses a very similar query:


UPDATE experimental.airroutes
SET the_geom =
(SELECT ST_Transform(ST_Segmentize(ST_MakeLine(
ST_Transform(a.the_geom, 953027),
ST_Transform(b.the_geom, 953027)
), 100000 ), 4326 )
FROM experimental.airports a, experimental.airports b
WHERE a.id = airroutes.source_id

AND b.id = airroutes.dest_id
);

But I'm having trouble pulling either of these apart. I'm decently conversant in MySQL but totally new to PostGIS and I would love help taking apart what is happening here and replacing source_id with a specific point.


I thought I could do something like


UPDATE my_table
SET the_geom =
(SELECT ST_Transform(ST_Segmentize(ST_MakeLine(
ST_Transform(my_table.the_geom, 953027),
ST_Transform(ST_GeomFromText('POINT(-71.060316 48.432044)', 4326), 953027)

), 100000 ), 4326 )
FROM my_table
);

But I know I need a "where" clause in there and I'm not even sure what it should be.



Answer



Short answer first, you are super close. Try this instead,


UPDATE my_table
SET the_geom =
ST_Transform(

ST_Segmentize(
ST_MakeLine(
ST_Transform(my_table.the_geom, 953027),
ST_Transform(CDB_LatLng(48.432044, -71.060316), 953027)
),
100000
),
4326
)


This would update your table from being points, to being a line from the original point to the point at -71, 48. I used the CDB_LatLng helper function. You didn't need to As_Text function, because that just turns geometries into human readable text.


Now the longer is that each of these functions has a documentation page, you can find those here,


ST_Transform: http://postgis.net/docs/ST_Transform.html


ST_Segmentize: http://postgis.net/docs/ST_Segmentize.html


ST_MakeLine: http://postgis.net/docs/ST_MakeLine.html


You are basically reprojecting the reference map to make it so that the straightest line between points follows more closely a curved globe (actually in this case I think it is conic, but minor detail). Next, you segmentize so that you have waypoints all along that line that follow the shortest path across the curved map. Then, your reproject back to WGS84. If you just grabbed the start and end points in the curved world, when you reproject back it would still just be a straight line. By segmentizing you get the shortest line plus waypoints along that line that use the curved map. When you reproject back to WGS84 then, it will still use the same way points which will now appear curved.


Hope that makes some sense.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...