Wednesday, 12 September 2018

carto - Where to start with great circle maps in CartoDB?


I've pulled together a series of maps in R based on Nathan Yau's excellent Great Circles tutorial and a nice zipcode package The maps show how many people from many zip codes descended on a single zip code. Here's one example:


sample map


I'd like to reproduce the maps in CartoDB (which uses PostGIS, I believe), but I'm not sure where to start.



My read of what I was doing in R was generating a matrix called "inter" that was full of lat/long points, drawing the line connecting those points and then starting over. So one solution that occurs to me is to have R store each line as a map line and take that series of lines to CartoDB. But I'm wondering if there's a reasonable way to just generate the maps in CartoDB directly, given an end point and a series of starting points.



Answer



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. Note that you have to also insert the CRS into the spatial_ref_sys table as explained in the same post.


enter image description here


http://cdb.io/1dEXBh9


Update: The actual query is:


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
);

and you'll need to have srid 953027 defined, so try:


SELECT * FROM spatial_ref_sys WHERE srid = 953027


And if it isn't set, set it with:


INSERT into spatial_ref_sys (srid, auth_name, auth_srid, proj4text, srtext) values ( 953027, 'esri', 53027, '+proj=eqdc +lat_0=0 +lon_0=0 +lat_1=60 +lat_2=60 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs ', 'PROJCS["Sphere_Equidistant_Conic",GEOGCS["GCS_Sphere",DATUM["Not_specified_based_on_Authalic_Sphere",SPHEROID["Sphere",6371000,0]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Equidistant_Conic"],PARAMETER["False_Easting",0],PARAMETER["False_Northing",0],PARAMETER["Central_Meridian",0],PARAMETER["Standard_Parallel_1",60],PARAMETER["Standard_Parallel_2",60],PARAMETER["Latitude_Of_Origin",0],UNIT["Meter",1],AUTHORITY["EPSG","53027"]]');

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