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