Thursday, 18 April 2013

MapGuide and pgRouting: A shortest path example

Ok folks, here's how you do this

For a little back-story, since reading Jo Cook's pgRouting example with MapServer, I wanted to do an equivalent one for MapGuide just to see for myself (and many others) if the MapGuide and pgRouting combination was indeed possible.

This blog post will show you a basic PHP proof of concept example of pgRouting integration with MapGuide using the following stack:
  • MapGuide Open Source 2.5
  • PostgreSQL 9.2
  • PostGIS 2.0.3
  • pgRouting 1.0.7
For different versions of PostGIS and/or pgRouting, you may have to adjust some of the code shown here.

This example has only been tested on Windows, but assuming you satsify the above pre-requisites I couldn't see why this example won't work for Linux as well.

You can find the code and sample package here for reference. This example works with both the AJAX and Fusion viewers

Getting the source data

This example uses the GeoBase "National Road Network" dataset for British Columbia (nrn_rrn_bc_shp_en.zip). This is the same reference dataset used in the pgRouting tutorial.

Creating the PostGIS database

Using pgAdmin or the command-line tools, create a new database "pgrouting_example" using "template_postgis_20". 

In case you're wondering, FDO Toolbox does not support creating a PostGIS 2.0 data store because the provider's data store and spatial index creation logic hasn't been updated for PostGIS 2.0. If you're using the PostgreSQL FDO provider that has this changeset compiled in, then you have this support. This change didn't make it for the PostgreSQL provider in FDO 3.8, so as it stands MapGuide Open Source 2.5 can connect to PostGIS 2.0 databases but not create them, which is still okay for this example as we're using external tools to set it up.

Enabling pgRouting

The pgrouting_example database is PostGIS-enabled, but not pgRouting-enabled. Execute this SQL on that database to enable pgRouting support:

CREATE EXTENSION pgrouting

Loading our Roads Shape File

We use the shp2pgsql utility to load our roads shape file into PostGIS

Why don't we eat our own dogfood and use FDO Toolbox to do this? Again, some shortcomings:
  • FDO Toolbox bulk copies geometries as-is, multi-linestrings don't get converted to linestrings. This is important for pgRouting in terms of building a routable network as you'll see in the next step
  • FDO Toolbox lacks PostGIS 2.0 support in its bundled FDO provider (as mentioned above).
So for now, we'll have to do with shp2pgsql

For our example dataset, the shape file in question is NRN_BC_9_0_ROADSEG.shp. The necessary command to load this into our PostGIS database (using the postgres db account) is:

shp2pgsql -s 4617 -I -S NRN_BC_9_0_ROADSEG.shp roads | psql -d pgrouting_example -U postgres

4617 is the SRID to load these geometries in. The "-S" switch instructs the tool to generate simple geometries instead of MULTI geometries, ensuring our lines are put into PostGIS as LINESTRING geometries and not MULTILINESTRING geometries. You'll see why this is important in the next step.

Making this data routable

We now need to extend our imported roads table to support pgRouting, we need to add the following columns:
  • source (integer, the source node id)
  • target (integer, the target node id)
  • length (double precision, the geometric length of the road line segment)
The following SQL creates these columns

alter table roads add column source integer;
alter table roads add column target integer;
alter table roads add column length double precision;

Then we use the pgr_assign_vertex_id function in pgRouting to populate the columns we just added. This function only works on linestrings and not multi-linestrings. This is why we use shp2pgsql, which can do this necessary conversion. This function was called assign_vertex_id in older versions of pgRouting. We call it like so:

select pgr_assign_vertex_id('roads', 0.00000001, 'geom', 'gid');

The 0.00000001 value is the "snap tolerance", since this dataset is in latitudes and longitudes (thus degree-based units) we need a really fine-grained value. This value varies from dataset to dataset.

This command will take a few minutes. Then we need to set the length column:

update roads set length = st_length(geom);

And finally, apply the necessary indexes for good performance (no need for a spatial index, as shp2pgsql already created one):

create index source_idx on roads(source);
create index target_idx on roads(target);

Consuming PostGIS from MapGuide

As mentioned, as of MapGuide Open Source 2.5, the bundled PostgreSQL FDO provider supports connecting to PostGIS 2.0 databases. It just can't create them, which is okay because the MapGuide API does not support Feature Source creation for RDBMS providers anyways. So setting up the Feature Source and layers is fairly straight-forward. 

But to save time and words on this post, just load in the referenced sample package, and tweak the credentials of Library://pg_routing/Data/postgis.FeatureSource to match your PostgreSQL credentials as required.

The sample package has two maps, a normal map (for consumption with the AJAX viewer) and a "Web Mercator" projected version (for consumption with Fusion with an OpenStreetMap backdrop).

Installing the sample application

Assuming the standard installation directory of C:\Program Files\OSGeo\MapGuide, extract the sample application to the C:\Program Files\OSGeo\MapGuide\Web\www directory.

The extracted directory should look like this

Running the sample application

The sample application works with both the AJAX viewer and Fusion viewer. Load them like so:

AJAX Viewer: 

http://localhost/mapguide/mapviewerajax/?WEBLAYOUT=Library://pg_routing/Layouts/pgrouting_basic.WebLayout&LOCALE=en

Fusion: 

http://localhost/mapguide/fusion/templates/mapguide/slate/index.html?ApplicationDefinition=Library://pg_routing/Layouts/pgrouting_fusion.ApplicationDefinition&locale=en

Plug in the port numbers as necessary into the URL for Apache configurations

The Demo

Regardless of the viewer, the Task Pane will be loaded with the main page of the sample application.


Now this is a proof of concept, so the UI is very simplistic and barebones and requires you to enter the source and target node ids to compute the shortest path. Ideally, you would record the start/end points with the digitizing viewer API, tap into PostGIS trickery server-side to "snap" to the nearest road segments, extract the relevant source and target IDs and calculate shortest path from there.

In our case, this has to be done manually, so select the road segment closest to your start point and observe its "source" property value.


Repeat for the selected road segment closest to your end point. But this time, observe its "target" property value.


Once you've noted both ids, plug them into the form, click the "Calculate" button and pgRouting will do its magic!


So what does the "Calculate" button do exactly? It does this:

1. Creates a temporary MapGuide Feature Source (SDF) to store the shortest path result

2. Creates a MgLayer pointing to the shortest path result Feature Source

3. Use MgFeatureService.ExecuteSqlQuery() to execute the following SQL:

select gid, geom from roads 
            join (select * from shortest_path('select gid as id, source::integer, target::integer, length::double precision as cost from roads', source_id, target_id, false, false)) as route 
            on roads.gid = route.edge_id

The key thing to take from this, is that any PostGIS geometry columns in that SQL query will be correctly interpreted as FDO geometry values in the MgSqlDataReader we get back. If what I just said didn't ignite any light bulbs, well then it should!

Because this means for cases where our standard MapGuide Data Access APIs doesn't support or cover, we can always use MgFeatureService.ExecuteSqlQuery() as an "escape hatch" API to tap into the full set of geospatial/routing functionality offered by pgRouting/PostGIS and if any such SQL queries involve geometry columns, we can trust the PostgreSQL FDO provider to properly translate such values back as FDO geometry values in the returning SQL reader (and not as unrecognizable BLOBs), allowing us to do geometry operations on such results with the MapGuide API.

Does this technique only work for the PostgresSQL provider? No, the King.Oracle provider is also known to properly translate FDO geometry values from raw SQL query results via FDO. MySQL, SQL Server and SQLite are currently unknown. Please do comment if you know if this technique also works on these providers.

Now we can't create a Layer Definition off of a SQL query result (wouldn't that be a nice feature to have?), so we do the next best thing, which is to:

4. Copy the SQL reader contents into our temporary SDF feature source (that's why we created it). On subsequent calculations, the temporary SDF feature source is cleared first before copying in the new results from the SQL reader.

5. Save the MgMap and output a client-side JS call to refresh the map.

So that's the jist of this example. You can read the PHP code in the sample for all the gory details.

Caveats (or: pgRouting-isms)

EDIT: Thanks to dkastl for pointing out the very high tolerance value I had previously in my pgr_assign_vertex_id() example. A proper tolerance value will ensure proper display of routing results. Consider what was previously written below as a case of what happens if you specify a tolerance value that's too high.

If I haven't made it clear, this is a proof of concept example. There's some data/pgRouting related things that baffle me. For example, zoom in closer to that shortest path. Why is it broken off in some parts?


This shortest path result is straight from the pgRouting shortest_path() function.

Is the source data bad? Did I set up the network routing information incorrectly? Is my pgRouting SQL query bad? Did I base this example off of semi-inaccurate/outdated tutorial documentation? Is it a bug in pgRouting itself?

These are questions a more seasoned expert on pgRouting could probably explain. But if you do know why such results would be the case, I'd like to know through the comments on this post.

But besides that, the root question of whether MapGuide be integrated with pgRouting has been answered.

With an emphatic YES!

13 comments:

dkastl said...

Your issue with the returned route is probably a data (network topology) issue.

When you run select pgr_assign_vertex_id('roads', 0.001, 'geom', 'gid'); to calculate source and target attributes, then "0.001" is the "snapping tolerance" in the unit of your data.

You may have bad quality data where roads are not connected, and then making this parameter bigger is a possibility to fix errors. But if your data is good, then there is no need to snap nearby nodes. And even worse, when you chose the wrong unit, then you create a wrong network topology.

I think this could be the reason for your issues: try to run assign_vertex_id function again with something much smaller (ie. 0.000000001) as your projection unit is in degree, right?

Jackie Ng said...

Good point. I did copypasta the tolerance value from pgrouting tutorial that I assumed was relevant to this dataset. Maybe that's not the case, I'll have to double check.

dkastl said...

Well, I didn't look at the tutorial in detail, but if 0.0001 is for a projection in meter then this means a snapping tolerance that is extremely small. In case of degree as unit it's relatively big.

But if you can assume that your data source is of good quality, it's totally fine to make this snapping value very small.

Jackie Ng said...

Yes it was indeed the high tolerance value. Using a more fine-grained value fixes the problem.

This blog post has been updated to reflect these findings. Thanks for the heads up.

khoa kieu said...

Hi guys,
When I came to run "select pgr_assign_vertex_id('roads', 0.001, 'geom', 'gid')", I got an ERROR: function pgr_assign_vertex_id(unknown, numeric, unknown, unknown) does not exist.
Could u guys help me out this problem? Thank so much.
I'm using PostgresSql 9.3, postgis 2.1, pgrouting 2.0

iqnaul said...

Hi, I'm a new user how can I find the "Library://pg_routing/Data/postgis.FeatureSource" ? Because I looked at the MapGuide installation directory there is no such file, this made the data from PostgreSQL cannot be displayed.

Jackie Ng said...

You missed the link to the necessary setup instructions that I already put in the blog post.

Vikram Santhanam said...

I would like to perform the same in MSSQL 2012.. Is it possible ?

Jackie Ng said...

Unless there's a routing extension for SQL Server 2012, no.

Vikram Santhanam said...

Is there any routing extension?

Parth Trivedi said...

I am not able to find sample application for mapguide. Please provide download link

Amatya Gupta said...

Can their are options to implement Linear Referencing Through MapGuide????
What are the steps and their options?

Diego Ivan said...

We have successfully tested this approach with MapGuide and AIMS, it works very well. Thanks

Now the question is how to perform the routing analysis from arbitrary points within the edges instead of specific nodes ( start and end point of each edge).

Thank you

Diego