PostGIS Manual

Some packaged distributions of PostGIS (in particular the Win32 installers for PostGIS >= 1.1.5) load the PostGIS functions into a template database called template_postgis. If the template_postgis database exists in your PostgreSQL installation then it is possible for users and/or applications to create spatially-enabled databases using a single command. Note that in both cases, the database user must have been granted the privilege to create new databases.

From the shell:

# createdb -T template_postgis my_spatial_db

From SQL:

postgres=# CREATE DATABASE my_spatial_db
        TEMPLATE=template_postgis

Upgrading existing spatial databases can be tricky as it requires replacement or introduction of new PostGIS object definitions.

Unfortunately not all definitions can be easily replaced in a live database, so sometimes your best bet is a dump/reload process.

PostGIS provides a SOFT UPGRADE procedure for minor or bugfix releases, and an HARD UPGRADE procedure for major releases.

Before attempting to upgrade postgis, it is always worth to backup your data. If you use the -Fc flag to pg_dump you will always be able to restore the dump with an HARD UPGRADE.

$ psql -f lwpostgis_upgrade.sql -d
          your_spatial_database

$ utils/postgis_proc_upgrade.pl lwpostgis.sql
            > lwpostgis_upgrade.sql

$ pg_dump -Fc olddb > olddb.dump

$ sh utils/postgis_restore.pl lwpostgis.sql newdb
          olddb.dump -E=UNICODE > restore.log

$ grep ^KEEPING restore.log | less

$ psql newdb -c "ALTER TABLE geometry_columns
          DROP attrelid" $ psql newdb -c "ALTER TABLE geometry_columns
          DROP varattnum" $ psql newdb -c "ALTER TABLE
          geometry_columns DROP stats"

$ psql newdb newdb=> delete from
          spatial_ref_sys; DROP newdb=> \i spatial_ref_sys.sql

There are several things to check when your installation or upgrade doesn't go as you expected.

Also check that you have made any necessary changes to the top of the Makefile.config. This includes:

The OpenGIS specification defines two standard ways of expressing spatial objects: the Well-Known Text (WKT) form and the Well-Known Binary (WKB) form. Both WKT and WKB include information about the type of the object and the coordinates which form the object.

Examples of the text representations (WKT) of the spatial objects of the features are as follows:

The OpenGIS specification also requires that the internal storage format of spatial objects include a spatial referencing system identifier (SRID). The SRID is required when creating spatial objects for insertion into the database.

Input/Output of these formats are available using the following interfaces:

bytea WKB = asBinary(geometry); text WKT =
        asText(geometry); geometry = GeomFromWKB(bytea WKB, SRID); geometry =
        GeometryFromText(text WKT, SRID);

For example, a valid insert statement to create and insert an OGC spatial object would be:

INSERT INTO SPATIALTABLE ( THE_GEOM, THE_NAME ) VALUES
        ( GeomFromText('POINT(-126.4 45.32)', 312), 'A Place'
        )

OGC formats only support 2d geometries, and the associated SRID is *never* embedded in the input/output representations.

PostGIS extended formats are currently superset of OGC one (every valid WKB/WKT is a valid EWKB/EWKT) but this might vary in the future, specifically if OGC comes out with a new format conflicting with our extensions. Thus you SHOULD NOT rely on this feature!

PostGIS EWKB/EWKT add 3dm,3dz,4d coordinates support and embedded SRID information.

Examples of the text representations (EWKT) of the extended spatial objects of the features are as follows:

Input/Output of these formats are available using the following interfaces:

bytea EWKB = asEWKB(geometry); text EWKT =
        asEWKT(geometry); geometry = GeomFromEWKB(bytea EWKB); geometry =
        GeomFromEWKT(text EWKT);

For example, a valid insert statement to create and insert a PostGIS spatial object would be:

INSERT INTO SPATIALTABLE ( THE_GEOM, THE_NAME ) VALUES
        ( GeomFromEWKT('SRID=312;POINTM(-126.4 45.32 15)'), 'A
        Place' )

The "canonical forms" of a PostgreSQL type are the representations you get with a simple query (without any function call) and the one which is guaranteed to be accepted with a simple insert, update or copy. For the postgis 'geometry' type these are:

 - Output - binary: EWKB ascii: HEXEWKB (EWKB in
        hex form) - Input - binary: EWKB ascii: HEXEWKB|EWKT 

For example this statement reads EWKT and returns HEXEWKB in the process of canonical ascii input/output:

=# SELECT 'SRID=4;POINT(0 0)'::geometry;
        geometry ----------------------------------------------------
        01010000200400000000000000000000000000000000000000 (1 row)

The SQL Multimedia Applications Spatial specification extends the simple features for SQL spec by defining a number of circularly interpolated curves.

The SQL-MM definitions include 3dm, 3dz and 4d coordinates, but do not allow the embedding of SRID information.

The well-known text extensions are not yet fully supported. Examples of some simple curved geometries are shown below:

The SPATIAL_REF_SYS table definition is as follows:

CREATE TABLE SPATIAL_REF_SYS ( SRID INTEGER NOT NULL
        PRIMARY KEY, AUTH_NAME VARCHAR(256), AUTH_SRID INTEGER, SRTEXT
        VARCHAR(2048), PROJ4TEXT VARCHAR(2048) )

The SPATIAL_REF_SYS columns are as follows:

The GEOMETRY_COLUMNS table definition is as follows:

CREATE TABLE GEOMETRY_COLUMNS ( F_TABLE_CATALOG
        VARCHAR(256) NOT NULL, F_TABLE_SCHEMA VARCHAR(256) NOT NULL,
        F_TABLE_NAME VARCHAR(256) NOT NULL, F_GEOMETRY_COLUMN VARCHAR(256) NOT
        NULL, COORD_DIMENSION INTEGER NOT NULL, SRID INTEGER NOT NULL, TYPE
        VARCHAR(30) NOT NULL )

The columns are as follows:

Creating a table with spatial data is done in two stages:

Here is an example of SQL used to create a table and add a spatial column (assuming that an SRID of 128 exists already):

CREATE TABLE parks ( PARK_ID int4, PARK_NAME
        varchar(128), PARK_DATE date, PARK_TYPE varchar(2) ); SELECT
        AddGeometryColumn('parks', 'park_geom', 128,
        'MULTIPOLYGON', 2 );

Here is another example, using the generic "geometry" type and the undefined SRID value of -1:

CREATE TABLE roads ( ROAD_ID int4, ROAD_NAME
        varchar(128) ); SELECT AddGeometryColumn( 'roads',
        'roads_geom', -1, 'GEOMETRY', 3 );

Most of the functions implemented by the GEOS library rely on the assumption that your geometries are valid as specified by the OpenGIS Simple Feature Specification. To check validity of geometries you can use the IsValid() function:

gisdb=# select isvalid('LINESTRING(0 0, 1
        1)'), isvalid('LINESTRING(0 0,0 0)'); isvalid | isvalid
        ---------+--------- t | f

By default, PostGIS does not apply this validity check on geometry input, because testing for validity needs lots of CPU time for complex geometries, especially polygons. If you do not trust your data sources, you can manually enforce such a check to your tables by adding a check constraint:

ALTER TABLE mytable ADD CONSTRAINT
        geometry_valid_check CHECK (isvalid(the_geom));

If you encounter any strange error messages such as "GEOS Intersection() threw an error!" or "JTS Intersection() threw an error!" when calling PostGIS functions with valid input geometries, you likely found an error in either PostGIS or one of the libraries it uses, and you should contact the PostGIS developers. The same is true if a PostGIS function returns an invalid geometry for valid input.

If you can convert your data to a text representation, then using formatted SQL might be the easiest way to get your data into PostGIS. As with Oracle and other SQL databases, data can be bulk loaded by piping a large text file full of SQL "INSERT" statements into the SQL terminal monitor.

A data upload file (roads.sql for example) might look like this:

BEGIN; INSERT INTO ROADS_GEOM (ID,GEOM,NAME ) VALUES
        (1,GeomFromText('LINESTRING(191232 243118,191108
        243242)',-1),'Jeff Rd'); INSERT INTO ROADS_GEOM
        (ID,GEOM,NAME ) VALUES (2,GeomFromText('LINESTRING(189141
        244158,189265 244817)',-1),'Geordie Rd'); INSERT INTO
        ROADS_GEOM (ID,GEOM,NAME ) VALUES
        (3,GeomFromText('LINESTRING(192783 228138,192612
        229814)',-1),'Paul St'); INSERT INTO ROADS_GEOM
        (ID,GEOM,NAME ) VALUES (4,GeomFromText('LINESTRING(189412
        252431,189631 259122)',-1),'Graeme Ave'); INSERT INTO
        ROADS_GEOM (ID,GEOM,NAME ) VALUES
        (5,GeomFromText('LINESTRING(190131 224148,190871
        228134)',-1),'Phil Tce'); INSERT INTO ROADS_GEOM
        (ID,GEOM,NAME ) VALUES (6,GeomFromText('LINESTRING(198231
        263418,198213 268322)',-1),'Dave Cres'); COMMIT;

The data file can be piped into PostgreSQL very easily using the "psql" SQL terminal monitor:

psql -d [database] -f roads.sql

The shp2pgsql data loader converts ESRI Shape files into SQL suitable for insertion into a PostGIS/PostgreSQL database. The loader has several operating modes distinguished by command line flags:

Note that -a, -c, -d and -p are mutually exclusive.

An example session using the loader to create an input file and uploading it might look like this:

# shp2pgsql shaperoads myschema.roadstable >
        roads.sql # psql -d roadsdb -f roads.sql

A conversion and upload can be done all in one step using UNIX pipes:

# shp2pgsql shaperoads myschema.roadstable | psql -d
        roadsdb

The most straightforward means of pulling data out of the database is to use a SQL select query and dump the resulting columns into a parsable text file:

db=# SELECT id, AsText(geom) AS geom, name FROM
        ROADS_GEOM; id | geom | name
        ---+-----------------------------------------+----------- 1 |
        LINESTRING(191232 243118,191108 243242) | Jeff Rd 2 |
        LINESTRING(189141 244158,189265 244817) | Geordie Rd 3 |
        LINESTRING(192783 228138,192612 229814) | Paul St 4 |
        LINESTRING(189412 252431,189631 259122) | Graeme Ave 5 |
        LINESTRING(190131 224148,190871 228134) | Phil Tce 6 |
        LINESTRING(198231 263418,198213 268322) | Dave Cres 7 |
        LINESTRING(218421 284121,224123 241231) | Chris Way (6 rows)

However, there will be times when some kind of restriction is necessary to cut down the number of fields returned. In the case of attribute-based restrictions, just use the same SQL syntax as normal with a non-spatial table. In the case of spatial restrictions, the following operators are available/useful:

Next, you can use these operators in queries. Note that when specifying geometries and boxes on the SQL command line, you must explicitly turn the string representations into geometries by using the "GeomFromText()" function. So, for example:

SELECT ID, NAME FROM ROADS_GEOM WHERE GEOM ~=
        GeomFromText('LINESTRING(191232 243118,191108 243242)',-1);

The above query would return the single record from the "ROADS_GEOM" table in which the geometry was equal to that value.

When using the "&&" operator, you can specify either a BOX3D as the comparison feature or a GEOMETRY. When you specify a GEOMETRY, however, its bounding box will be used for the comparison.

SELECT ID, NAME FROM ROADS_GEOM WHERE GEOM &&
        GeomFromText('POLYGON((191232 243117,191232 243119,191234
        243117,191232 243117))',-1);

The above query will use the bounding box of the polygon for comparison purposes.

The most common spatial query will probably be a "frame-based" query, used by client software, like data browsers and web mappers, to grab a "map frame" worth of data for display. Using a "BOX3D" object for the frame, such a query looks like this:

SELECT AsText(GEOM) AS GEOM FROM ROADS_GEOM WHERE GEOM
        && SetSRID('BOX3D(191232 243117,191232
        243119)'::box3d,-1);

Note the use of the SRID, to specify the projection of the BOX3D. The value -1 is used to indicate no specified SRID.

The pgsql2shp table dumper connects directly to the database and converts a table (possibly defined by a query) into a shape file. The basic syntax is:

pgsql2shp [<options>] <database>
        [<schema>.]<table>

pgsql2shp [<options>] <database>
        <query>

The commandline options are:

GiST stands for "Generalized Search Tree" and is a generic form of indexing. In addition to GIS indexing, GiST is used to speed up searches on all kinds of irregular data structures (integer arrays, spectral data, etc) which are not amenable to normal B-Tree indexing.

Once a GIS data table exceeds a few thousand rows, you will want to build an index to speed up spatial searches of the data (unless all your searches are based on attributes, in which case you'll want to build a normal index on the attribute fields).

The syntax for building a GiST index on a "geometry" column is as follows:

CREATE INDEX [indexname] ON [tablename] USING
        GIST ( [geometryfield] GIST_GEOMETRY_OPS ); 

Building a spatial index is a computationally intensive exercise: on tables of around 1 million rows, on a 300MHz Solaris machine, we have found building a GiST index takes about 1 hour. After building an index, it is important to force PostgreSQL to collect table statistics, which are used to optimize query plans:

VACUUM ANALYZE [table_name] [column_name]; --
        This is only needed for PostgreSQL 7.4 installations and below SELECT
        UPDATE_GEOMETRY_STATS([table_name], [column_name]);

GiST indexes have two advantages over R-Tree indexes in PostgreSQL. Firstly, GiST indexes are "null safe", meaning they can index columns which include null values. Secondly, GiST indexes support the concept of "lossiness" which is important when dealing with GIS objects larger than the PostgreSQL 8K page size. Lossiness allows PostgreSQL to store only the "important" part of an object in an index -- in the case of GIS objects, just the bounding box. GIS objects larger than 8K will cause R-Tree indexes to fail in the process of being built.

Ordinarily, indexes invisibly speed up data access: once the index is built, the query planner transparently decides when to use index information to speed up a query plan. Unfortunately, the PostgreSQL query planner does not optimize the use of GiST indexes well, so sometimes searches which should use a spatial index instead default to a sequence scan of the whole table.

If you find your spatial indexes are not being used (or your attribute indexes, for that matter) there are a couple things you can do:

When constructing a query it is important to remember that only the bounding-box-based operators such as && can take advantage of the GiST spatial index. Functions such as distance() cannot use the index to optimize their operation. For example, the following query would be quite slow on a large table:

SELECT the_geom FROM geom_table WHERE distance(
        the_geom, GeomFromText( 'POINT(100000 200000)', -1 ) ) <
        100

This query is selecting all the geometries in geom_table which are within 100 units of the point (100000, 200000). It will be slow because it is calculating the distance between each point in the table and our specified point, ie. one distance() calculation for each row in the table. We can avoid this by using the && operator to reduce the number of distance calculations required:

SELECT the_geom FROM geom_table WHERE the_geom
        && 'BOX3D(90900 190900, 100100 200100)'::box3d AND
        distance( the_geom, GeomFromText( 'POINT(100000 200000)', -1 )
        ) < 100

This query selects the same geometries, but it does it in a more efficient way. Assuming there is a GiST index on the_geom, the query planner will recognize that it can use the index to reduce the number of rows before calculating the result of the distance() function. Notice that the BOX3D geometry which is used in the && operation is a 200 unit square box centered on the original point - this is our "query box". The && operator uses the index to quickly reduce the result set down to only those geometries which have bounding boxes that overlap the "query box". Assuming that our query box is much smaller than the extents of the entire geometry table, this will drastically reduce the number of distance calculations that need to be done.

The examples in this section will make use of two tables, a table of linear roads, and a table of polygonal municipality boundaries. The table definitions for the bc_roads table is:

Column | Type | Description
        ------------+-------------------+------------------- gid | integer |
        Unique ID name | character varying | Road Name the_geom | geometry |
        Location Geometry (Linestring)

The table definition for the bc_municipality table is:

Column | Type | Description
        -----------+-------------------+------------------- gid | integer |
        Unique ID code | integer | Unique ID name | character varying | City /
        Town Name the_geom | geometry | Location Geometry (Polygon)

postgis=# SELECT sum(length(the_geom))/1000 AS
              km_roads FROM bc_roads; km_roads ------------------
              70842.1243039643 (1 row)

postgis=# SELECT area(the_geom)/10000 AS
              hectares FROM bc_municipality WHERE name = 'PRINCE
              GEORGE'; hectares ------------------ 32657.9103824927 (1
              row)

postgis=# SELECT name, area(the_geom)/10000 AS
              hectares FROM bc_municipality ORDER BY hectares DESC LIMIT 1;
              name | hectares ---------------+----------------- TUMBLER RIDGE
              | 155020.02556131 (1 row)

postgis=# SELECT m.name,
              sum(length(r.the_geom))/1000 as roads_km FROM bc_roads AS
              r,bc_municipality AS m WHERE r.the_geom && m.the_geom
              AND contains(m.the_geom,r.the_geom) GROUP BY m.name ORDER BY
              roads_km; name | roads_km
              ----------------------------+------------------ SURREY |
              1539.47553551242 VANCOUVER | 1450.33093486576 LANGLEY DISTRICT |
              833.793392535662 BURNABY | 773.769091404338 PRINCE GEORGE |
              694.37554369147 ...

postgis=# CREATE TABLE pg_roads as SELECT
              intersection(r.the_geom, m.the_geom) AS intersection_geom,
              length(r.the_geom) AS rd_orig_length, r.* FROM bc_roads AS r,
              bc_municipality AS m WHERE r.the_geom && m.the_geom AND
              intersects(r.the_geom, m.the_geom) AND m.name = 'PRINCE
              GEORGE';

postgis=# SELECT sum(length(r.the_geom))/1000 AS
              kilometers FROM bc_roads r, bc_municipality m WHERE r.the_geom
              && m.the_geom AND r.name = 'Douglas St' AND
              m.name = 'VICTORIA'; kilometers ------------------
              4.89151904172838 (1 row)

postgis=# SELECT gid, name, area(the_geom) AS
              area FROM bc_municipality WHERE nrings(the_geom) > 1 ORDER
              BY area DESC LIMIT 1; gid | name | area
              -----+--------------+------------------ 12 | SPALLUMCHEEN |
              257374619.430216 (1 row)

To use PostGIS with Mapserver, you will need to know about how to configure Mapserver, which is beyond the scope of this documentation. This section will cover specific PostGIS issues and configuration details.

To use PostGIS with Mapserver, you will need:

Mapserver accesses PostGIS/PostgreSQL data like any other PostgreSQL client -- using libpq. This means that Mapserver can be installed on any machine with network access to the PostGIS server, as long as the system has the libpq PostgreSQL client libraries.

The USING pseudo-SQL clause is used to add some information to help mapserver understand the results of more complex queries. More specifically, when either a view or a subselect is used as the source table (the thing to the right of "FROM" in a DATA definition) it is more difficult for mapserver to automatically determine a unique identifier for each row and also the SRID for the table. The USING clause can provide mapserver with these two pieces of information as follows:

DATA "the_geom FROM (SELECT table1.the_geom AS
        the_geom, table1.oid AS oid, table2.data AS data FROM table1 LEFT JOIN
        table2 ON table1.id = table2.id) AS new_table USING UNIQUE oid USING
        SRID=-1"

Lets start with a simple example and work our way up. Consider the following Mapserver layer definition:

LAYER CONNECTIONTYPE postgis NAME "roads"
        CONNECTION "user=theuser password=thepass dbname=thedb
        host=theserver" DATA "the_geom FROM roads" STATUS ON TYPE
        LINE CLASS COLOR 0 0 0 END END

This layer will display all the road geometries in the roads table as black lines.

Now lets say we want to show only the highways until we get zoomed in to at least a 1:100000 scale - the next two layers will achieve this effect:

LAYER CONNECTION "user=theuser password=thepass
        dbname=thedb host=theserver" DATA "the_geom FROM roads"
        MINSCALE 100000 STATUS ON TYPE LINE FILTER "road_type =
        'highway'" CLASS COLOR 0 0 0 END END LAYER CONNECTION
        "user=theuser password=thepass dbname=thedb host=theserver"
        DATA "the_geom FROM roads" MAXSCALE 100000 STATUS ON TYPE LINE
        CLASSITEM road_type CLASS EXPRESSION "highway" SIZE 2 COLOR
        255 0 0 END CLASS COLOR 0 0 0 END END

The first layer is used when the scale is greater than 1:100000, and displays only the roads of type "highway" as black lines. The FILTER option causes only roads of type "highway" to be displayed.

The second layer is used when the scale is less than 1:100000, and will display highways as double-thick red lines, and other roads as regular black lines.

So, we have done a couple of interesting things using only mapserver functionality, but our DATA SQL statement has remained simple. Suppose that the name of the road is stored in another table (for whatever reason) and we need to do a join to get it and label our roads.

LAYER CONNECTION "user=theuser password=thepass
        dbname=thedb host=theserver" DATA "the_geom FROM (SELECT
        roads.oid AS oid, roads.the_geom AS the_geom, road_names.name as name
        FROM roads LEFT JOIN road_names ON roads.road_name_id =
        road_names.road_name_id) AS named_roads USING UNIQUE oid USING
        SRID=-1" MAXSCALE 20000 STATUS ON TYPE ANNOTATION LABELITEM name
        CLASS LABEL ANGLE auto SIZE 8 COLOR 0 192 0 TYPE truetype FONT arial
        END END END

This annotation layer adds green labels to all the roads when the scale gets down to 1:20000 or less. It also demonstrates how to use an SQL join in a DATA definition.

...

...

Current PostgreSQL versions (including 8.0) suffer from a query optimizer weakness regarding TOAST tables. TOAST tables are a kind of "extension room" used to store large (in the sense of data size) values that do not fit into normal data pages (like long texts, images or complex geometries with lots of vertices), see http://www.postgresql.org/docs/8.0/static/storage-toast.html for more information).

The problem appears if you happen to have a table with rather large geometries, but not too much rows of them (like a table containing the boundaries of all European countries in high resolution). Then the table itself is small, but it uses lots of TOAST space. In our example case, the table itself had about 80 rows and used only 3 data pages, but the TOAST table used 8225 pages.

Now issue a query where you use the geometry operator && to search for a bounding box that matches only very few of those rows. Now the query optimizer sees that the table has only 3 pages and 80 rows. He estimates that a sequential scan on such a small table is much faster than using an index. And so he decides to ignore the GIST index. Usually, this estimation is correct. But in our case, the && operator has to fetch every geometry from disk to compare the bounding boxes, thus reading all TOAST pages, too.

To see whether your suffer from this bug, use the "EXPLAIN ANALYZE" postgresql command. For more information and the technical details, you can read the thread on the postgres performance mailing list: http://archives.postgresql.org/pgsql-performance/2005-02/msg00030.php

The PostgreSQL people are trying to solve this issue by making the query estimation TOAST-aware. For now, here are two workarounds:

The first workaround is to force the query planner to use the index. Send "SET enable_seqscan TO off;" to the server before issuing the query. This basically forces the query planner to avoid sequential scans whenever possible. So it uses the GIST index as usual. But this flag has to be set on every connection, and it causes the query planner to make misestimations in other cases, so you should "SET enable_seqscan TO on;" after the query.

The second workaround is to make the sequential scan as fast as the query planner thinks. This can be achieved by creating an additional column that "caches" the bbox, and matching against this. In our example, the commands are like:

SELECT
        addGeometryColumn('myschema','mytable','bbox','4326','GEOMETRY','2');
        UPDATE mytable set bbox = Envelope(Force_2d(the_geom));

Now change your query to use the && operator against bbox instead of geom_column, like:

SELECT geom_column FROM mytable WHERE bbox &&
        SetSrid('BOX3D(0 0,1 1)'::box3d,4326);

Of course, if you change or add rows to mytable, you have to keep the bbox "in sync". The most transparent way to do this would be triggers, but you also can modify your application to keep the bbox column current or run the UPDATE query above after every modification.