Thursday, December 6, 2012

Importing OpenStreetMap data into Postgres

I have already set up PostGres 9.1 and PostGIS 2.0 installed. Create a new Postgres DB named GIS.
Run CREATE EXTENSION postgis on this database to add all of the required postgis functions and the spatial reference data.

 Run legacy.sql and rtpostgis_legacy.sql on your new GIS database to add in the legacy function names. OSM2PSQL has not been compiled recently, and still utilizes the old functions...

 You also have to recreate one function in the version I am using... supposedly this function is added back in the latest release... 2.1 I believe...
 CREATE OPERATOR CLASS gist_geometry_ops  
      FOR TYPE geometry USING GIST AS  
      STORAGE box2df,  
      OPERATOR    1    << ,  
      OPERATOR    2    &<      ,  
      OPERATOR    3    && ,  
      OPERATOR    4    &>      ,  
      OPERATOR    5    >>      ,  
      OPERATOR    6    ~=      ,  
      OPERATOR    7    ~      ,  
      OPERATOR    8    @      ,  
      OPERATOR    9    &<| ,  
      OPERATOR    10    <<| ,  
      OPERATOR    11    |>> ,  
      OPERATOR    12    |&> ,  
      OPERATOR    13    <-> FOR ORDER BY pg_catalog.float_ops,  
      OPERATOR    14    <#> FOR ORDER BY pg_catalog.float_ops,  
      FUNCTION    8    geometry_gist_distance_2d (internal, geometry, int4),  
      FUNCTION    1    geometry_gist_consistent_2d (internal, geometry, int4),  
      FUNCTION    2    geometry_gist_union_2d (bytea, internal),  
      FUNCTION    3    geometry_gist_compress_2d (internal),  
      FUNCTION    4    geometry_gist_decompress_2d (internal),  
      FUNCTION    5    geometry_gist_penalty_2d (internal, internal, internal),  
      FUNCTION    6    geometry_gist_picksplit_2d (internal, internal),  
      FUNCTION    7    geometry_gist_same_2d (geom1 geometry, geom2 geometry, internal);  

This is referenced in the following bug: http://trac.osgeo.org/postgis/ticket/1287 Then just run the osm2psql importer using the following command and you will be good.
 osm2pgsql -c -d gis -U postgres -W -H localhost -P 5432 -S default.style C:\gisdatafiles\2012-12-05_new-brunswick.osm  

Note that this imports the tables using spherical mercator projections... SRID=900913 There are command line switches that are supposed to allow you to change this, but I could not get them to work. For me this is just a temporary table anyway, so to add a translation field and then translate to another SRID is very easy to do so I did not bother wasting any more time trying to get that to work. (See my other blog post for instructions)

Tuesday, December 4, 2012

Where in the world? Adventures in GeoSpatial coding

I have recently been experimenting with GeoSpatial coding... there are some excellent and really cool technologies out there to help with the support of GeoSpatial coding.

I have been playing around with a little application to map as many geo-coded things as I can within my little corner of the world.  This is really just a geeky hobby, but has allowed me to touch on some REALLY cool technologies...

This blog is to capture some of the learnings I have had, for my future reference, and may even help others voyaging into this area.

My current task is to create a local mapping portlet showing all of the buildings that are mapped by the OpenStreetMaps project (openstreetmaps.org).  This dataset is a crowd sourced open geospatial dataset.  Depending on where you are in the world, there will be more or less accurate data available.  I am currently testing for Saint John, New Brunswick, Canada... and there is pretty decent data for this city.

The technologies I am using for this task are a PostGres 9.1 Database, installed with PostGIS extensions.
I have brought down the OpenStreetMap.org dataset for New Brunswick from the GEOFABRIK mirror site:  http://download.geofabrik.de/openstreetmap/north-america/canada/
And have imported it into my PostGres database using the OSM2PGSQL tools:
http://wiki.openstreetmap.org/wiki/Osm2pgsql

Note:  This did not work first try... as it turns out there have been some changes to the latest version of the PostGIS library which actually break this application since it has not been re-compiled against the new postGIS library.  There is good information on the PostGIS FAQ for this issue... basically you have to run a series of PostGIS .sql files which are included with the installation to make the functions backwards compatible... then do your import, and then revert the changes back with some other scripts....  there was one other problem I encountered too... which they did not have a script for, but for which there was a bug report on their site when I did a google search of the error... I copied the sql from the bug report 'fix' and ran it, and it worked after that.  I am sorry I do not remember specifically what this was, it was a couple of weeks ago when I ran this step... also, this utility may be updated by the time you do the import and/or the postGIS may have gone up a version so these instructions may be moot.  Either way, it was not that difficult to get running....

So at this point, I had a database called GIS which contained all of the data from the OpenMap dump for New Brunswick.  I have created my own application, with my own domain objects, so the next step was to bring the data over to my domain.  When you import the openmap data, all data is stored according to the type of geospatial data.   So for example, polygon data is stored in the planet_osm_polygon table.  In my domain I had created an 'open_street_map_building' table with corresponding JPA EJB Domain object.
To populate the data, it was just a matter of doing a simple PostGres dblink query (since my tables are in different databases:

 INSERT into open_map_building (id,amenity,leisure,name,om_id,shop,sport,tourism,way) 
     select nextval('building_seq'),amenity,leisure,name,om_id,shop,sport,tourism,way 
     from dblink('dbname=gis password=passwordhere', 
     'select amenity,leisure,name,osm_id,shop,sport,tourism,way from planet_osm_polygon where building <>''''') 
     as t1(amenity varchar(255),leisure varchar(255),name varchar(255),om_id bigint,shop varchar(255),sport varchar(255),tourism varchar(255),way geometry); 


The openstreetmap data is flattened out in the imported schema to have a variety of optional 'tags' inserted as columns.  This is fine, and as you can see I brought in a few of them that I was interested in for buildings...
I also imported the local business data into another table... basically POINT data and the name of the business...

Once this was brought in, I was able to load up inside of my EJBs as JPA objects.  I used the Hibernate Spatial library to provide the link between my Glassfish v3.1.2 container and my PostGres database.  By using the Hibernate Spatial library, you will be able to use the Spatial Datatypes in your db queries, and can map them in your JPA domain objects ... VERY COOL!  This was a pain in the ass to get working, I must admit.  You have to be very careful to put the hibernate spatial jars in the correct place or the EJB3 classloader will not be able to locate them... put them in the glassfish/domain/domain1/lib folder and make sure there are not other copies upstream (in the glassfish/lib).

my glassfish/domain/domain1/lib folder contains the following:
antlr-2.7.7.jar
commons-collections-3.1.jar
dom4j-1.6.1.jar
hibernate-commons-annotations-4.0.1.Final.jar
hibernate-core-4.0.0.Final.jar
hibernate-entitymanager-4.0.0.Final.jar
hibernate-spatial-4.0-M1.jar
javassist-3.15.0-GA.jar
jboss-logging-3.1.0.GA.jar
jboss-transaction-api_1.1_spec-1.0.0.Final.jar
jts-1.12.jar
log4j-1.2.17.jar

I then created a Liferay Portlet using Vaadin 6 and the OpenLayers add-in
https://vaadin.com/directory#addon/openlayers-wrapper
And map the OpenMap building onto it...


COOL.  I have a map dynamically displaying the business locations....
Now I want to retrieve a collection of buildings within the perimeter of the map being displayed...

Here I ran into a bit of a snag initially.  PostGIS provides some excellent spatial functions (and they seem to run really fast too...)

So I wrote an EJB method that would accept the lat/long of the bottom left and upper right corners of my map.  I added a column to my open_map_buildings table called 'centerpoint' which I populated by calling the postgis method UPDATE open_map_building SET centerpoint = ST_Centroid(way);
This worked great.  Then I attempted to so a query:

Select * from open_map_building where ST_Contains(ST_MakeBox2D(ST_Point(bottom, left),ST_Point(top,right)),centerpoint);
This however always returned no results... strange...

So I decided to look at the points in the database.  To make a human readable version of the GEOMETRY fields, just use the ST_AsText() function:  Select ST_AsText(centerpoint) FROM open_map_building.
I noticed that the coordinates looked strange compared to my lat/longs.  That is when I realized that the spatial_reference system being used were different.  I am not going to go into the gory details, do a google search to get lots of great information... but basically, there are a few standard spatial reference systems that are used.  Most GPS systems use WGS84 (srid=4326 in postGIS).... when the data is imported from the openstreetmap data, it is in Sperical Mercator format though... srid=3857 ... I believe...

Once I realized this, it was a fairly simple matter to create a new field and populate it with new coordinates translated into  WGS84.  To do this, I actually added two fields, one for the POLYGON containing the actual building perimeter, and one for the POINT containing the center of the building...

The postgis function to create the columns are:

 SELECT AddGeometryColumn('open_map_building','bldg_perimeter',4326,'POLYGON',2);

 SELECT AddGeometryColumn('open_map_building','bldg_center',4326,'POINT',2);


Then UPDATE was simply:

 UPDATE open_map_building SET bldg_perimeter = ST_Transform(way,4326); 

 UPDATE open_map_building SET bldg_center = ST_Transform(centerpoint,4326); 


Then to retrieve the buildings was a piece of cake:

SELECT id, name,ST_AsText(bldg_center) as way FROM open_map_building where ST_Contains(ST_SetSRID(ST_MakeBox2D(ST_Point(-66.06097291707987,45.27056897229072),ST_Point(-66.0572178244588,45.27283408390478)),4326),bldg_center);