A global routing network
Appendum to ``Imperial roots of Global Trade’’, Journal of Economic Growth
Wessel N Vermeulen
January, 2020
1 Introduction
In this document, I intend to explain the steps that I went through to create a global network of ‘walking routes’, that is PGrouting ready network to calculate the time and distance it would take to ‘walk’ or ‘sail’ between two points on earth.
A PDF version of this document is provided here, https://www.wnvermeulen.com/empires/routing_tutorial.pdf
This projected is part of a paper with Pierre-Louis Vézina and Gunes Gokmen on the Imperial Routes of Global trade, published in Journal of Economic Growth.
I’m an autodidact in GIS methods, and while I had some basic knowledge in SQL when I started, I learned a lot on how to do things efficiently within this project. Nevertheless, I’m sure there are still a lot of inefficient parts in my code. When coding, everyone makes a continuous trade-off between (a) quick and dirty coding but potentially slow execution time and (b) slow and elegant coding with fast execution. Most of the time, I go for fast solutions only when the execution is unbearably long, or when I have epiphanies. With every block of code I learned a bit more, and am better able to implement efficient solutions early on, such that the difference between (a) and (b) becomes smaller over time.
Many solutions were adapted from online sources, in particular stackexchange. However, I did not keep track of all these sources, so references for the implementations that I used are very scarce. I’m sorry for that, since the online community has been invaluable and people deserve a huge credit for offering their solutions publicly. Perhaps my full solution posted here can be a form of payback. My implementation is also quite similar to http://orbis.stanford.edu and probably I used pieces of code, or strategies from their implementation.
I did this in PostGIS rather than the seemingly more commonly used Python with ArcGIS or QGIS. I needed routing, and at the time that I started with this, the most readily available solution seemed to be PGRouting. I use QGIS to create maps and to check that spatial manipulation in PostGIS is doing what I intent, but not to actually do any ‘real’ analysis or manipulation. The only exception, as will be noted below, is if sometimes I have to add or remove individual edges based on their specific location on the globe.
1.1 Game plan
The basic sequence is as follows
- create database
- load relevant external data, e.g. global altitude data
- create a globe with dots, clip by land.
- connect dots to a mesh, attach altitude information for land-based dots
- calculate distance and time-cost for walking, given altitude difference
- repeat the above for water, with lower resolution and different costs of traversing water based on winds.
- calculate optimal routes for current day cities, to each nearest city.
- network of routes becomes global routing network.
There is a lot of housekeeping for each, and in between, steps. Removal and creation of indices, removal and recreation of (temporary) tables, etc. As with other types of coding and data projects, over 90% of the code is housekeeping.
One issue with Postgress, is that it is not straight forward to use multiple treads to do time-consuming operations. For instance, when you have millions of rows, and a process on each row takes a few milliseconds, you’ll be waiting a long-time. However, when the operation is completely independent across rows, such an operation could be easily split between processes. I do this manually, by starting new threads in new consoles. Therefore, the script below should probably not be run in one go. Each block needs to be separately called in and waited for completion before starting the next. I looked into ways to automate multi-threading, but found it to much of a hassle to get it going so left it out.
2 Prerequisites
You’ll need an installed Postgres system, with Postgis and PGrouting installed. I followed this person’s instructions: http://www.kyngchaos.com/software/postgres/
You need a large and fast harddisk to do this. I coded everything on a 2014 iMac with 16GB of RAM, but had to buy a 512GB ssd drive to put the database on. SSD because otherwise things take way too long, and additional drive because especially the global raster on altitudes is space intensive. When I interact with the system, e.g. when loading data into the database, I only put code that I used for my system. A windows version would probably be similar.
My database runs on my local machine, not online. I use a standard set of username/password, and the database is called ‘empires’ after the main project for which this was done. Any paths to the local file system to load or store data are indicated as <path>, and need to be replaced with something suitable.
This tutorial is not an introduction to PostgreSQL and PostGIS. So, to really understand what’s going on, some prior knowledge is useful. Perhaps some minimum required knowledge is the following
- SQL commands ends with
;. So when copying, this is how to recognize how commands break over lines. - when copying from the PDF version of this tutorial, things might go wrong quickly. While I use indentation for clarity of code, these indents should be spaces, not tabs. At least that’s my experience.
- The PostgreSQL are not case sensitive. There is an old convention to have all commands and terms in capitals, e.g.
SELECT * FROM cities WHERE capital = TRUE, but I’m too lazy to implement this structurally and in any decent editor with syntax highlighting the colour schemes replaced the issue of making code more readable. Yet the capitalisation remained, albeit not perfectly. - I cut lines to make things more readable. The character for continuation on the next line is
/in python and\in bash.
3 Setting up the database and loading data
Create database with required features, then log in using the postgres account,in Bash:
createdb -U wessel empires
psql -U postgresthen in the postgres console:
GRANT ALL PRIVILEGES ON DATABASE empires to wessel;Get out of psql console (cntl+d) and then log back in using your own user account, in bash:
psql -U wessel -d empires(in Postgres) Attach the required extensions to the empires database.
CREATE EXTENSION postgis;
CREATE EXTENSION postgis_topology;
CREATE EXTENSION pgrouting;Given the global scale of the project, we cannot use a projection without significantly distorting distance calculations later. Therefore, we use the standard EPSG:4326 projection, which saves all data as long/lat coordinates, and apply functions where necessary that account for the spherical dimension. Take-away: load or save all external data in 4326.
3.1 Global altitude data
Global altitude data was taken from http://www.viewfinderpanoramas.org/Coverage%20map%20viewfinderpanoramas_org3.htm. I used two python files to download all the tiles and load them as a raster table in the database. This is a lot of data, around 16GB in zip files, and in the database the table is (with indices) 25GB.
Downloading
from lxml import html
import requests
import urllib2
from os import listdir
from sys import stdout
url = 'http://www.viewfinderpanoramas.org/' /
+ 'Coverage%20map%20viewfinderpanoramas_org3.htm'
page = requests.get(url)
tree = html.fromstring(page.text)
links = tree.xpath('/html/body/map/area/@href')
# build in check for file existence
from os import listdir
from os.path import isfile, join
mypath = '<path>/dem'
files_in_dir = [ f for f in listdir(mypath) if isfile(join(mypath,f)) ]
i=0
for link in links:
filename = link.split('/')[-1]
if filename not in files_in_dir:
if filename not in ['01-15.zip', '16-30.zip',
'31-45.zip', '46-60.zip']:
stdout.write(filename+'\n')
stdout.flush()
file = urllib2.urlopen(link)
with open(join(mypath,filename), 'wb') as output:
while True:
data = file.read()
if data:
output.write(data)
else:
breaktransferring to the database
from sys import stdout
import zipfile
from os import listdir, remove
from os.path import isfile, join, basename
from subprocess import call
import shutil
def unzip(source_filename, dest_dir):
with zipfile.ZipFile(source_filename) as zf:
for member in zf.infolist()[1::]:
file = member.filename.split('/')[1::]
filename = basename(member.filename)
#print filename
if not filename: continue
zf.extract(member, dest_dir)
shutil.move(join(dest_dir,member.filename),
join(dest_dir,filename))
mypath = '<path>/dem/'
onlyfiles = [ f for f in listdir(mypath) if isfile(join(mypath,f)) ]
i=0
for file in listdir(mypath):
if file.endswith(".zip"):
stdout.write(str(i)+" "+file+'\n')
stdout.flush()
os.makedirs(join('<path>'+'temp'))
unzip(join(mypath,file), join('/Volumes/WD/','temp'))
if i==0: # prepare table on first run
command = "raster2pgsql -s 4326 -p -M -t 100x100 -N -32768" /
+ "<path>/temp/*.hgt dem3 | psql -U wessel -d empires"
else: # insert data
command = "raster2pgsql -s 4326 -a -t 100x100 -N -32768" /
+ "<path>/temp/*.hgt dem3 | psql -U wessel -d empires"
return_code = call(command, shell=True)
return_code = call("rm -R <path>/temp/" , shell=True)
i+=1
# vacuum and index afterwardsYou could probably do it in one go, but for my development I chose to do it in steps. Keeping the raw altitude data can also become useful when creating tables for location specific maps, although in principle it could also be extracted from the table that is created here.
Housekeeping
CREATE INDEX dem3_gist_idx ON dem3 USING GIST (ST_ConvexHull(rast));3.2 External shapefiles
These shapefiles are obtained from https://www.naturalearthdata.com/downloads/
From bash, load land mass shapefile, current day countries, lakes and rivers. (The \ signs are just for clarity and to allow for lines to break. Each command starts with shp2pgsql.)
shp2pgsql -I -s 4326 -W "LATIN1" <path>/ne_50m_land/ne_50m_land.shp \
land | psql -d empires -U wessel
shp2pgsql -I -s 4326 -W "LATIN1" <path>/ne_50m_admin_0_countries.shp \
countries | psql -d empires -U wessel
shp2pgsql -d -s 4326 -I -W "LATIN1" <path>/ne_10m_lakes.shp \
public.lakes | psql -d empires -U wessel
shp2pgsql -c -s 4326 -I -W "LATIN1" \
<project path>/ne_10m_rivers_lake_centerlines.shp \
public.rivers | psql -d empires -U wessel4 The Grid: a globe with dots
Imagine the earth. Now we would like to have dots at equally spaced distances from each other, we’ll call this the grid. This only works when keep either latitude or longitude constant. Let’s keep latitude (the vertical dimension) constant, and choose the equator (lat=0) to start. Then, given the circumspection of the earth (360 degrees, and around 44000 km), and a desired distance between two dots (e.g. 5 km), you can create a string (or ring) of dots at equally spaced distances.
Moving to the next latitude, about 5 km north, we’d need to calculate the number of required dots given that the distance required to go the full 360 degrees will have slightly reduced. I work this backwards. First I create an index table that holds the information on the latitudes that the rings of dots need to be created and I calculate the distance in longitudal degrees for each latitude. Then I have all the information to create the actual grid.
drop table if exists grid_index cascade;
create table grid_index (
id serial primary key,
latt real,
latt_d real, -- latitude geography km distance from 0 degrees
latt_d_57 real, -- latitude geography km distance from 57 degrees
long_d real -- how many degrees is 1 (or other defined) km
);
create index grid_index_latt_idx on grid_index(latt);
create index grid_index_latt_d_idx on grid_index(latt_d);
create index grid_index_long_d_idx on grid_index(long_d);From point 0,0 on earth, looking towards the North pole, gives me the distance points that are equally spaced at 5000 meters (5km). A more precise routing mesh would require here to set the distance to something smaller, but this will also increase the disk space and computer power required to do everything below. I took 5km as reasonable resolution that still allows for sufficient detail of global altitude, yet manageable with my home system. I don’t need to go the full 90 degrees, since there is no land that high. So I stop at 81 degrees north.
insert into grid_index (latt_d) (
select y
from
generate_series(0,
floor(
st_distance(
geography(st_point(0,0)), geography(st_point(0,81)))
)::int, 5000) as y
);Do the same, going south. Now I can stop already at -57 degrees, because we do not care about Antarctica.
insert into grid_index (latt_d) (
select -1*y
from
generate_series(1000,
floor(
st_distance(
geography(st_point(0,0)), geography(st_point(0,-57)))
)::int, 5000) as y
);Then given the distances, I need to get the degrees of these points. I do this by creating a buffer from the origin point with a distance as large as defined in the above. Given the buffer, obtain the maximum altitude point, which is the number I want. However, buffers don’t work nice after some extreme distance. So I take new reference point after 57 degrees
update grid_index set
latt = st_ymax(st_buffer(st_point(0,0)::geography,latt_d)::geometry)
where latt_d>0
and
latt_d<=st_distance(st_point(0,0)::geography,st_point(0,57)::geography);
update grid_index set
latt = st_ymin(st_buffer(st_point(0,0)::geography,-1*latt_d)::geometry)
where latt_d<0;
-- update for equator only
update grid_index set latt = 0 where latt_d = 0;
update grid_index set latt_d_57 = floor(latt_d-st_distance(
geography(st_point(0,0)), geography(st_point(0,57))))
where
latt_d>st_distance(st_point(0,0)::geography,st_point(0,57)::geography);
update grid_index set
latt =
st_ymax(st_buffer(st_point(0,57)::geography,latt_d_57)::geometry)::real
where
latt_d>st_distance(st_point(0,0)::geography,st_point(0,57)::geography);
update grid_index set
long_d =
st_xmax(st_buffer(st_point(0,latt)::geography,5000)::geometry)::real;The final preparatory step is to realise that we only need coordinates where there is land. The obvious solution is to go ahead with a full set of coordinates and then clip this against geometries of land. However, this turned out to be a very time-consuming and inefficient procedure, because of the huge number of dots in play. So instead, for each latitude that we have defined above, I create a line across the globe and intersect this with the land geometry, so I have vertical lines that only exist where there is land, and exactly overlay the rings of dots that will be created later. Again this is a sort of index table, where each row holds a combination of a latitude and longitudal limits that define the extent (through min and max) of the grid to be created.
-- define min and max long_d for each latt in land for each polygon.
drop table if exists long_limits;
create table long_limits as(
SELECT
latt::real,
st_xmin(clipped_geom)::real long_min ,
st_xmax(clipped_geom)::real long_max
FROM (
SELECT
latt,
gid,
(ST_Dump(ST_Intersection(box.geom, land.geom))).geom As clipped_geom
FROM (
select
latt ,
st_setsrid(
st_makeline(st_point(-180, latt),St_point(180, latt)),
4326) as geom
FROM grid_index ) as box
INNER JOIN land
ON ST_Intersects(box.geom, land.geom)
) As clipped
);
create index long_limits_latt on long_limits(latt);Now I’m ready to create rows of x coordinates for each latitude, the grid. This is the first time that I split the code over four processes that I execute separately by opening four windows of the postgres console, based on latitudes. Each process creates its own table, after which the four tables are joined in a view. Since the grid is not my final objective, the objective is a mesh that connects each grid-point, this was my solution to deal with scale and speed issues.
First, prepare the four tables.
drop table if exists grid_1;
create table grid_1 (
latt real,
long real
);
drop table if exists grid_2;
create table grid_2 (
latt real,
long real
);
drop table if exists grid_3;
create table grid_3 (
latt real,
long real
);
drop table if exists grid_4;
create table grid_4 (
latt real,
long real
);
Then, the code for each table is identical, except for the table reference and the limits of the coordinates that it deals with. This will create 134,495,648 points. In order to make the id’s unique between the tables, set the start of the sequences of the serials manually to 0; 200,000,000; 300,000,000; 400,000,000.
I first create a full ring of dots, after which the relevant ones are filtered out using the long_limits that were defined above.
Table 1
insert into grid_1 (latt, long) (
with idx as (
select *
from grid_index
where
grid_index.latt<(-57 + (57+81)/4) --the latitude limit for Table 1
)
select
idx.latt,
generate_series((-180*10^5)::int, (180*10^5)::int,
floor(idx.long_d*10^5)::int)::real/(10^5) as long
from idx); -- this creates the longitudal coordiantes, by going from
-- -180 degrees to 180 degrees in steps of long_d. The
-- multiplication by 10^5 is to avoid computer rouding
-- along the way
-- filter the points above to those that are relevant,
-- these therefore filter out all coordinates that are in oceans
drop table if exists grid_1_filtered;
create table grid_1_filtered as (
select
grid_1.latt, grid_1.long
from
grid_1
join
long_limits
on (grid_1.latt=long_limits.latt and
grid_1.long >= long_limits.long_min and
grid_1.long<=long_limits.long_max
)
);
--here we create the actual geospatial grid points.
alter table grid_1_filtered add column geom GEOMETRY;
update grid_1_filtered set
geom = st_setsrid(ST_POINT(long::real, latt::real),4326);
-- housekeeping
drop table grid_1;
alter table grid_1_filtered rename to grid_1;
create index grid_1_geom on grid_1 using gist(geom);
-- note the adjustment of the pkey sequence,
-- to allow for unique pkey when the four tables are combined.
create sequence grid_1_pkey_seq start with 1;
alter table grid_1
add column pkey INTEGER DEFAULT nextval('grid_1_pkey_seq');
alter sequence grid_1_pkey_seq owned by grid_1.pkey;
create index grid_1_pkey_idx on grid_1(pkey);
alter table grid_1 add constraint grid_1_pkey_chk CHECK(pkey<200000000);
-- adding altitude by looking up the value in dem3 raster table.
alter table grid_1 add column alt real;
update grid_1 set alt = st_value(rast, geom)
from dem3 where st_intersects(rast, geom);
This creates the required points starting at 57 degrees south.
A cut out from around Tasmania
Table 2
insert into grid_2 (latt, long) (
with idx as (
select * from grid_index
where
grid_index.latt>=(-57 + (57+81)/4) AND
grid_index.latt<(-57 + 2*(57+81)/4)
)
select idx.latt, generate_series((-180*10^5)::int, (180*10^5)::int,
floor(idx.long_d*10^5)::int)::real/(10^5) as long from idx);
drop table if exists grid_2_filtered;
create table grid_2_filtered as (
select grid_2.latt, grid_2.long
from
grid_2
join
long_limits on
(grid_2.latt=long_limits.latt and
grid_2.long >= long_limits.long_min and
grid_2.long<=long_limits.long_max
)
);
alter table grid_2_filtered add column geom GEOMETRY;
update grid_2_filtered
set geom = st_setsrid(ST_POINT(long::real, latt::real),4326);
drop table grid_2;
alter table grid_2_filtered rename to grid_2;
create index grid_2_geom on grid_2 using gist(geom);
create sequence grid_2_pkey_seq start with 200000001;
alter table grid_2
add column pkey INTEGER DEFAULT nextval('grid_2_pkey_seq');
alter sequence grid_2_pkey_seq owned by grid_2.pkey;
create index grid_2_pkey_idx on grid_2(pkey);
alter table grid_2
add constraint grid_2_pkey_chk CHECK(pkey>200000000 and pkey<300000000);
alter table grid_2
add column alt real;
update grid_2
set alt = st_value(rast, geom) from dem3 where st_intersects(rast, geom);
Table 3
insert into grid_3 (latt, long) (
with idx as (
select * from grid_index
where grid_index.latt>=(-57 + 2*(57+81)/4) AND
grid_index.latt<(-57 + 3*(57+81)/4)
)
select
idx.latt,
generate_series((-180*10^5)::int, (180*10^5)::int,
floor(idx.long_d*10^5)::int)::real/(10^5)
as long from idx
);
drop table if exists grid_3_filtered;
create table grid_3_filtered as (
select grid_3.latt, grid_3.long from
grid_3
join
long_limits
on (grid_3.latt=long_limits.latt and
grid_3.long >= long_limits.long_min and
grid_3.long<=long_limits.long_max)
);
alter table grid_3_filtered add column geom GEOMETRY;
update grid_3_filtered
set geom = st_setsrid(ST_POINT(long::real, latt::real),4326);
drop table grid_3;
alter table grid_3_filtered rename to grid_3;
create index grid_3_geom on grid_3 using gist(geom);
create sequence grid_3_pkey_seq start with 300000001;
alter table grid_3
add column pkey INTEGER DEFAULT nextval('grid_3_pkey_seq');
alter sequence grid_3_pkey_seq owned by grid_3.pkey;
create index grid_3_pkey_idx on grid_3(pkey);
alter table grid_3
add constraint grid_3_pkey_chk CHECK(pkey>300000000 and pkey<400000000);
alter table grid_3 add column alt real;
update grid_3
set alt = st_value(rast, geom) from dem3 where st_intersects(rast, geom);
and finally, Table 4
insert into grid_4 (latt, long) (
with idx as (
select * from grid_index where grid_index.latt>=(-57 + 3*(57+81)/4)
)
select idx.latt,
generate_series(
(-180*10^5)::int, (180*10^5)::int,floor(idx.long_d*10^5)::int
)::real/(10^5) as long from idx);
drop table if exists grid_4_filtered;
create table grid_4_filtered as (
select grid_4.latt, grid_4.long
from grid_4
join long_limits
on (grid_4.latt=long_limits.latt and
grid_4.long >= long_limits.long_min and
grid_4.long<=long_limits.long_max
)
);
alter table grid_4_filtered add column geom GEOMETRY;
update grid_4_filtered
set geom = st_setsrid(ST_POINT(long::real, latt::real),4326);
drop table grid_4;
alter table grid_4_filtered rename to grid_4;
create index grid_4_geom on grid_4 using gist(geom);
create sequence grid_4_pkey_seq start with 400000001;
alter table grid_4
add column pkey INTEGER DEFAULT nextval('grid_4_pkey_seq');
alter sequence grid_4_pkey_seq owned by grid_4.pkey;
create index grid_4_pkey_idx on grid_4(pkey);
alter table grid_4 add constraint grid_4_pkey_chk CHECK(pkey>400000000);
alter table grid_4 add column alt real;
update grid_4
set alt = st_value(rast, geom)
from dem3 where st_intersects(rast, geom);We can check that this is correct, by taking two sequential points and calculate their distance.
-- distance in degrees between two points:
select st_distance(p1.geom, p2.geom) dist
from grid_4 as p1, grid_4 as p2
where p1.pkey=401316438 and p2.pkey=401316439;
-- in meters
select st_distance(geography(p1.geom), geography(p2.geom)) dist
from grid_4 as p1, grid_4 as p2
where p1.pkey=401316438 and p2.pkey=401316439;Combine the four tables in a view.
create view grid as (
select * from grid_1
union
select * from grid_2
union
select * from grid_3
union
select * from grid_4);5 Connect the dots to create a mesh
The idea is to connect each dot with its nearest neighbour. This is fine when keeping latitude constant, but when going up or down it can be tricky. The grid is not rectangular, since we’re dealing with a globe.
First, prepare the mesh table. It’s called o_mesh, since I borrowed the code from somewhere in stackexchange.
drop table if exists o_mesh CASCADE;
create table o_mesh
(
-- id BIGSERIAL PRIMARY KEY,
"source" int,
target int,
direction smallint, -- varchar(2),
geom GEOMETRY('LINESTRING', 4326),
dist real
);The following is repeated for each of the four grid tables, by just replacing grid_1 with grid_2, grid_3 and grid_4 sequentially.
The idea is to connect each two points together once. This sounds quite straightforward, but what we have is a large set of dots, with no information on which two dots are closest. Moreover, the dots are only equally spaced longitudinal given a latitude, at 5km. We also know that these rings are 5km spaced from each other in latitude. But these do not line-up exactly. So we know that the maximum distance is sqrt(52+52), but this is in km, where we need it in degrees, which will vary as we move up the globe (i.e. at different latitudes). Therefore, the maximum allowed distance in degrees must be calculated for each latitude separately. (I guess people that understand spherical calculations better than I do, will have a closed form formula for this. I do not.)
We also only want each link once, so we need to limit links to only a quarter of the quadrant, pure east, pure north, and anything north-east. This is implemented at the very end.
INSERT INTO o_mesh (source, target, direction, geom, dist)
with indeks as (
select
idx1.*,
idx2.latt+0.000001 as next_latt,
idx3.latt-0.000001 as prev_latt,
st_distance(st_point(0,idx1.latt), st_point(idx1.long_d,idx2.latt))
from
(
select row_number() over () id1, *
from grid_index order by latt
) idx1
left join
(
select (row_number() over () -1) id2, *
from grid_index order by latt
) idx2
on id1=id2
left join
(
select (row_number() over () +1) id3, *
from grid_index order by latt
) idx3
on id1=id3
)
SELECT
source.pkey as source, target.pkey as target, (
case
when source.long = target.long AND source.latt < target.latt then 1 --'N'
when source.long < target.long AND source.latt < target.latt then 2 --'NE'
when source.long < target.long AND source.latt = target.latt then 3 --'E'
when source.long < target.long AND source.latt > target.latt then 4 --'SE'
when source.long = target.long AND source.latt > target.latt then 5 --'S'
when source.long > target.long AND source.latt > target.latt then 6 --'SW'
when source.long > target.long AND source.latt = target.latt then 7 --'W'
when source.long > target.long AND source.latt < target.latt then 8 --'NW'
else 9 --'ER'
end ) as direction,
-- the line that represent the connection
st_setsrid(
ST_MakeLine(
source.geom::GEOMETRY, target.geom::GEOMETRY
), 4326) AS geom,
-- the distance in meters between two points
st_distance(geography(source.geom), geography(target.geom)) AS dist
FROM
-- distance should match what is made in the grid.
(select
grid_1.*, dist_limit, next_latt, prev_latt
FROM
grid_1
left join
indeks on grid_1.latt=indeks.latt order by grid_1.latt
) AS source
JOIN
(select *, FROM grid_1) AS target
ON
-- one way links only
(
target.long > source.long
AND target.latt between (source.prev_latt) and (source.next_latt)
OR
(target.long = source.long
AND target.latt > source.latt AND target.latt <= next_latt)
)
and ST_DWithin(target.geom, source.geom, dist_limit);
A cut out from around Tasmania
The preceding, when done for each of the 4 grid tables, neglects that we need links between them as well. I fix this separately. The reason that I did not do the previous using the view has probably something to do with the fact that the indexes on the view do not work well, and having one large table with all grid points was also unwieldy.
Note that while I create an actual geospatial line that I can plot in a map, this is not required for any of the calculation. The table is nothing more than an index that indicates which two nodes (grid points) are connected, with added information on the ‘cost’ to traverse the length. This cost can be what ever we want it to be, geographical distance, time, or anything else.
-- insert where connections overlap partial grid tables.
create index grid_1_latt_idx on grid_1(latt);
create index grid_2_latt_idx on grid_2(latt);
create index grid_3_latt_idx on grid_3(latt);
create index grid_4_latt_idx on grid_4(latt);
-- fixing missing links. this is slower than it should be. there are just
-- 64772 grid points, but the 'with'-query probably looses the index.
-- maybe a temporary table with index on geom would be better.
drop table if exists border_grid;
create table border_grid as (
select * from grid_1 where latt= (select max(latt) from grid_1)
union
select * from grid_2
where
latt= (select min(latt) from grid_2 )
OR
latt= (select max(latt) from grid_2)
union
select * from grid_3
where
latt= (select min(latt) from grid_3 )
OR
latt= (select max(latt) from grid_3)
union
select * from grid_4
where latt= (select min(latt) from grid_4 )
);
create index border_grid_geom_idx on border_grid using gist(geom);
with indeks as (
select
idx1.*,
idx2.latt+0.000001 as next_latt,
idx3.latt-0.000001 as prev_latt,
st_distance(
st_point(0,idx1.latt), st_point(idx1.long_d,idx2.latt)
) dist_limit
from
(select row_number() over () id1, *
from grid_index order by latt) as idx1
left join
(select (row_number() over () -1) id2, *
from grid_index order by latt) as idx2
on id1=id2
left join
(select (row_number() over ()+1) id3, *
from grid_index order by latt) as idx3
on id1=id3
)
INSERT INTO o_mesh (source, target, direction, geom, dist)
SELECT
source.pkey as source, target.pkey as target, (
case
when source.long = target.long AND source.latt < target.latt then 1 --'N'
when source.long < target.long AND source.latt < target.latt then 2 --'NE'
when source.long < target.long AND source.latt = target.latt then 3 --'E'
when source.long < target.long AND source.latt > target.latt then 4 --'SE'
when source.long = target.long AND source.latt > target.latt then 5 --'S'
when source.long > target.long AND source.latt > target.latt then 6 --'SW'
when source.long > target.long AND source.latt = target.latt then 7 --'W'
when source.long > target.long AND source.latt < target.latt then 8 --'NW'
else 9 --'ER'
end ) as direction,
st_setsrid(
ST_MakeLine(source.geom::GEOMETRY,
target.geom::GEOMETRY), 4326) AS geom,
st_distance(geography(source.geom), geography(target.geom)) AS dist
FROM
-- distance should match what is made in the grid.
(select border_grid.*, dist_limit, next_latt , prev_latt
FROM border_grid left join indeks
on border_grid.latt=indeks.latt order by border_grid.latt
) AS source
JOIN
(select * FROM border_grid) AS target
ON
-- one way links only
(
target.long > source.long
AND target.latt between (source.prev_latt) and (source.next_latt)
OR
(target.long = source.long
AND target.latt > source.latt AND target.latt <= next_latt)
)
and ST_DWithin(target.geom, source.geom, dist_limit);
drop table border_grid;Housekeeping.
-- vacuum after entering so much data
VACUUM ANALYZE;
-- creating indices: do this after insertions.
CREATE INDEX o_mesh_source_idx ON o_mesh("source");
CREATE INDEX o_mesh_target_idx ON o_mesh("target");
-- CREATE INDEX o_mesh_geog_idx ON o_mesh using GIST(geog);
CREATE INDEX o_mesh_geom_idx ON o_mesh using GIST(geom);
alter table o_mesh add column pkey serial primary key;
-- vacuum after entering so much data
VACUUM ANALYZE;6 Walking cost
We already have the distance between each two points, around 5km. Next is to calculate the walking time it requires to go from one point to the other. For this we take into account the altitude difference, and the direction.
First prepare the mesh table for the additional information.
alter table o_mesh add column cost_elev real;
alter table o_mesh add column source_alt real; --altitude source
alter table o_mesh add column target_alt real; -- altitude origin
alter table o_mesh add column angle real; -- angle
alter table o_mesh add column hours real; -- travel hours
alter table o_mesh add column rev_hours real; -- " reverse direction
alter table o_mesh add column "type" smallint; -- land, water, cross-over
Add information on the altitude.
update o_mesh set source_alt = alt
from grid where o_mesh.source=grid.pkey;
update o_mesh set target_alt = alt
from grid where o_mesh.target=grid.pkey;I use Naismith’s rule https://en.wikipedia.org/wiki/Naismith%27s_rule with additions, which states the following - Allow 1 hour for every 5 kilometres (3.1 mi) forward, plus 1 hour for every 600 metres (2,000 ft) of ascent. - Aitken - Langmuir corrections – Subtract 10 minutes for every 300 meters of descent for slopes between 5 degrees and 12 degrees – add 10 minutes for every 300 meters of descent for slopes greater than 12 degrees.
Actually, this is a proper function: \[\text{velocity} = 6 * \exp ( -3.5 * | \tan (angle)+0.05|)\]
I add to this that the angle cannot be more than 12 degrees. If the slope is higher than that, a walker will lengthen its path to reach destination; making z-shapes to ascend a mountain.
So, the length of slope = sqrt(geodesic length^2 + height^2) =height/sin(angle) If angle>0.12, geodesic distance: height/tan(0.12), slope length: height/sin(0.12)
So, then cost = distance / speed = (height/sin(angle)) / (1000 * 6EXP(-3.5ABS(TAN(angle)+0.05)))
To implement this I distinguish 3 cases, no angle (flat surface), inclined surfaces up to 12 degrees, and inclination over 12 degrees. The the whole thing needs repetition for when the surface slopes down.
-- note transformation from radients to degrees and back.
update o_mesh
set angle = atan2((target_alt-source_alt),dist)
where pkey>=0E6 and pkey<4E6; -- determines range for processes
update o_mesh set hours = (
case
when angle=0. THEN dist/ 5036.742
when abs(angle) <= 12.0/180.0*pi()
THEN (target_alt-source_alt)/sin(angle) /
(1000.0 * 6.0*EXP(-3.5*ABS(TAN(angle)+0.05)))
else
(abs(target_alt-source_alt)/sin(12.0/180.0*pi())) /
(1000.0 * 6.0*EXP(-3.5*ABS(TAN(sign(angle)*12.0/180.0*pi())+0.05)))
end
),
rev_hours = ( -- only difference is the minus in the dominator
-- adjusting speed, distance is the same up or down.
case
when angle = 0 THEN dist/ 5036.742
when abs(angle) <= 12.0/180.0*pi()
THEN (target_alt-source_alt)/sin(angle) /
(1000.0 * 6.0*EXP(-3.5*ABS(TAN(-1*angle)+0.05)))
-- alternative: (dist/1000)/COS(angle)
else (abs(target_alt-source_alt)/sin(12.0/180.0*pi())) /
(1000.0 * 6.0*EXP(-3.5*ABS(TAN(sign(angle)*-1*12.0/180.0*pi())+0.05)))
end
),
"type" = 1 where pkey>=0 and pkey<4E6; -- needs to be same as above
-- cleaning up the table
delete from o_mesh where
source is null or
target is null or
source_alt is null or
target_alt is null or
dist is null; Next is to add time costs of crossing water-bodies, like rivers. The arbitrary assumption is 3 hours, which is quite a lot when recognising that the normal distance between two neighbouring points on a flat surface would be crossed in around 1 hour.
First prepare a table, tw for traversing water.
create table upd_tw (pkey int);Then split the following over multiple processes by pkey to speed up the insertions.
insert into upd_tw
select pkey from o_mesh where exists(select 1 from
(select gid, geom from rivers
union all
select gid, geom from lakes
) as water
where st_intersects(o_mesh.geom, water.geom))
and pkey>=0E6 and pkey<2E6; -- determines range for processes
create index upd_tw_pkey on upd_tw(pkey);Then update the costs in one go.
begin;
set local work_mem = '5000 MB';
-- create new table with new column filled with default values
create table o_mesh_new as
select source, target, direction, geom, dist,
pkey, cost_elev, source_alt, target_alt,
angle hours, rev_hours, 1 as "type",
0 as travers_waters, hours as hours_tw,
rev_hours as rev_hours_tw
from o_mesh;
create index o_mesh_new_pkey on o_mesh_new(pkey);
-- update new columns where necessary
-- set the crossing time for lakes and rivers here as original hours + X
update o_mesh_new set
travers_waters=1,
hours_tw=hours + 3,
rev_hours_tw=rev_hours + 3
where pkey in (select pkey from upd_tw);
-- replace old table with the new
drop table o_mesh;
alter table o_mesh_new rename to o_mesh;
-- recreate primary key
alter table o_mesh
drop column pkey,
add column pkey bigserial primary key;
--alter table o_mesh add constraint
create index o_mesh_geom_idx on o_mesh using gist(geom);
create index o_mesh_source_idx on o_mesh(source);
create index o_mesh_target_idx on o_mesh(target);
end;Housekeeping
vacuum full analyze o_mesh;Finally, it is quite straightforward to extent the cost travelling between any two points with additional information. Say, with information from a raster on specific geology zones (jungle, desert, savannah, steppe, etc.), one could add this information by associating such data to the mesh, and adding some multiplication constant to the hours for each type of zone.
7 Water
Adding a water-based network is to a large extent identical to the land based network. There is one major simplification: dots are not spaced 5km apart but 50km. That reduces the size of the network substantially, and therefore allows it to be created quite a bit faster as well. A complication is that the date line (where -180 meets +180) needs to be dealt with structurally. This was sidestepped in the land mesh, but will still have to solved later.
The relevant grid_index_water is just a subset from grid_index, because we do not use every 5km, but every 50km. the mod function is useful for this.
drop view if exists grid_index_water;
create view grid_index_water as (
select id, latt, latt_d latt_d_57,
st_xmax(st_buffer(st_point(0,latt)::geography,50000)::geometry)::real
as long_d
from grid_index
where
mod(latt_d::int,50000::int)=0 or
mod(-latt_d::int,50000::int)=1000
);The straightforward solution would be to create all the points. Cut out those that overlap land. In this case this is feasible, because the number of points that we deal with is much smaller. An alternative solution would have followed the one for land.
The water grid:
drop table if exists grid_water;
create table grid_water (
latt real,
long real
);
insert into grid_water (latt, long) (
with idx as (
select * from grid_index_water
where grid_index_water.latt>-57 and grid_index_water.latt<81
)
select idx.latt, generate_series((-180*10^5)::int, (180*10^5)::int,
floor(idx.long_d*10^5)::int)::real/(10^5) as long
from idx);
South and Western Europe
Adding manually some points in Sea of Marmara, to facilitate navigation from Mediterranean to the Black Sea. This issue comes up because of the low resolution, but I don’t want to increase the resolution for the globe for this small water body.
insert into grid_water(latt, long)
VALUES (40.889, 27.624), (40.889, 28.319), (40.889, 28.984);
Then create all points, and delete those that cross land.
alter table grid_water add column geom GEOMETRY('POINT', 4326);
update grid_water
set geom = st_setsrid(ST_POINT(long::real, latt::real),4326);
create index grid_water_geom on grid_water using gist(geom);
delete from grid_water using land
where st_intersects(grid_water.geom, land.geom);
-- this is slow because we're comparing each of 69k rows with 133k rows.
alter table grid_water add column pkey serial primary key;
South and Western Europe
Additionally, I would like to know which dots are within 200 km of a coastline. This 200 km is used in the analysis to make a distinction between crossing small water bodies, such as the Mediterranean or connect islands, but prohibit ocean crossings. These buffers will interfere with the date-line, so we need to deal with that separately.
Starting from the longitudinal limits we created for land, take 200km buffers around these, with spatial union where necessary and use this to identify all water dots within 200km of a coast line.
drop table if exists land_buffers;
create table land_buffers as (
with points_1 as (
select latt, long_min, long_max,
(st_dump(
st_split(
st_buffer(
st_point(long_min, latt)::geography, 200000
)::geometry('POLYGON', 4326),
st_setsrid(st_makeline(st_point(-174, -57),
st_point(-174, 81)), 4326)
)
)).geom p1_geom,
(st_dump(
st_split(
st_buffer(
st_point(long_max, latt)::geography, 200000
)::geometry('POLYGON', 4326),
st_setsrid(st_makeline(st_point(174, -57),
st_point(174, 81)), 4326)
)
)).geom p2_geom from long_limits
)
select latt, long_min long,
p1_geom::GEOMETRY('POLYGON', 4326) geom
from points_1
union
select latt, long_min long,
p2_geom::GEOMETRY('POLYGON', 4326) geom
from points_1
);alter table land_buffers add column id serial primary key;
create index land_buffers_geom on land_buffers using gist(geom);
-- identify buffers crossing dateline
alter table land_buffers add column dl int;
with points as (
select row_number() over (PARTITION BY id) r, *
from (
select id, (st_dumppoints(st_envelope(geom))).geom
from land_buffers lb
where
(lb.geom &&
st_setsrid(st_makeline(ST_MakePoint(174,90),
ST_MakePoint(174,-90)),4326)
)
or
(lb.geom &&
st_setsrid(st_makeline(ST_MakePoint(-174,90),
ST_MakePoint(-174,-90)),4326)
)
) foo),
points_dl as (select
p1.id,
st_astext(p1.geom) as geom_bl,
sign(st_x(p1.geom)) sign,
st_astext(p2.geom) as geom_tr,
sign(st_x(p2.geom)) sign
from points p1 join points p2 on (p1.id=p2.id and p1.r=1 and p2.r=3)
where sign(st_x(p1.geom))=-1*sign(st_x(p2.geom))
)
update land_buffers set dl = 1
from points_dl where points_dl.id = land_buffers.id;
update land_buffers set dl = 0 where dl is null;
South and Western Europe
Then if a point falls in the buffer, it should be denoted as such, using the indicator is_lb (is land border).
-- intersects land_buffers
alter table grid_water add column is_lb int;
update grid_water set is_lb = foo.is_lb
from
(select g.pkey,
bool_or(st_intersects(g.geom,lb.geom))::int as is_lb
from grid_water g, land_buffers lb
where g.geom && lb.geom and lb.dl=0 group by g.pkey) as foo
where grid_water.pkey = foo.pkey;
South and Western Europe
While this works well generally, around the +/-180 degrees longitude this is a bit of mess. So, we use a the geography method there, and work specifically on all points in the 174 to 180 and -174 to -180 degrees range.
-- grid points to be added back
drop table if exists grid_water_dl;
create table grid_water_dl as (
select * from grid_water where (long<-174.0 or long > 174.0)
);
create index grid_water_dl_geom_idx on grid_water_dl using gist(geom);Check distance of each point from land for within 200km?
alter table grid_water_dl add column is_lb int;
update grid_water_dl set is_lb = chk
from (
select g.pkey, bool_or(st_dwithin(geography(g.geom),
geography(land.geom), 200000))::int chk
from grid_water_dl g, land group by pkey) as foo
where grid_water_dl.pkey = foo.pkey; Update points in grid_water with dateline fix
update grid_water
set is_lb = 1 where pkey in (
select pkey from grid_water_dl where is_lb=1
);
Alaska
Alaska
Then creating the mesh is the same as before, but I allow for more cross-connections. The reason is a combination of the larger resolution, which would create some large shipping detours, and the constraint of winds, which cannot be directed against at certain angels. So I need to allow for more flexibility.
drop table o_mesh_water;
create table o_mesh_water (
"source" int,
target int,
direction smallint, -- varchar(2),
geom GEOMETRY('LINESTRING', 4326),
dist real
);
INSERT INTO o_mesh_water (source, target, direction, geom, dist)
-- grid_index_water with next and previous latt value
with indeks as (
select
idx1.*,
idx2.latt+0.000001 as next_latt,
idx3.latt-0.000001 as prev_latt,
st_distance(st_point(0,idx1.latt),
st_point(2*idx1.long_d,idx2.latt)) dist_limit
from
(select row_number() over (order by latt) id1, *
from grid_index_water ) idx1
left join
(select (row_number() over (order by latt) -1) id2, *
from grid_index_water ) idx2
on id1=id2
left join
(select (row_number() over (order by latt) +1) id3, *
from grid_index_water ) idx3
on id1=id3
)
SELECT
source.pkey as source, target.pkey as target, (
case
when source.long = target.long AND source.latt < target.latt then 1 --'N'
when source.long < target.long AND source.latt < target.latt then 2 --'NE'
when source.long < target.long AND source.latt = target.latt then 3 --'E'
when source.long < target.long AND source.latt > target.latt then 4 --'SE'
when source.long = target.long AND source.latt > target.latt then 5 --'S'
when source.long > target.long AND source.latt > target.latt then 6 --'SW'
when source.long > target.long AND source.latt = target.latt then 7 --'W'
when source.long > target.long AND source.latt < target.latt then 8 --'NW'
else 9 --'ER'
end ) as direction,
st_setsrid(ST_MakeLine(source.geom::GEOMETRY,
target.geom::GEOMETRY), 4326) AS geom,
st_distance(geography(source.geom), geography(target.geom)) AS dist
FROM
-- distance should match what is made in the grid.
(select grid_water.*, dist_limit, next_latt , prev_latt
FROM
grid_water
left join
indeks
on grid_water.latt=indeks.latt
order by grid_water.latt
) AS source
JOIN
(select * FROM grid_water ) AS target
ON
-- one way links only
(
( target.long > source.long -- next long, up or down
AND target.latt between (source.prev_latt) and (source.next_latt)
AND target.latt != source.latt
AND ST_DWithin(target.geom, source.geom, dist_limit)
)
OR
( target.long > source.long -- next long straight
AND target.latt = source.latt
AND ST_DWithin(target.geom, source.geom, dist_limit/2)
)
OR
(target.long = source.long -- same long, up or down
AND target.latt > source.latt AND target.latt <= next_latt
AND ST_DWithin(target.geom, source.geom, dist_limit))
);
This creates a nice mesh. However, there are also some issues, as are apparent in the yellow lines on the map. Some lines cut across way more land than I’m comfortable with to ignore. However, I did not want to remove all lines that cross land, because that would remove too many lines. So I ended up just selecting those lines manually in QGIS, copying them to a separate layer, and removing them from the original. Not all is bad, however. Around Gibraltar the Atlantic is beautifully connected with the Mediterranean. But in other parts of the world additional links should be required. For instance, around Istanbul, and Denmark.
Southern Europe, North Africa
Black Sea
Denmark and Kattegat
This is fixed manually in the following. First I copied the crossing lines to a new shapefile in QGIS, my manual selection. Then these lines are deleted using a simple spatial match on the original. This shapefile that includes these is provided at http://www.wnvermeulen.com/empires/shp/o_mesh_water_rm.zip
Load that first (in bash console)
shp2pgsql -I -s 4326 -W "LATIN1" <path>/o_mesh_water_rm.shp \
o_mesh_water_rm | psql -d empires -U wesselThen remove from original table using a spatial match.
delete from o_mesh_water o using o_mesh_water_rm rm
where st_within(o.geom, rm.geom);Next, I add manually connections that I want to exist but did not appear due to the constraints set in the algorithm above.
-- additional connections to connect see of marmara with black sea
insert into o_mesh_water (source, target, direction, geom, dist)
select source, target,
3 as direction,
st_setsrid(ST_MakeLine(s.the_geom::GEOMETRY,
t.the_geom::GEOMETRY), 4326) AS geom,
st_distance(geography(s.the_geom),
geography(t.the_geom)) AS dist
from
(select id as source, the_geom from o_mesh_water_vertices_pgr o
order by o.the_geom <#> st_setsrid(
ST_POINT(27.624::real, 40.889::real), 4326) limit 1) s,
(select id as target, the_geom from o_mesh_water_vertices_pgr o
order by o.the_geom <#> st_setsrid(
ST_POINT(28.319::real, 40.889::real),4326) limit 1) t;
insert into o_mesh_water (source, target, direction, geom, dist)
select source, target,
4 as direction,
st_setsrid(ST_MakeLine(s.the_geom::GEOMETRY,
t.the_geom::GEOMETRY), 4326) AS geom,
st_distance(geography(s.the_geom),
geography(t.the_geom)) AS dist
from
(select id as source, the_geom from o_mesh_water_vertices_pgr o
order by o.the_geom <#> st_setsrid(
ST_POINT(27.624::real, 40.889::real),4326) limit 1) s,
(select id as target, the_geom from o_mesh_water_vertices_pgr o
order by o.the_geom <#> st_setsrid(
ST_POINT(27.758::real, 40.548::real),4326) limit 1) t;
insert into o_mesh_water (source, target, direction, geom, dist)
select source, target,
4 as direction,
st_setsrid(ST_MakeLine(s.the_geom::GEOMETRY,
t.the_geom::GEOMETRY), 4326) AS geom,
st_distance(geography(s.the_geom),
geography(t.the_geom)) AS dist
from
(select id as source, the_geom from o_mesh_water_vertices_pgr o
order by o.the_geom <#> st_setsrid(
ST_POINT(27.624::real, 40.889::real),4326) limit 1) s,
(select id as target, the_geom from o_mesh_water_vertices_pgr o
order by o.the_geom <#> st_setsrid(
ST_POINT(28.351::real, 40.548::real),4326) limit 1) t;
insert into o_mesh_water (source, target, direction, geom, dist)
select source, target,
3 as direction,
st_setsrid(ST_MakeLine(s.the_geom::GEOMETRY,
t.the_geom::GEOMETRY), 4326) AS geom,
st_distance(geography(s.the_geom),
geography(t.the_geom)) AS dist
from
(select id as source, the_geom from o_mesh_water_vertices_pgr o
order by o.the_geom <#> st_setsrid(
ST_POINT(28.319::real, 40.889::real),4326) limit 1) s,
(select id as target, the_geom from o_mesh_water_vertices_pgr o
order by o.the_geom <#> st_setsrid(
ST_POINT(28.984::real, 40.889::real),4326) limit 1) t;
insert into o_mesh_water (source, target, direction, geom, dist)
select source, target,
4 as direction,
st_setsrid(ST_MakeLine(s.the_geom::GEOMETRY,
t.the_geom::GEOMETRY), 4326) AS geom,
st_distance(geography(s.the_geom),
geography(t.the_geom)) AS dist
from
(select id as source, the_geom from o_mesh_water_vertices_pgr o
order by o.the_geom <#> st_setsrid(
ST_POINT(28.319::real, 40.889::real),4326) limit 1) s,
(select id as target, the_geom from o_mesh_water_vertices_pgr o
order by o.the_geom <#> st_setsrid(
ST_POINT(28.351::real, 40.548::real),4326) limit 1) t;
insert into o_mesh_water (source, target, direction, geom, dist)
select source, target,
2 as direction,
st_setsrid(ST_MakeLine(s.the_geom::GEOMETRY,
t.the_geom::GEOMETRY), 4326) AS geom,
st_distance(geography(s.the_geom),
geography(t.the_geom)) AS dist
from
(select id as source, the_geom from o_mesh_water_vertices_pgr o
order by o.the_geom <#> st_setsrid(
ST_POINT(28.984::real, 40.889::real),4326) limit 1) s,
(select id as target, the_geom from o_mesh_water_vertices_pgr o
order by o.the_geom <#> st_setsrid(
ST_POINT(29.406::real, 41.481::real),4326) limit 1) t;
insert into o_mesh_water (source, target, direction, geom, dist)
select source, target,
2 as direction,
st_setsrid(ST_MakeLine(s.the_geom::GEOMETRY,
t.the_geom::GEOMETRY), 4326) AS geom,
st_distance(geography(s.the_geom),
geography(t.the_geom)) AS dist
from
(select id as source, the_geom from o_mesh_water_vertices_pgr o
order by o.the_geom <#> st_setsrid(
ST_POINT(28.984::real, 40.889::real),4326) limit 1) s,
(select id as target, the_geom from o_mesh_water_vertices_pgr o
order by o.the_geom <#> st_setsrid(
ST_POINT(28.810::real, 41.481::real),4326) limit 1) t;
-- baltic sea with north sea
insert into o_mesh_water (source, target, direction, geom, dist)
select source, target,
2 as direction,
st_setsrid(ST_MakeLine(s.the_geom::GEOMETRY,
t.the_geom::GEOMETRY), 4326) AS geom,
st_distance(geography(s.the_geom),
geography(t.the_geom)) AS dist
from
(select id as source, the_geom from o_mesh_water_vertices_pgr o
order by o.the_geom <#> st_setsrid(
ST_POINT(10.813::real, 55.420::real),4326) limit 1) s,
(select id as target, the_geom from o_mesh_water_vertices_pgr o
order by o.the_geom <#> st_setsrid(
ST_POINT(11.119::real, 56.289::real),4326) limit 1) t;
insert into o_mesh_water (source, target, direction, geom, dist)
select source, target,
2 as direction,
st_setsrid(ST_MakeLine(s.the_geom::GEOMETRY,
t.the_geom::GEOMETRY), 4326) AS geom,
st_distance(geography(s.the_geom),
geography(t.the_geom)) AS dist
from
(select id as source, the_geom from o_mesh_water_vertices_pgr o
order by o.the_geom <#> st_setsrid(
ST_POINT(10.813::real, 55.420::real),4326) limit 1) s,
(select id as target, the_geom from o_mesh_water_vertices_pgr o
order by o.the_geom <#> st_setsrid(
ST_POINT(11.119::real, 56.289::real),4326) limit 1) t;To keep track of the different ‘modes’ of traveling later, we create a new integer column where each number stands for a mode: 1=land, 2=water, 3=land-water.
alter table o_mesh_water add column "type" smallint;Some housekeeping and other preparations
CREATE INDEX o_mesh_water_geom_idx ON o_mesh_water using GIST(geom);
alter table o_mesh_water add column pkey serial primary key;The following command is a preparation for the routing commands. It creates new tables of nodes, called o_mesh_water_vertices_pgr.
select pgr_createTopology('o_mesh_water', 0.0001, 'geom', 'pkey');Here I add the indicator on the dots within 200km of a coast.
alter table o_mesh_water alter column within200km set default FALSE;
update o_mesh_water set within200km = TRUE from
(select o.id from o_mesh_water_vertices_pgr o , grid_water gw
where o.the_geom=gw.geom and gw.is_lb=1) idx
where o_mesh_water.source=idx.id OR o_mesh_water.target=idx.id;7.1 Winds
For the cost of sailing we require two pieces of information. wind direction and speed across the oceans, and information on sailing speed given the wind and sailing direction.
The wind information originally comes from the QuickSCAT satellites: https://podaac-opendap.jpl.nasa.gov/opendap/allData/quikscat/L3/wind_1deg_1mo/ We used the file: sfcWind_QuikSCAT_L2B_v20110531_199908-200910.nc. This holds daily data from 1999-2009. Assuming that sailing predominantly occurs in certain months of the year we take the average speed and direction for these months and this is the information that we work with. The assumption is similar to Pascali (2017)1
The second piece of information also comes from Pascali (2017), but the working paper version, which documents polar diagram of typical sailing ship that gives the speed of sailing given the wind speed and the directions the wind is coming from relative to the direction the ship is taking.
Both files are attached here: <windsSFC.csv>, <polardiag_vec.csv>
drop table if exists winds;
create table winds (
sfcWind numeric, lat numeric, lon numeric, vas numeric,
uas numeric, direction numeric, speed2 numeric
);
\copy winds FROM 'windsSFC.csv' DELIMITERS ',' CSV HEADER
alter table winds
add column geom GEOMETRY('POINT', 4326),
add id SERIAL PRIMARY KEY;
update winds set geom = st_setsrid(st_point(lon, lat), 4326);
create index winds_geom on winds using gist(geom);
alter table winds add column speed numeric;
update winds set speed = speed2*3.6; -- from m/s to km/h
drop table if exists polar;
create table polar (
angle numeric, windspeed_raw numeric, sailspeed numeric
);
-- speeds as kmh = knots * 1.852001
\copy polar FROM 'polardiag_vec.csv' DELIMITERS ',' CSV HEADER
alter table polar add id SERIAL PRIMARY KEY;
alter table polar add column windspeed numeric;
insert into polar
select 180+(180-angle) angle, windspeed_raw, sailspeed
from polar where angle!=0;
update polar set windspeed = round(ceiling(windspeed_raw)/6)*6;
create index polar_angle_speed on polar(angle, windspeed);First things first. The winds data, with some rendering in QGIS already allows for some awesome figures.
North Atlantic
Global winds
The resolution of the wind is not as great as the one we used for the mesh, and the data is also using a regular sequence in degrees, such that it looks like a perfect rectangular grid in the 4236 projection, and makes that the resolution of the data varies in terms of real distances across the globe, the resolution is higher towards the poles, and lowest on the equator. Also note that some parts of the globe are not well covered, and some association with by distance will become imprecise in such places.
The first step is to associate wind points to the mesh, then for each line we can calculate the sailing speed given the wind information and the direction of travel. So we have the grid points of the mesh, and the locations of the winds data, the tasks is to find the closest wind point to each grid point. This wind point may differ depending on whether you are going from the grid point, or moving away from it. Since wind speed and direction change gradually from one data point to another, I’m no so concerned about this extrapolation that is essentially going on.
Southern Europe and North Africa
alter table o_mesh_water add column sailing_angle numeric,
add column sailing_angle_rev numeric,
add column nearest_wind_point integer,
add column nearest_wind_point_rev integer;
-- associate mesh points to wind table to get info for later
update o_mesh_water
set nearest_wind_point = foo.id,
nearest_wind_point_rev = foo.idr
from (
SELECT o.pkey,
(SELECT w1.id id
FROM winds w1
ORDER BY st_startPoint(o.geom) <#> w1.geom LIMIT 1),
(SELECT w2.id idr
FROM winds w2
ORDER BY st_endPoint(o.geom) <#> w2.geom LIMIT 1)
FROM o_mesh_water o
) as foo
where o_mesh_water.pkey = foo.pkey;
create index o_mesh_water_wp on o_mesh_water(nearest_wind_point);
create index o_mesh_water_wp_rev on o_mesh_water(nearest_wind_point_rev);Now we have the wind angle in absolute degrees from north, but what matters is the wind angle relative to the sailing angle. So for each line I need the angle from north as well, and then calculate the difference between the two. That took some drawing of circles on paper with arrows to get it right.
update o_mesh_water
set sailing_angle = st_azimuth(st_startPoint(geom),
st_EndPoint(geom))/pi()*180;
update o_mesh_water set sailing_angle_rev =
CASE
when sailing_angle + 180 < 360 then sailing_angle + 180
else sailing_angle - 180
end;
alter table o_mesh_water
add column windangle numeric,
add column windspeed numeric,
add column real_windangle numeric,
add column sailingspeed numeric,
add column windangle_rev numeric,
add column real_windangle_rev numeric,
add column sailingspeed_rev numeric;
update o_mesh_water
set windangle = winds.direction, windspeed =winds.speed
from winds
where o_mesh_water.nearest_wind_point = winds.id;
-- windangle_rev does not inverse wind direction,
-- but it allows to take the wind info from a different point.
update o_mesh_water
set windangle_rev = winds.direction,
windspeed =winds.speed
from winds
where o_mesh_water.nearest_wind_point_rev = winds.id;
-- using sailing angle reverse from windangle, since if ship goes north,
-- and wind goes north, the relevant angle is 180, not zero.
update o_mesh_water set real_windangle = windangle - sailing_angle_rev;
update o_mesh_water set real_windangle = real_windangle + 360
where real_windangle<0;
update o_mesh_water
set real_windangle_rev = windangle_rev - sailing_angle;
update o_mesh_water
set real_windangle_rev = real_windangle_rev + 360
where real_windangle_rev<0;Next is to translate the winds and direction data to sailing speed. For this we use the information of the polar diagram that we loaded before. Since the polar data is just a table with wind angle in 10 degrees increments, and wind speed at various knots (0, 6, 12, 24, …) we need some conditions to match the more detailed data to this information on the speed of sailing.
-- first for 3 knots
update o_mesh_water o set sailingspeed = sailspeed from
polar p
where
ceiling(round(o.windspeed/1.852)/3)*3=round(p.windspeed_raw/1.852)
and
round(round(o.real_windangle/10))*10=p.angle
and
ceiling(round(o.windspeed/1.852)/3)*3<=6;
-- then the others
update o_mesh_water o set sailingspeed = sailspeed*1.852 from
polar p
where
ceiling(round(o.windspeed/1.852)/6)*6=round(p.windspeed_raw/1.852)
and
round(floor(o.real_windangle/10))*10=p.angle
and
sailingspeed is null;
-- reverse
-- first for 3 knots
update o_mesh_water o set sailingspeed_rev = sailspeed from
polar p
where
ceiling(round(o.windspeed/1.852)/3)*3=round(p.windspeed_raw/1.852)
and
round(floor(o.real_windangle_rev/10))*10=p.angle
and
floor(round(o.windspeed/1.852)/3)*3<=6;
-- then the others
update o_mesh_water o set sailingspeed_rev = sailspeed*1.852 from
polar p
where
ceiling(round(o.windspeed/1.852)/6)*6=round(p.windspeed_raw/1.852)
and
round(floor(o.real_windangle_rev/10))*10=p.angle
and
sailingspeed_rev is null;
alter table o_mesh_water
add column hours_sailing numeric,
add column rev_hours_sailing numeric;
update o_mesh_water
set hours_sailing = (dist/1000)/sailingspeed
where sailingspeed>0;
update o_mesh_water
set hours_sailing = -1
where sailingspeed = 0;
update o_mesh_water
set rev_hours_sailing = -1
where sailingspeed_rev = 0;
update o_mesh_water
set rev_hours_sailing = (dist/1000)/sailingspeed_rev
where sailingspeed_rev>0;
-- limit area down south
update o_mesh_water
set hours_sailing = -1,
rev_hours_sailing= -1
where st_intersects(o_mesh_water.geom,
(select st_union(geom) geom
from
(select 1 as id ,
ST_SetSRID(ST_MakeBox2D(
ST_Point(-40, -70),
ST_Point(155, -45)),4326) geom
union
select 2 as id ,
ST_SetSRID(ST_MakeBox2D(
ST_Point(155, -70),
ST_Point(179.9, -50)),4326) geom
union
select 3 as id,
ST_SetSRID(ST_MakeBox2D(
ST_Point(-179.9, -70),
ST_Point(-100, -45)),4326) geom
union
select 4 as id,
ST_SetSRID(ST_MakeBox2D(
ST_Point(-100, -70),
ST_Point(-40, -60)),4326) geom
) as foo
));Finally, some combination of winds with narrow straights do not allow for two-way navigation. Again this is due to the coarseness of the water grid. So it’s important to manually check some of these straights and ensure two-way navigation where necessary. Gibraltar is the case here, but also check the around Denmark, the Strait of Hormuz, around Indonesia, and the see of Marmara to the Black Sea (recall we created additional grid points there, which was exactly to fix a problem that would otherwise occur here.)2
-- make strait of gibrarlater 2 way navigatable.
update o_mesh_water set rev_hours = hours where
rev_hours=-1 and
st_contains('0102000020E610000002000000000000E08CE818C00000006096D1414000000000F0F614C0000000A02B0B4240',
geom)8 Adding cities
While we have a full global scale of dots from which we can create routes, it would be cumbersome to create routes from each dot in one country to each dot in another country and calculate some average of that.
Instead, I reduce the problems in two steps. 1. limit the routes to only major cities 2. only route from 1 city to the nearest city in four directions.
Then the full route from Netherlands to poland would route from the main cities in the Netherlands and Germany to the main cities in Poland.
The first step is to prepare the cities table, and associate a dot to each city in order to identify their locations on the mesh.
alter table cities alter worldcity type integer;
alter table cities alter adm0cap type integer;
alter table cities rename adm0cap to cap;
alter table cities rename iso_a2 to cntry_code;
alter table cities rename "name" to city;
alter table cities add column node bigint;
alter table cities add column alt real;
select set_srid(geom, 4326) from cities;
delete from cities where latitude<-57 or latitude > 81;
update cities set node = idx.node from (
select gid, node from (
select row_number() over w as r, gid, o.id as node
from cities, o_mesh_vertices_pgr o
where
st_dwithin(cities.geom, o.the_geom, 0.15) -- distance in degrees!
window w AS (
partition by cities.gid order by cities.geom<->o.the_geom
)
) as x
where x.r=1
) as idx
where cities.gid = idx.gid;Next, I add port cities. These are cities that are closely positioned to a coast and can be used to connect to the water grid. In this way the water grid does not just connect at any point to the land grid, but only in specific location of where there are major cities today.
alter table cities add column water_node int;
-- find all cities with maximum 50km away from water grid
with all_coastal_cities as (
select row_number() over w as r, c.gid, c.geom,
c.sov_a3 country_code, g.id node,
st_distance(geography(c.geom), geography(g.the_geom)) dist
from cities c
join
o_mesh_water_vertices_pgr g
on st_dwithin(c.geom, g.the_geom, 1)
window w as (
partition by gid
order by st_distance(geography(c.geom), geography(g.the_geom))
)
),
-- filters out cities close to each other connected to the same node,
-- only taking the closest.
filter_coastal_cities as (
select row_number() over w as s,
gid, geom, country_code, node, dist
from all_coastal_cities a
where a.r=1
window w as (partition by node,country_code order by dist)
),
land_border as (
select gid, st_exteriorRing(geom_land) geom_land
from (select gid, (st_dump(geom)).geom geom_land from land) dmp_land
),
filter_coastal_cities_2 as (
select distinct on (d.gid) d.gid, d.geom, country_code, node, dist
from filter_coastal_cities d
join land_border on
st_dwithin(d.geom, land_border.geom_land, 0.3)
where d.s=1 and dist<=50000
)
update cities set water_node = e.node
from filter_coastal_cities_2 e
where cities.gid=e.gid;
-- add caspian sea cities (they are filtered out due to lake.
-- find all cities with maximum 50km away from water grid
with all_coastal_cities as (
select row_number() over w as r,
c.gid, c.geom, c.sov_a3 country_code, g.id node,
st_distance(geography(c.geom), geography(g.the_geom)) dist
from cities c
join
o_mesh_water_vertices_pgr g on st_dwithin(c.geom, g.the_geom, 1)
window w as (partition by gid
order by st_distance(geography(c.geom), geography(g.the_geom))
)
),
-- filters out cities close to each other connected to the same node,
-- only taking the closest.
filter_coastal_cities as (
select row_number() over w as s,
gid, geom, country_code, node, dist
from all_coastal_cities a
where a.r=1
window w as (partition by node,country_code order by dist)
)
update cities set water_node = e.node
from filter_coastal_cities e
where cities.gid=e.gid and water_node is null;
-- add land-node to these cities
update cities set node = idx.node from (
select gid, node from (
select row_number() over w as r, gid, o.id as node
from cities, o_mesh_vertices_pgr o
where st_dwithin(cities.geom, o.the_geom, 2) -- distance in degrees!
window w AS (partition by cities.gid
order by cities.geom<->o.the_geom
)
) as x where x.r=1
) as idx
where cities.gid = idx.gid and water_node is not null;
-- add mode information
update o_mesh_water set "type"=2;While the idea is nice, and fully automated, some cases, such as Amsterdam, Rotterdam, London (and possibly others) don’t come out very nicely in this method.
The next step is to create explicit links between land nodes and water nodes. These are the links that rule the ‘boarding of a ship’, mode number 3. The cost of this boarding, is set at 24 hours, independent of distance or winds. It’s just assumed that if you want to board or unboard a ship, you’ll have to wait 3 days. This is useful since sailing is usually faster than walking, but a start-up cost for sailing seems like a good way to approximate that sailing is more useful for long voyages.
with land_water_nodes as (
select water_node, water.the_geom water_geom,
node as land_node, land.the_geom as land_geom
from cities
join o_mesh_water_vertices_pgr water on
(cities.water_node = water.id)
join o_mesh_vertices_pgr land on
(cities.node = land.id)
),
water_land_geom as (
select
water_geom, land_geom,
st_setsrid(ST_MakeLine(land_geom::GEOMETRY,
water_geom::GEOMETRY), 4326) AS geom,
st_x(st_startpoint( st_setsrid(ST_MakeLine(land_geom::GEOMETRY,
water_geom::GEOMETRY),
4326))) source_long,
st_y(st_startpoint( st_setsrid(ST_MakeLine(land_geom::GEOMETRY,
water_geom::GEOMETRY),
4326))) source_latt,
st_x(st_endpoint( st_setsrid(ST_MakeLine(land_geom::GEOMETRY,
water_geom::GEOMETRY),
4326))) target_long,
st_y(st_endpoint( st_setsrid(ST_MakeLine(land_geom::GEOMETRY,
water_geom::GEOMETRY),
4326))) target_latt
from land_water_nodes
)
insert into o_mesh_water (
geom, hours, rev_hours, "type", direction, dist
)
select geom, 24 as hours, 24 as rev_hours, 3 as "type",
(
case
when source_long = target_long AND source_latt < target_latt then 1 --'N'
when source_long < target_long AND source_latt < target_latt then 2 --'NE'
when source_long < target_long AND source_latt = target_latt then 3 --'E'
when source_long < target_long AND source_latt > target_latt then 4 --'SE'
when source_long = target_long AND source_latt > target_latt then 5 --'S'
when source_long > target_long AND source_latt > target_latt then 6 --'SW'
when source_long > target_long AND source_latt = target_latt then 7 --'W'
when source_long > target_long AND source_latt < target_latt then 8 --'NW'
else 9 --'ER'
end ) as direction,
st_distance(geography(water_geom), geography(land_geom)) AS dist
from water_land_geom;Then add the entire water mesh into the land mesh.
alter table o_mesh add column "within200km" boolean;
insert into o_mesh (
geom, hours, rev_hours, dist, "type", direction, within200km
)
select
geom,
hours_sailing hours,
rev_hours_sailing rev_hours,
dist,
"type",
direction,
within200km
from o_mesh_water;
create index o_mesh_type_idx on o_mesh("type");
-- somehow cost of type 3 didn't get through, perhaps it does.
-- Otherwise enforce with the following line.
update o_mesh set hours=24, rev_hours=24 where type=3;9 Fix for date line
The date line (+180 / -180 degrees longitude), creates a bit of an issues, as we saw already before. The water mesh and land mesh are not properly connecting the two sides. So, the following is just a fix for this. Essentially it repeats a lot of what is done above, but specifically on the selection of grid points right at the date line.
-- fix for dateline
drop table if exists o_mesh_dl CASCADE;
create table o_mesh_dl
(
-- id BIGSERIAL PRIMARY KEY,
"source" int,
target int,
direction smallint, -- varchar(2),
geom GEOMETRY('LINESTRING', 4326),
dist real,
"type" integer
);
-- grid points concerned Land:
with dl_grid_pos_full as (
--select * from grid where long>179.95
select row_number() over w r, *
from grid where long>179
window w as (partition by latt order by latt, long)
),
dl_grid_pos as (
select * from dl_grid_pos_full where r=1
),
dl_grid_neg as (
select * from grid where long<-179.95)
INSERT INTO o_mesh_dl (source, target, direction, geom, dist, "type")
-- grid_index with next and previous latt value
SELECT
source.pkey as source, target.pkey as target, 7 as direction,
st_setsrid(ST_MakeLine(source.geom::GEOMETRY,
target.geom::GEOMETRY), 4326) AS geom,
st_distance(geography(source.geom),
geography(target.geom)) AS dist,
1 as "type"
FROM
-- distance should match what is made in the grid.
(select * from dl_grid_pos) AS source
JOIN
(select * FROM dl_grid_neg) AS target
ON
-- one way links only
(
target.latt =source.latt
);
-- grid points concerned Water:
with dl_grid_pos as (
select * from grid_water_dl where long > 179.3),
dl_grid_neg as (
select * from grid_water_dl where long< -179.95)
INSERT INTO o_mesh_dl (source, target, direction, geom, dist, "type")
-- grid_index with next and previous latt value
SELECT
source.pkey as source, target.pkey as target, 7 as direction,
st_setsrid(ST_MakeLine(source.geom::GEOMETRY,
target.geom::GEOMETRY), 4326) AS geom,
st_distance(geography(source.geom),
geography(target.geom)) AS dist,
2 as "type"
FROM
-- distance should match what is made in the grid.
(select * from dl_grid_pos) AS source
JOIN
(select * FROM dl_grid_neg) AS target
ON
-- one way links only
(
target.latt =source.latt
);
alter table o_mesh_dl add column new_source int;
alter table o_mesh_dl add column new_target int;
alter table o_mesh_dl add column source_alt real,
add column target_alt real;
update o_mesh_dl set new_source = id
from o_mesh_vertices_pgr
where st_startpoint(geom) && the_geom;
update o_mesh_dl set new_target = id
from o_mesh_vertices_pgr
where st_endpoint(geom) && the_geom;
update o_mesh_dl set source_alt = o.source_alt
from o_mesh o
where o.source = new_source;
update o_mesh_dl set target_alt = o.target_alt
from o_mesh o
where o.target = new_target;
update o_mesh_dl set source_alt = o.target_alt
from o_mesh o
where o.target = new_source and
o_mesh_dl.source_alt is null and
o_mesh_dl."type"=1;
update o_mesh_dl set target_alt = o.source_alt
from o_mesh o
where o.source = new_target and
o_mesh_dl.target_alt is null and
o_mesh_dl."type"=1;
alter table o_mesh_dl add column angle real;
update o_mesh_dl set angle = atan2((target_alt-source_alt),dist);
alter table o_mesh_dl
add column hours real,
add column rev_hours real,
add column hours_tw real,
add column rev_hours_tw real,
add column travers_waters integer;
update o_mesh_dl set hours = (
case
when angle=0. THEN dist/ 5036.742
when abs(angle) <= 12.0/180.0*pi()
THEN (target_alt-source_alt)/sin(angle) /
(1000.0 * 6.0*EXP(-3.5*ABS(TAN(angle)+0.05)))
else (abs(target_alt-source_alt)/sin(12.0/180.0*pi())) /
(1000.0 * 6.0*EXP(-3.5*ABS(TAN(sign(angle)*12.0/180.0*pi())+0.05)))
end
) ,
rev_hours = ( -- only difference is the minus in the dominator
--adjusting speed, distance is the same going up or down.
case
when angle = 0 THEN dist/ 5036.742
when abs(angle) <= 12.0/180.0*pi()
THEN (target_alt-source_alt)/sin(angle) /
(1000.0 * 6.0*EXP(-3.5*ABS(TAN(-1*angle)+0.05)))
-- alternative: (dist/1000)/COS(angle)
else (abs(target_alt-source_alt)/sin(12.0/180.0*pi())) /
(1000.0 * 6.0*EXP(-3.5*ABS(TAN(sign(angle)*-1*12.0/180.0*pi())+0.05)
))
end);
update o_mesh_dl
set hours_tw = hours, rev_hours_tw = rev_hours, travers_waters=0;
update o_mesh_dl
set hours = (dist/1000.) / 15.0, hours_tw =(dist/1000.) / 15.0,
rev_hours = (dist/1000.) / 15.0, rev_hours_tw =(dist/1000.) / 15.0
where "type"=2;
-- indicator column for if something goes wrong
alter table o_mesh add column dl int;
insert into o_mesh (source, target, direction, geom, dist,
source_alt, target_alt, hours, rev_hours,
type, travers_waters, hours_tw, rev_hours_tw,
angle, dl)
select new_source as source, new_target as target,
direction, geom, dist, source_alt, target_alt,
hours, rev_hours, type, travers_waters, hours_tw,
rev_hours_tw, angle, 1
from o_mesh_dl; Housekeeping
-- indicator column can be deleted.
vacuum full analyze o_mesh;
drop table if exists o_mesh_vertices_pgr cascade;
select pgr_createTopology('o_mesh', 0.0001, 'geom', 'pkey');
vacuum full analyze o_mesh;
vacuum analyze;10 Routing
Basically our mesh is now ready to route through. Choose an origin and destination and the PGRouting algorithm will give you the optimal route, with costs (travel time). But we need this for many combinations, so it requires automation.
10.1 Preparation
First update the information in the cities table to match the nodes information. This identification number may have changed with the pgr_createTopology call before.
update cities set node = idx.node from (
select gid, node from (
select row_number() over w as r, gid, o.id as node from cities,
o_mesh_vertices_pgr o
where st_dwithin(cities.geom, o.the_geom, 2) -- distance in degrees!
window w AS (partition by cities.gid
order by cities.geom<->o.the_geom)) as x
where x.r=1 ) as idx
where cities.gid = idx.gid;Next, in preparation for what’s coming, we create a new column type, called route_rec (route record), which has three pieces of information in one, the cost in hours, the cost in distance and the actual geospatial information of the route.
create type route_rec as (
hours double precision,
dist double precision,
geom geometry
); Then we prepare a table that will hold all routes that we wish to calculate. For each city we want information on the origin, frm (from), and to.
drop table if exists routes cascade;
create table routes (
id serial primary key,
frm_gid int, --from city id
to_gid int, -- to city id
frm_node int, -- node ide
to_node int,
frm_city varchar(100), -- city name
to_city varchar(100),
frm_cntry_code varchar(5), -- country code
to_cntry_code varchar(5),
frm_cap int, -- city is national capital
to_cap int,
dist_geodesic double precision, -- geodesic distance
dist_walking double precision, -- distance walking/sailing
hours_walking double precision, -- hours waling/sailing
geom_geodesic GEOMETRY, -- geospatial geodesic
geom_walking GEOMETRY, -- geospatial routing
direct int,
pair_id int,
res_route route_rec
);
create index routes_frm_gid_idx on routes(frm_gid);
create index routes_to_gid_idx on routes(to_gid);
create index routes_geom_geodesic_idx on routes using gist(geom_geodesic);
create index routes_geom_walking_idx on routes using gist(geom_walking);
create index routes_pair_id_idx on routes(pair_id);
alter table routes add constraint unq_gid UNIQUE(frm_gid,to_gid);The main selection for included cities is that they are of specific size, as given in the cities table through scalerank or a national capital. This gives 1138 cities (out of 7280). If we’d routes from each city to any other cities we’d get 1138\(\times\)(1138-1) \(\times\) 2 = 2.59 million routes, many of which will be greatly overlapping. So in order to avoid this, we route only to the nearest cities in four direction (north, east, south west).
with all_combinations as (
select row_number() over () id,
frm.gid frm_gid, "to".gid to_gid,
frm.node frm_node, "to".node to_node,
frm.city frm_city, "to".city to_city,
frm.cntry_code frm_cntry_code, "to".cntry_code to_cntry_code,
frm.longitude frm_long, frm.latitude frm_lat,
"to".longitude to_long, "to".latitude to_lat,
frm.geom frm_geom, "to".geom to_geom,
frm.cap frm_cap, "to".cap to_cap
from cities "to"
join cities "frm" on
(("frm".cap = 1
or frm.scalerank <5)
and ("to".cap = 1
or "to".scalerank <5)
and frm.gid<>"to".gid)),
-- for negative long and lat these quadrants are flipped,
-- but that shouldn't matter for the final combined selection
quadrant_north as (
select row_number() over w as r, *,
st_distance(geography(frm_geom), geography(to_geom))
from all_combinations
where
(abs(to_lat)-abs(frm_lat))>= abs(abs(to_long)-abs(frm_long)) and
st_dwithin(geography(frm_geom), geography(to_geom), 1500000)
window w as (partition by frm_gid
order by st_distance(geography(frm_geom), geography(to_geom))
)
),
quadrant_east as (
select row_number() over w as r, *,
st_distance(geography(frm_geom), geography(to_geom))
from all_combinations
where
(abs(to_long)-abs(frm_long))>= abs(abs(to_lat)-abs(frm_lat)) and
st_dwithin(geography(frm_geom), geography(to_geom), 1500000)
window w as (partition by frm_gid
order by st_distance(geography(frm_geom), geography(to_geom))
)
),
quadrant_south as (
select row_number() over w as r, *,
st_distance(geography(frm_geom), geography(to_geom))
from all_combinations
where
(abs(frm_lat)-abs(to_lat))>= abs(abs(to_long)-abs(frm_long)) and
st_dwithin(geography(frm_geom), geography(to_geom), 1500000)
window w as (partition by frm_gid
order by st_distance(geography(frm_geom), geography(to_geom))
)
),
quadrant_west as (
select row_number() over w as r, *,
st_distance(geography(frm_geom), geography(to_geom))
from all_combinations
where
(abs(frm_long)-abs(to_long))>= abs(abs(to_lat)-abs(frm_lat)) and
st_dwithin(geography(frm_geom), geography(to_geom), 1500000)
window w as (partition by frm_gid
order by st_distance(geography(frm_geom), geography(to_geom))
)
),
combination as (
select * from quadrant_north where quadrant_north.r<3
union
select * from quadrant_east where quadrant_east.r<3
union
select * from quadrant_south where quadrant_south.r=1
union
select * from quadrant_west where quadrant_west.r=1)
insert into routes (frm_gid, to_gid, frm_node, to_node,
frm_city, to_city, frm_cntry_code, to_cntry_code, frm_cap, to_cap)
select distinct on (frm_gid, to_gid) frm_gid, to_gid,
frm_node, to_node, frm_city, to_city,
frm_cntry_code, to_cntry_code, frm_cap, to_cap
from combination;
-- match only when they frm and to are both optimal?
with idx as (
select distinct on (greatest(frm_gid, to_gid), least(frm_gid, to_gid))
row_number() over () as id, frm_gid, to_gid
from routes
)
update routes r set pair_id = idx.id
from idx
where (r.frm_gid=idx.frm_gid and r.to_gid=idx.to_gid) or
(r.to_gid=idx.frm_gid and r.frm_gid=idx.to_gid);
-- then select routes with where clause
delete from routes where
pair_id not in (
select pair_id from (
select pair_id, count(*) cnt from routes group by pair_id
) cnt where cnt.cnt=2
);While a convenient starting point, it turns out that quite a few crucial links are not included. This could be because of the limit of nearest city that was set above (150000 meter = 150km), or some other issues that just need some manual fixing.
-- delete some routes with missing data
delete from routes where frm_node is null or to_node is null;
-- add some
insert into routes (frm_gid, to_gid, frm_node, to_node,
frm_city, to_city, frm_cntry_code, to_cntry_code, frm_cap, to_cap)
select
frm.gid frm_gid, "to".gid to_gid, frm.node frm_node,
"to".node to_node, frm.city frm_city, "to".city to_city,
frm.cntry_code frm_cntry_code, "to".cntry_code to_cntry_code,
frm.cap frm_cap, "to".cap to_cap
from cities "to", cities "frm"
where
("to".gid=6579 and frm.gid=3087) OR (frm.gid=6579 and "to".gid=3087)
OR ("to".gid=3087 and frm.gid=5196) OR (frm.gid=3087 and "to".gid=5196)
OR ("to".gid=5196 and frm.gid=6951) OR (frm.gid=5196 and "to".gid=6951)
OR ("to".gid=6951 and frm.gid=2158) OR (frm.gid=6951 and "to".gid=2158)
OR ("to".gid=2158 and frm.gid=6836) OR (frm.gid=2158 and "to".gid=6836)
-- baring street
OR ("to".gid=7215 and frm.gid=7281) OR (frm.gid=7215 and "to".gid=7281)
-- Japan Hiroshima to Kobe
OR ("to".gid=7157 and frm.gid=7166) OR (frm.gid=7157 and "to".gid=7166)
-- Arabian Peninsula (sanaa - mecca)
OR ("to".gid=7157 and frm.gid=6914) OR (frm.gid=7157 and "to".gid=6914)
-- Arabian Peninsula (sanaa - muscat)
OR ("to".gid=6180 and frm.gid=5833) OR (frm.gid=6180 and "to".gid=5833)
-- Italy (Monaco - Genua)
OR ("to".gid=5833 and frm.gid=7182) OR (frm.gid=5833 and "to".gid=7182)
-- Italy (Genua - Milan)
OR ("to".gid=5833 and frm.gid=3379) OR (frm.gid=5833 and "to".gid=3379)
-- Italy (Genua - Pisa)
OR ("to".gid=3379 and frm.gid=6485) OR (frm.gid=3379 and "to".gid=6485)
-- Italy (Pisa - Florence)
OR ("to".gid=6315 and frm.gid=6314) OR (frm.gid=6315 and "to".gid=6314)
-- Turkey (Adana - Konya)
OR ("to".gid=7300 and frm.gid=7022) OR (frm.gid=7300 and "to".gid=7022)
-- Turkey (Istanbul - Tessaloniki)
OR ("to".gid=7180 and frm.gid=7237) OR (frm.gid=7180 and "to".gid=7237)
-- West African Coast (Monrovia - Abidjan)
OR ("to".gid=7237 and frm.gid=7221) OR (frm.gid=7237 and "to".gid=7221)
-- West African Coast (Abidjan - Accra)
OR ("to".gid=7006 and frm.gid=6671) OR (frm.gid=7006 and "to".gid=6671)
-- West African Coast (Port Gentil - Point-Noire)
OR ("to".gid=6671 and frm.gid=7250) OR (frm.gid=6671 and "to".gid=7250)
-- West African Coast (Point-Noire - Luanda)
OR ("to".gid=6375 and frm.gid=7170) OR (frm.gid=6375 and "to".gid=7170)
-- South African Coast (East-London - Durban)
OR ("to".gid=7170 and frm.gid=6899) OR (frm.gid=7170 and "to".gid=6899)
-- South African Coast (Durban - Maputo)
OR ("to".gid=6374 and frm.gid=6887) OR (frm.gid=6374 and "to".gid=6887)
-- East African Coast (Quelimane - Nacala)
OR ("to".gid=6966 and frm.gid=6913) OR (frm.gid=6966 and "to".gid=6913)
-- East African Coast (Mombasa - Mogadishu)
OR ("to".gid=6482 and frm.gid=6484) OR (frm.gid=6482 and "to".gid=6484)
-- Indonesia (Ambon - Jayapura)
OR ("to".gid=6482 and frm.gid=5827) OR (frm.gid=6482 and "to".gid=5827)
-- Indonesia (Ambon - merauke)
OR ("to".gid=6484 and frm.gid=5827) OR (frm.gid=6484 and "to".gid=5827)
-- Indonesia (Jayapura - merauke)
OR ("to".gid=6484 and frm.gid=5827) OR (frm.gid=6484 and "to".gid=5827)
-- Indonesia (Jayapura - merauke)
OR ("to".gid=5827 and frm.gid=5622) OR (frm.gid=5827 and "to".gid=5622)
-- Indonesia (merauke - weipa)
OR ("to".gid=6645 and frm.gid=6297) OR (frm.gid=6645 and "to".gid=6297)
-- Caspian Sea (Atyrau - Nukus)
OR ("to".gid=7096 and frm.gid=7024) OR (frm.gid=7096 and "to".gid=7024)
-- Caspian Sea (Baku - Tblisi)
OR ("to".gid=7024 and frm.gid=5375) OR (frm.gid=7024 and "to".gid=5375)
-- Black Sea (Tblisi - Batumi)
OR ("to".gid=5375 and frm.gid=6389) OR (frm.gid=5375 and "to".gid=6389)
-- Black Sea (Batumi - Sochi)
OR ("to".gid=5375 and frm.gid=6313) OR (frm.gid=5375 and "to".gid=6313)
-- Black Sea (Batumi - Samsun)
OR ("to".gid=6295 and frm.gid=6645) OR (frm.gid=6295 and "to".gid=6645)
-- Caspian Sea (Turkmenbase - Atyrau)
;
delete from routes where
(to_gid=6756 and frm_gid = 6852) -- australia
OR (to_gid=7153 and frm_gid = 6302) -- thailand
;Next, we need to allow for global sailing. Essentially this allows to go from any city to the nearest port and then to the final destination using the nearest port there. This will create a mess again, so we use the assumption that to go from/to Europe to/from asia one has to dock in Cape Good Hope in South Africa.
Using this reasoning we add the following combinations that will be used for global sailing.
-- sailing routes:
with sel_cities as (
select row_number() over w as r, gid, city, cntry_code, node, cap
from cities where (cap=1 or scalerank<7) and water_node is not null
window w as (partition by cntry_code order by pop_max DESC)
),
combinations as (
-- europe with americas
select frm.gid frm_gid, "to".gid to_gid,
frm.node frm_node, "to".node to_node,
frm.city frm_city, "to".city to_city,
frm.cntry_code frm_cntry_code, "to".cntry_code to_cntry_code,
frm.cap frm_cap, "to".cap to_cap
from sel_cities "frm", sel_cities "to"
where
(
-- europe to west africa
(
"frm".cntry_code in ('NL', 'GB', 'FR', 'ES', 'IT', 'PT',
'IE', 'DK', 'NO') and
"to".cntry_code in ('MA', 'MR', 'SM', 'BJ', 'GW', 'GN', 'SL',
'LR', 'GH', 'TG', 'BJ', 'NG', 'CQ', 'GA',
'CG', 'AO', 'NA', 'ZA')
)
OR
(
"to".cntry_code in ('NL', 'GB', 'FR', 'ES', 'IT', 'PT',
'IE', 'DK', 'NO') and
"frm".cntry_code in ('MA', 'MR', 'SM', 'BJ', 'GW', 'GN', 'SL',
'LR', 'GH', 'TG', 'BJ', 'NG', 'CQ', 'GA',
'CG', 'AO', 'NA', 'ZA')
)
OR
-- east africa/middle east to south asia
(
"frm".cntry_code in ('ZA', 'MZ', 'TZ', 'KE', 'SO', 'DJ',
'ER', 'SD', 'SA', 'YE', 'MG', 'BH',
'QA', 'AE', 'OM', 'IR') and
"to".cntry_code in ('IN', 'LK', 'BD', 'MM', 'MY', 'SG',
'ID', 'PG', 'AU', 'NZ')
)
OR
(
"to".cntry_code in ('ZA', 'MZ', 'TZ', 'KE', 'SO', 'DJ',
'ER', 'SD', 'SA', 'YE', 'MG', 'BH',
'QA', 'AE', 'OM', 'IR') and
"frm".cntry_code in ('IN', 'LK', 'BD', 'MM', 'MY', 'SG',
'ID', 'BN', 'PG', 'AU', 'NZ')
)
OR
-- south to east asia
(
"frm".cntry_code in ('ID') and
"to".cntry_code in ('VN', 'CN', 'TW', 'JP', 'KR', 'RU')
)
OR
(
"to".cntry_code in ('ID') and
"frm".cntry_code in ('VN', 'CN', 'TW', 'JP', 'KR', 'RU')
)
OR
-- europe to east americas
(
"frm".cntry_code in ('NL', 'GB', 'FR', 'ES', 'IT', 'PT',
'IE', 'DK', 'NO') and
"to".cntry_code in ('US', 'CU', 'HT', 'GT', 'CR', 'CO', 'VE',
'TT', 'GY', 'BR', 'UY', 'AR')
)
OR
(
"to".cntry_code in ('NL', 'GB', 'FR', 'ES', 'IT', 'PT', 'IE',
'DK', 'NO') and
"frm".cntry_code in ('US', 'CU', 'HT', 'GT', 'CR', 'CO',
'VE', 'TT', 'GY', 'BR', 'UY', 'AR')
)
OR
-- west americas to east asia
(
"frm".cntry_code in ('CA', 'MX', 'NI', 'EC', 'PE', 'CL') and
"to".cntry_code in ('VN', 'CN', 'TW', 'JP', 'KR', 'RU')
)
OR
(
"to".cntry_code in ('CA', 'MX', 'NI', 'EC', 'PE', 'CL') and
"frm".cntry_code in ('VN', 'CN', 'TW', 'JP', 'KR', 'RU')
)
)
and
"frm".r=1 and "to".r=1
)
insert into routes (frm_gid, to_gid, frm_node, to_node,
frm_city, to_city, frm_cntry_code, to_cntry_code,
frm_cap, to_cap, with_sailing)
select frm_gid, to_gid, frm_node, to_node, frm_city, to_city,
frm_cntry_code, to_cntry_code, frm_cap, to_cap, TRUE as with_sailing
from combinations;
-- update pair_id
with index as (
select distinct on (greatest(frm_gid, to_gid), least(frm_gid, to_gid))
row_number() over () as id, frm_gid, to_gid from routes
)
update routes r set pair_id = index.id from index
where
(r.frm_gid=index.frm_gid and r.to_gid=index.to_gid) or
(r.to_gid=index.frm_gid and r.frm_gid=index.to_gid);Housekeeping
alter table routes
drop column id,
add column id serial primary key;We finally have 970,000+ routes to calculate.
10.2 Geodesic distance
Geodesic distance is the easiest to calculate and incorporate, so we’ll start with this.
update routes set dist_geodesic=dist, geom_geodesic=geom from (
select
frm.gid as frm_gid, "to".gid as to_gid,
st_distance(geography(frm.geom), geography("to".geom)) dist,
ST_Segmentize(st_makeline(frm.geom ,"to".geom)::geography,
greatest(1, round(st_distance(geography(frm.geom),
geography("to".geom))/10000))
)::geometry geom
from
cities as frm,
cities as "to"
) idx
where idx.frm_gid = routes.frm_gid and idx.to_gid=routes.to_gid;
create index routes_geom_geod on routes using gist(geom_geodesic);The distance measure will be a useful benchmark to compare the walking/sailing routes against. The geometry that is also created to show the route is cute, but not very useful otherwise.
10.3 Routing
As indicated quite early, we want to make a distinction between two era’s, one where global sailing wasn’t possible, and one where it was. The technical distinction is in cost of travel in the routing network. The cost are exactly the same, except that points more than 200km from a coastline are inaccessible. Instead of completely duplicating the routing table, we create a view where this condition is incorporated.
drop view if exists o_mesh_ancient;
create view o_mesh_ancient as
(
select
source, target, dist, geom,
(case
when "type"=1 THEN hours_tw
when "type"=2 AND within200km=TRUE THEN hours
when "type"=3 THEN hours
else -1
end) as hours,
(case
when "type"=1 THEN rev_hours_tw
when "type"=2 AND within200km=TRUE THEN rev_hours
when "type"=3 THEN rev_hours
else -1
end) as rev_hours,
pkey
from o_mesh);Next we create two functions, that perform the actual routing algorithm and throw back the required information in the route_rec, given the two nodes that we defined in the routes table.
The reason that we need two functions is again due to the dateline and an optimisation that is incorporated in the first function that causes a issues with that. The optimisation is that given the location of the origin node and destination node, we draw a box around it and find an optimal route within that box. This speeds-up the optimisation considerably because going from Amsterdam to Rotterdam, we do not need to consider any possible route around the globe. However, this ‘box’ can create issues when there are large water bodies, where the sailing constraints kick in, or again with the date line.
So the first function deals with the easy routes quickly, and the second function does the remaining.
-- function that calculates the distance, hours (cost) and
-- geometry between 2 nodes
create or replace function full_route(n1 integer, n2 integer)
returns route_rec as
$$
DECLARE
result_record route_rec;
BEGIN
with full_route AS (
SELECT
route.seq AS seq,
route.id1 AS node,
route.id2 AS edge,
route."cost",
o_m.dist,
o_m.geom
FROM pgr_dijkstra(
format('
SELECT r.pkey as id,
r.source::INT,
r.target::INT,
r.hours::DOUBLE PRECISION AS "cost",
r.rev_hours::double precision AS reverse_cost
FROM o_mesh_ancient as r
JOIN
-- the box will cause the Mediterranean/Baltic sea etc
-- in particular to be an obstacle for walking.
-- add minimum height and width.
(SELECT ST_Expand(ST_Extent(the_geom),2) as box
FROM o_mesh_vertices_pgr v WHERE v.id=%1$s
OR v.id=%2$s ) as box on box.box && r.geom'
, n1, n2
)::TEXT,
n1,n2, true , true) AS route
LEFT JOIN
o_mesh_ancient o_m on route.id2=o_m.pkey
)
SELECT
sum(cost) ,
sum(dist) ,
st_multi(st_union(geom))
INTO
result_record.hours, result_record.dist, result_record.geom
FROM full_route;
RETURN result_record;
END $$
LANGUAGE plpgsql;Note, that this functions refers to o_mesh_ancient. So will not route across oceans. The function needs to be redefined to use o_mesh for the whole network, including sailing to be used.
Second function just omits the box part.
create or replace function full_route_dl(n1 integer, n2 integer)
returns route_rec as
$$
DECLARE
result_record route_rec;
BEGIN
with full_route AS (
SELECT
route.seq AS seq,
route.id1 AS node,
route.id2 AS edge,
route."cost",
o_mesh.dist,
o_mesh.geom
FROM pgr_dijkstra(
'
SELECT r.pkey as id,
r.source::INT,
r.target::INT,
r.hours::DOUBLE PRECISION AS "cost",
r.rev_hours::double precision AS reverse_cost
FROM o_mesh_ancient as r
where type!=1 -- only sailing
-- where pkey in (select pkey from o_mesh_dl_sel)
'::TEXT,
n1,n2, true , true) AS route
LEFT JOIN
o_mesh on route.id2=o_mesh.pkey
)
SELECT
sum(cost) ,
sum(dist) ,
st_multi(st_union(geom))
INTO
result_record.hours, result_record.dist, result_record.geom
FROM full_route;
RETURN result_record;
END $$
LANGUAGE plpgsql;then run:
-- actual query can be split between processes, in sets of 800
update routes set res_route = full_route(frm_node, to_node)
where id>=(1-1)*800 and id<1*800 and with_sailing=FALSE;
-- mop up the remaining
update routes set res_route = full_route_dl(frm_node, to_node)
where id in (select id from routes where res_route is null);This saves all information for each row in the res_route column. This can be quickly unpacked into the costs and geometry columns:
update routes
set dist_walking = (res_route).dist,
hours_walking = (res_route).hours,
geom_walking = (res_route).geom::GEOMETRY('MULTILINESTRING', 4326);Then repeat the above for sailing. Basically this can be limited to the a subset of the routes that require sailing only, since the routes based on land only would not have changed. Before running this command, the second function needs to be adjusted to use o_mesh instead of o_mesh_ancient.
update routes set res_route = full_route_dl(frm_node, to_node)
where with_sailing=TRUE;The assumption of the above is that no pair of cities would have the possibility to have a route without cross-ocean sailing, and one with. So we only have to keep track of the city pairs that require cross-ocean sailing, and bring that forward in the following.
Housekeeping
create index routes_frm_node_idx on routes(frm_node);
create index routes_to_node_idx on routes(to_node);
create index routes_geom_walking_idx on routes using gist(geom_walking);The ratio of geodesic over optimal routes can give a measure of detour required by route. This can be useful for descriptive statistics and for plotting. The benchmark speed is 5036.742 m/hour, which is the number when the angle of walking is 0. Sailing will usually be much faster then walking, so the sailing routes will have a ratio below 1.
alter table routes
add column rat_dist numeric,
add column rat_hours numeric;
update routes
set rat_dist = dist_walking/dist_geodesic,
rat_hours = hours_walking/ (dist_walking/5036.742);Every route is there twice. We can take the one with the highest value for hours, since presumably that’s the harder route which may be easier in the other direction.
drop view if exists routes_sel;
create view routes_sel as (
select pair_id, id, frm_gid, to_gid, frm_node, to_node,
frm_city, to_city frm_cntry_code, to_cntry_code,
frm_cap, to_cap, dist_geodesic, dist_walking, hours_walking,
geom_walking, geom_geodesic, direct, rat_dist, rat_hours from (
select row_number() over w as r, *
from routes
where geom_walking is not null
window w as (partition by pair_id order by hours_walking DESC)
) as foo
where r=1
);
alter table routes add column sel boolean;
update routes set sel =FALSE;
update routes set sel =TRUE where id in (select id from routes_sel);The routing network is now done, and can be represented in a global map, or output as a shapefile for further processing.
pgsql2shp -f '<file+path to save to>' -u wessel empires /
"select id, frm_gid, to_gid, frm_city, to_city, frm_cntry_code,
to_cntry_code, dist_geodesic, dist_walking,
geom_walking, hours_walking, rat_dist, rat_hours, with_sailing
from routes
where geom_walking is not null and sel=TRUE"11 Actual routing
The city to city routes now represent the actual network. The nodes links to each other, and each route has a cost associated to it. So now to route from Amsterdam to Being, the PGrouting algorithm can use the city-to-city network to find the optimal route across the globe, either walking, or with sailing.
Note, there is no need to create topology (i.e. pgr_createTopology), it will only mess things up.
Add routes nodes to the cities table
alter table cities add column routes_node int;
update cities set routes_node = idx.node from (
select gid, node from (
select row_number() over w as rid, c.gid, r.id as node
from
cities c join o_mesh_vertices_pgr o on (c.node = o.id),
routes_vertices_pgr r
where st_dwithin(o.the_geom, r.the_geom, 0.1) -- distance in degrees!
window w AS (partition by o.id order by o.the_geom<->r.the_geom)) as x
where x.rid=1 ) as idx
where cities.gid = idx.gid;Create a table for all city combinations between countries. See the where constraint, where the city combinations must come from different countries.
drop table if exists between_all;
create table between_all as (
with gids as (
select distinct on (gid) gid from
(select frm_gid as gid from
routes union select to_gid as gid from routes) as foo
)
select
f_c.gid frm_gid,
f_c.city frm_city,
f_c.cntry_code frm_cntry_code,
f_c.node::int frm_node,
t_c.gid to_gid,
t_c.city to_city,
t_c.cntry_code to_cntry_code,
t_c.node::int to_node
from cities f_c, cities t_c
where
f_c.gid <> t_c.gid and
f_c.gid in (select gid from gids) and
t_c.gid in (select gid from gids)
and f_c.cntry_code <> t_c.cntry_code
);
alter table between_all
add column result_routing route_rec,
add column id serial primary key;Then a new function to execute these routes using this network
create or replace function full_route_reduced(n1 int, n2 int)
returns route_rec as
$$
DECLARE
result_record route_rec;
BEGIN
with full_route as (
SELECT
route.seq AS seq,
route.id1 AS node,
route.id2 AS edge,
route."cost" as cost,
routes.dist_walking as dist
FROM pgr_dijkstra('
SELECT r.id as id,
r.frm_node::INT source,
r.to_node::INT target,
r.hours_walking::DOUBLE PRECISION AS "cost",
-1::double precision AS reverse_cost
FROM routes as r
--where with_sailing=FALSE
'
, n1::int,n2::int, FALSE , FALSE) AS route
LEFT JOIN routes on (route.id2=routes.id)
)
select sum(cost) hours, sum(dist)
into result_record.hours, result_record.dist
from full_route;
RETURN result_record;
EXCEPTION when others then
raise notice ' % - % fails', n1, n2;
RETURN NULL;
END $$
LANGUAGE plpgsql;The function now allows for cross-ocean sailing. Simply uncomment the where clause that has the condition with_sailing=FALSE, to prevent sailing routes in the network.
and roll (again using multiple processes and increment the id numbers with a certain amount:
UPDATE between_all
set result_routing=full_route_reduced(frm_node, to_node)
where id>=(1-1)*140000 and id<1*140000; --range for multi-processOn my machine, takes 4500 rows take ~30 minutes (164999 = 51527352.88ms).
Housekeeping
alter table between_all
add column dist double precision,
add column hours double precision;
update between_all
set dist = (result_routing).dist,
hours = (result_routing).hours;
alter table between_all add column res_rout_sail route_rec;
-- change function above to include sailing
UPDATE between_all
set res_rout_sail=full_route_reduced(frm_node, to_node)
where id>=(1-1)*140000 and id<1*140000; -- range for multi proces
alter table between_all
add column dist_sail double precision,
add column hours_sail double precision;
update between_all
set dist_sail = (res_rout_sail).dist,
hours_sail = (res_rout_sail).hours;
alter table between_all add column dist_bird double precision;
alter table between_all
add foreign key(frm_node) references o_mesh_vertices_pgr(id),
add foreign key(to_node) references o_mesh_vertices_pgr(id);Add geodesic distance, called dist_bird.
with d as (
select b.id, b.frm_node, b.to_node,
st_distance(GEOGRAPHY(f.the_geom), GEOGRAPHY(t.the_geom)) dist
from
between_all b left join o_mesh_vertices_pgr f on frm_node=f.id
left join o_mesh_vertices_pgr t on to_node=t.id
)
update between_all set dist_bird = d.dist
from d where between_all.id=d.id;Create an indicator for capital-to-capital distance by matching against the original cities table and taking the capital flag from there.
alter table between_all
add column frm_cap boolean,
add column to_cap boolean;
with d as (
select b.id, b.frm_gid, b.to_gid,
f.cap::boolean frm_cap, t.cap::boolean to_cap,
f.sov_a3 frm_cntry_code, t.sov_a3 to_cntry_code
from
between_all b left join cities f on frm_gid=f.gid
left join cities t on to_gid=t.gid
)
update between_all set frm_cap = d.frm_cap, to_cap = d.to_cap ,
frm_cntry_code=d.frm_cntry_code, to_cntry_code=d.to_cntry_code
from d where between_all.id=d.id;Graphical representations of these routes are provided in the journal article.
The tabular results of the distances can be exported to a csv file for further processing and inclusion in your favourite statistics package
\copy between_all (id, frm_gid, to_gid, frm_city,
frm_cap, frm_cntry_code, to_city,
to_cap, to_cntry_code, dist, hours,
dist_sail, hours_sail, dist_bird)
to '<path+file>' WITH CSV DELIMITER ',' HEADEROr if interested in country-to-country average only, you can collapse in Postgres already.
drop table if exists between_all_cntry;
create table between_all_cntry as (
select frm_cntry_code, to_cntry_code,
avg(hours) hours,
avg(dist) dist,
avg(hours_sail) hours_sail,
avg(dist_sail) dist_sail,
avg(dist_bird) dist_bird
from between_all
group by frm_cntry_code, to_cntry_code
);and output that
\copy between_all_cntry (frm_cntry_code, to_cntry_code,
dist, hours, dist_sail, hours_sail, dist_bird)
to '<path+file>' WITH CSV DELIMITER ',' HEADERDone.