HOME

GRASS 5.7 Usage Examples


Tutorial HOME | Table of contents
README
GRASS 5.7 is currently under development. In case the examples described here do not work properly, you are kindly invited to send us further examples and/or code bugfixes/enhancements.

Example: Writing the GRASS 5.0 Spearfish 'soils' map into PostGIS with MapServer access
In this example we want to take the Spearfish 'soils' vector map and write it to PostGIS. This map will be also accessible by MapServer then.

First we create a new PostGIS database:

#you need the file postgis.sql from PostGIS in the current directory:
pg.postgisdb -c spearfish -h pgserver.itc.it -U postgres -G ./

#You may also create it directly in PostgreSQL.

We are ready to import the data now. As we are using the GRASS 5.0 spearfish data base, we have to convert the data to the new format (alternativly you can grab the Spearfish5.7 data set and skip v.convert):
v.convert input=soils output=soils

Note that the map is topologically correct in the new Spearfish5.7 data set.

The vector geometry and the attributes table are linked automatically. To verify that, we run:

v.info soils
 [...]
 Map format: native
 [...]
 Number of dblinks:      1
 [...]
 
#see history:
v.info -h soils

Now we'll write out the map to PostGIS:
#Export area from GRASS vector map directly to PostGIS:
v.out.ogr input=soils type=area dsn="PG:host=localhost dbname=postgis user=postgres" layer=soils format=PostgreSQL

Next we'll link it virtually back to GRASS:
v.external dsn="PG:host=localhost user=postgres dbname=postgis" layer=soils output=soilsPG

Let us look at the map now:
d.mon x0
g.region vect=soilsPG
d.vect soilsPG col=blue
v.db.connect -p soilsPG
d.what.vect soilsPG

To see the new tables in PostGRASS, we run:
db.connect driver=pg database="host=pgserver.itc.it,dbname=spearfish,user=postgres"
db.tables -p
db.describe -c soils

MapServer preparations

To make MapServer accept the PostGIS data, we need to fill the table 'GEOMETRY_COLUMNS' which links the geometry to the attributes within PostgreSQL.
One row has to be added:

db.connect driver=pg database="host=pgserver.itc.it,dbname=spearfish,user=postgres"
db.tables -p
echo "INSERT INTO geometry_columns VALUES ('','spearfish','streams_geom','geom','2','-1','GEOMETRY');" | db.execute
#verify:
echo "select * from geometry_columns" | db.select

Finally MapServer requires a mapfile entry similar to this (note: polygons are not supported yet in PostGRASS/MapServer and draw as boundaries only):

LAYER
  NAME "streams_geom"
  CONNECTIONTYPE postgis
  CONNECTION "user=postgres dbname=spearfish host=127.0.0.1"
  #      Column    geom_table
  DATA "geom from streams_geom"
  TYPE LINE
  STATUS DEFAULT
  CLASS
    SYMBOL "pattern1"
    NAME "label"
    OUTLINECOLOR 100 50 0
    COLOR 255 100 100
    LABEL
      POSITION CC
      SIZE SMALL
      COLOR 255 0 0
    END
  END
END
This 'streams' map does not contain attributes because of the related Spearfish map.

And here you can see how it looks like:
Demo Online MapServer reading GRASS raster, SHAPE, GeoTIFF and PostGRASS/PostGIS data

EXTRA: Add a table 'view'

Since we still have the geometry and attributes in different tables, we can add a 'view' onto both tables (sort of virtual table merging both). Note that for a 'view' MapServer 3.7 is required. If you use a version < 3.7 you should use 'CREATE TABLE' instead:

CREATE VIEW test AS select streams_geom.id, streams_cat.cat, streams_geom.type, streams_geom.geom \
 from streams_geom,streams_cat WHERE streams_geom.id = streams_cat.id AND \
 streams_cat.field = 1;
INSERT into geometry_columns values ( ' ', ' ', 'test2', 'geom', 2, -1, 'GEOMETRY' );
Then change the table name in the 'mapfile' (DATA entry) to the view. Should work... As the 'streams' map does not have exciting attributes, you may calculate the length of each vector, assign it to the attributes table and enjoy the map queries.

Adding stream length labels

Since the Spearfish 'streams' map does not contain attributes, we show how to add stream length labels. To have useful stream lengths, we first polish the map and generate polylines from lines (changing nodes along a vector to vertices, whereever possible):

v.build.polylines in=streams_catsnew out=streams_pl
d.vect streams_pl col=red

Now the existing 'cat' values must be removed as we need the vector lines numbered individually:
v.category input=streams_pl option=del out=streams_nocat
v.category input=streams_nocat option=report

The second command should report nothing. Now we add individual category values to each vector (increasing number):
v.category input=streams_nocat option=add step=1 cat=1 out=streams_length
v.category input=streams_length option=report
Now categories should be displayed.

Vector geometry statistics such as line length can be calculated with 'v.to.db' (print only first):

v.to.db -p streams_length option=length units=me col1=cat

The length values we want to assign in a new attribute column. Up to now there is no attribute table for the new map. We generate it therefore and add the length attributes. The link between vectors and attribute table rows is established by the 'cat' values. Note that uploading attributes of a vector map to the DBMS has to be done column-wise:
echo "create table streams_length (cat int, length float)" | db.execute
db.describe -c streams_length
v.db.connect -p streams_length

#set DBMS link for new map to new table:
v.db.connect driver=pg database="host=pgserver.itc.it,dbname=spearfish,user=postgres" table=streams_length
v.to.db streams_length option=cat col1=cat
v.to.db streams_length option=length col1=length

To verify, we can query the table and look at the map in the GRASS monitor:
echo "select * from streams_length" | db.select
d.vect streams_length col=brown display=shape,attr attrcol=length
To use the new map 'streams_length' also in MapServer, in table 'geometry_columns' a new row must be added:
echo "INSERT INTO geometry_columns VALUES ('','spearfish','streams_length_geom','geom','2','-1','GEOMETRY');" | db.execute
And finally the map has to be added to the 'mapfile' of MapServer.

Example: Conversion of SHAPE map to GRASS 5.7 vector map
This time we are not talking about simply registering a SHAPE map, but about really converting it into another spatial format.
WARNING: don't try this at home. Still pretty unstable.

v.out.ascii markveggy.shp out=test1
v.in.ascii test1 out=newmap
v.info newmap

#NOTE: this cannot work properly as SHAPE contains closed polygons
#v.clean should be used etc...

v.clean newmap out=newmap_clean tool=break,rmdupl
d.vect newmap_clean ty=area,boundary col=red

#we connect this map to the old DBF file in the 'dbln' file:
v.db.connect newmap_clean dr=dbf database=/ssi0/ssi/neteler/grassdata/botswanaLL/PERMANENT/shp table=markveggy key=VEGCNP_ID

#check if it is connected:
v.db.connect -p newmap_clean

#draw & query
d.erase
d.vect newmap_clean
d.what.vect newmap_clean
TODO: check if meanwhile g.copy does the same, probably followed by v.clean with rmdupl.


© 2002-2003 Markus Neteler
Comments about this page | FAQ | Download | Support | Docs | Programming | Back 5.7 Tutorial Home
Last change: $Date: 2008-03-27 21:31:14 +0000 (Thu, 27 Mar 2008) $