These are my notes on using the OSM data as imported by osm2pgsql in a PostgreSQL database with Grass GIS.
While Grass is a very powerful tool, it's Unixy feel is a little disconcerting. If you're a programmer like me, it helps to realize the Grass command line shell is in fact a preconfigured bash session:
Here's a map of pharmacies in my neighborhood as rendered by QGIS' export to PNG feature:
While Grass is a very powerful tool, it's Unixy feel is a little disconcerting. If you're a programmer like me, it helps to realize the Grass command line shell is in fact a preconfigured bash session:
GRASS 6.4.1 (OSM_900913):~ > ps ax | grep "$$" | grep -v grep 98804 s000 S 0:00.06 /bin/bashThe things you get from Grass are a set of binaries that interact with the data you have (organized in locations that are part of a mapset stored inside your filesystem). The other thing that is somewhat unexpected for a PostGIS user is that, within a Grass location, all the data has to share a single projection.
Creating a location for your OSM data
The first step required when using OSM data is to create a Grass location for the OSM data. I've used the same projection osm2pgsql uses when importing, the "fake" projection with SRID 900913. Grass doesn't know it, so you can define it as a custom PROJ4 string based on the definition you get from PostGIS:
osm=# select proj4text from spatial_ref_sys where srid = 900913;
proj4text
------------------------------------------------------------------------------------------------------------------
+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +units=m +k=1.0 +nadgrids=@null +no_defs
(1 row)
Here are the relevant screen shots when defining the location using Grass Python WXGUI:
Creating temporary tables in OSM
To reduce the size of the data set, I've created a couple of tables in PostGIS by selecting only those features I planned to use for my map by spatially selecting those features within a certain administrative boundary:create table bucharest_s1_streets as select * from planet_osm_line where st_contains((select way from planet_osm_polygon where boundary = 'administrative' and name ='Sector 1' limit 1), way) and highway is not null; create table bucharest_s1_pois as select * from planet_osm_point where st_contains((select way from planet_osm_polygon where boundary = 'administrative' and name ='Sector 1' limit 1), way);I then defined these tables in PostGIS' geometry_tables (otherwise Grass will refuse to import them):
insert into geometry_columns(f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, type) values('','public','bucharest_s1_pois','way',2,900913,'POINT');
insert into geometry_columns(f_table_catalog, f_table_schema, f_table_name, f_geometry_column, coord_dimension, srid, type) values('','public','bucharest_s1_streets','way',2,900913,'LINESTRING');
Importing the data in Grass
For some reason, my Grass setup fails to set the environment variables below, so I have to re-export them:GRASS 6.4.1 (OSM_900913):~ > export GISDBASE=/Users/diciu/grassdata GRASS 6.4.1 (OSM_900913):~ > export LOCATION_NAME=OSM_900913 GRASS 6.4.1 (OSM_900913):~ > export MAPSET=PERMANENTThe PostGIS database parameters have to be defined next:
export HOST=localhost
export DB=osm
export USER=gisuser
export PASS=somepassword
export DSN="PG:host=${HOST} dbname=${DB} user=${USER} password=${PASS}"
You can now run the v.ogr.in commands to import your data (I'm importing roads and pharmacies):
v.in.ogr dsn="${DSN}" layer=bucharest_s1_streets out=streets
v.in.ogr dsn="${DSN}" layer=bucharest_s1_pois out=pharmacy where="amenity='pharmacy'"
If everything went well, not you have two vector layers in your location:
GRASS 6.4.1 (OSM_900913):~ > g.list type=vect ---------------------------------------------- vector files available in mapset: pharmacy streets ----------------------------------------------
Computing distance to pharmacies
The Vector networking tutorial describes a number of very nice examples of processing tools available in Grass. I am using v.net.iso below to compute the distance from pharmacies:v.distance -p from=pharmacy to=streets output=pharmacy_conn_streets upload=dist column=dist v.patch in=pharmacy,streets,pharmacy_conn_streets out=pharmacy_net v.clean in=pharmacy_net out=pharmacy_net_clean tool=snap,break thresh=1 v.net.iso input=pharmacy_net_clean output=pharmacy_iso ccats=1-100 costs=200,500,1000 nlayer=1The command's output are a number of vector layers, pharmacy_iso being the one that stores our road network categorized by distance to pharmacies (less than 200 meters, between 200 and 500 meters and so on).
Rendering a map using the data
It's much easier and the results are probably better if you use QGIS instead of Grass' UI. There's a plugin for QGIS that asks you to specify the Grass mapset, location and loads a given layer. As an added bonus, QGIS is able to re-project on the fly so you can mix layers in different projections.
