View here on youtube: https://www.youtube.com/watch?v=oOhbbEkl4Kg

PostGIS 2.0 is out and the awesomness continues! You can install PostGIS 2.0 on Ubuntu using packages which is exactly what I am going to show you here. Read on for details on how to get up and running and do your first simple raster analysis! Note: You should make good backups first!   Before we... Read more »

Having a good world dataset is kind of essential for most GIS users to have as a standard reference dataset, give context to maps and so on. Anita Graser’s recent blog post about the Natural Earth Data project added another nice dataset to my world dataset collection. The Natural Earth Data set (I got the... Read more »

Just a quickie because I always forget this and need to look it up again. If you have a query in postgresql and you want to include a sequential identifier in the output, you can use the very handy ROW_NUMBER() function. Note, you need Postgresql 8.4+ I believe. Here is a quick example: create view [...]

Here is a little batch script I wrote to backup my postgresql databases all in one go. I use similar scripts on my servers but this on is a bit different in that it has an outer loop to backup the three different instances of postgres I have running on my system, each in turn [...]

I love bash and the gnu tools provided on a Linux system. I am always writing little oneliners that help me do my work. Today I was generating some fixtures for django (yes I am finally properly learning to use unit testing in django). I wanted to know how many records where in each table in my database. Here is the little one liner I came up with:

for TABLE in $(psql foo-test -c "\dt" | awk '{print $3}'); do \
psql foo-test -c "select '$TABLE', count(*) from $TABLE;" | grep -A1 "\-\-\-\-" | tail -1; done

It produces output that looks something like this:

auth_group |     0
auth_group_permissions |     0
auth_message |     0
auth_permission |   273
auth_user |   366
auth_user_groups |     0
auth_user_user_permissions |     0
etc.

The PostgreSQL tablespace system is a convenient way to manage storage space allocation in PostgreSQL and perhaps improve performance. Tablespaces allow you to spread your databases or parts of your databases around different parts of your filesystem, and more importantly onto different physical storage devices. This is important because under ubuntu / debian, the default is to put the database cluster into the /var directory under the root of the filesystem heirarchy. Of course you can always mount the postgres subdirectory under var on a different filesystem, but that doesnt offer much flexibility. Here is a little example of how I spread my ‘gis’ postgis enabled database across two tablespaces – one located on my /home directory (which is /dev/sdb1) and one on my /opt directory (which is /dev/sda1). I’ll put database indexes in the /opt subdir and actual data in the /home subdirectory.

First we create two subdirectories and change their ownership to postgres:

sudo mkdir /opt/pg9-indexes
sudo mkdir /home/pg9-data
sudo chown postgres /opt/pg9-indexes
sudo chown postgres /home/pg9-data

Now let’s go ahead and create the tablespaces in postgres (I used template1, but you can be in any database). You need superuser rights for this:

psql template1
CREATE TABLESPACE optspace LOCATION '/opt/pg9-indexes';
CREATE TABLESPACE homespace LOCATION '/home/pg9-data';

Now let’s create a database, storing its data in /home/pg9-data:

createdb --tablespace=homespace gis

If you look in the pg9-data directory, you will see a new subdirectory has been made:

sudo ls /home/pg9-data/PG_9.0_201008051
17332

Then go ahead and load your data into postgres. Now lets create a spatial index for a table in my ‘world’ schema, setting the tablespace to the /opt partition:

CREATE INDEX idx_cities_geom ON world.cities \
USING GIST ( the_geom ) TABLESPACE optspace;

Now a quick look in our /opt directory should show you the index existing on the separate drive:

sudo ls /opt/pg9-indexes/PG_9.0_201008051
17332

If you want to automate the updating of all indexes for a given schema you can do it in a few steps like this (inspired by comments in this blog article):

psql gis
\o /tmp/updateindexes.sql
select 'ALTER INDEX world.'||indexname||' SET TABLESPACE optspace ;' \
from pg_catalog.pg_indexes where schemaname = 'world' order by tablename;
\q

In this case I am generating a little sql that will update all the indexes in my ‘world’ schema to use the optspace tablespace, and writing it out to a file (/tmp/updateindexes.sql). Next we can run  this on our database (from bash):

cat /tmp/updateindexes.sql | grep ALTER | psql gis

Thats it! One thing to note, if you are a windows user you probably can’t benefit from table spaces as the docs seem to indicate that you need to have a filesystem that supports symlinking….though I haven’t verified this one way or the other myself.

Ok so I’m feeling a bit left behind – everyone else seems to have jumped over to 9.0 and I have resolutely been sticking to the stable version shipped with ubuntu or debian (depending on the machine I am working on). However a client has a bunch of gis data in pg 9.0 / postgis 1.5 and I needed to be able to replicate their environment. In this article I will describe how I got Postgresql 8.4 running side by side with Postgresql 9.0, both with Postgis installed.

I will assume you already have postgresql 8.4 and postgis installed via Ubuntu apt…

First install postgres 9.0 from ppa

Add this to your /etc/apt/sources.list (as root):

# ppa for postgis 9.x
deb http://ppa.launchpad.net/pitti/postgresql/ubuntu natty main
deb-src http://ppa.launchpad.net/pitti/postgresql/ubuntu natty main

Then do:

sudo apt-get update

Note: Although pg9.1 beta packages are available, postgis 1.5  wont work with it so use 9.0

sudo apt-get install postgresql-contrib-9.0 postgresql-client-9.0 \
postgresql-server-dev-9.0 postgresql-9.0

You should check that you have the correct pg_config in your path

pg_config --version

If it is not 9.0, link it:

cd /usr/bin
sudo mv pg_config pg_config.old
sudo ln -s /usr/lib/postgresql/9.0/bin/pg_config .

Download and install postgis

There are no packages available for postgresql 9 so we need to build from source.

cd ~/dev/cpp
wget -c http://postgis.refractions.net/download/postgis-1.5.2.tar.gz
tar xfz postgis-1.5.2.tar.gz
cd postgis-1.5.2
./configure

You should get something like this (note building against postgresql 9)

PostGIS is now configured for x86_64-unknown-linux-gnu
 -------------- Compiler Info -------------  
 C compiler:           gcc -g -O2 
C++ compiler:         g++ -g -O2
 -------------- Dependencies --------------  
GEOS config:          /usr/bin/geos-config 
GEOS version:         3.2.0 
PostgreSQL config:    /usr/bin/pg_config 
PostgreSQL version:   PostgreSQL 9.0.4 
PROJ4 version:        47 
Libxml2 config:       /usr/bin/xml2-config 
Libxml2 version:      2.7.8 
PostGIS debug level:  0
 -------- Documentation Generation --------  
xsltproc:             /usr/bin/xsltproc 
xsl style sheets:      
dblatex:               
convert:              /usr/bin/convert

Now build and install

make
sudo make install

Using 8.4 and 9.0 side by side

Apt will install 9.0 side by side with 8.x without issue – it just sets it to run on a different port (5433 as opposed to the default 5432). This means you can comfortably use both (at the expense of some disk and processor load on your system). In order to make sure you are using the correct database, always add a -p 5433 when you wish to connect to 9.0. Your 8.4 should be running unchanged. The 9.0 install will be ‘blank’. to begin with (no user created tables). I always make myself a superuser on the database and create and account that matches my unix login account so that I can use ident authentication for my day to day work:

sudo su - postgres
createuser -p 5433 -d -E -i -l -P -r -s timlinux

Note the -p option to specify that I am connecting to postgresql 9.0! After you have entered your password twice type “`exit“` to exit the postgres user account.

Install postgis tables and functions into template1

You can create a new template, install postgis functions into each database that needs them separately, or (as I prefer), install them to template1 so that every database you creates will automatically be spatially enabled. I do this because I very seldom need a non-spatial database, and I like to have postgis there as default.

createlang -p 5433 plpgsql template1
psql -p 5433 template1 -f \
/usr/share/postgresql/9.0/contrib/postgis-1.5/postgis.sql
psql -p 5433 template1 -f \
/usr/share/postgresql/9.0/contrib/postgis-1.5/spatial_ref_sys.sql

Validating that postgis is installed

To do this we will connect to the template1 database and run a simple query:

psql -p 5433 template1
select postgis_lib_version();
\q

Should produce output like this:

WARNING: psql version 8.4, server version 9.0.         
Some psql features might not work.Type "help" for help.
template1=# select postgis_lib_version(); 
postgis_lib_version 
--------------------- 1.5.2(1 row)

Creating a new database

The last thing to do is to create a new database. As you may have noticed above, the 8.4 command line tools are still default as they have PATH preference, so when working explicitly with pg 9.0, its good to do:

export PATH=/usr/lib/postgresql/9.0/bin:$PATH

in your shell to put postgresql first in your path. A quick test will confirm that it has worked:

createdb --version
createdb (PostgreSQL) 9.0.4

With that setup nicely, lets create our new database:

createdb -p 5433 foo
psql -p 5433 -l

You can verify the database was created correctly by using psql -l (once again also explicitly setting the port no.).
Thats all there is to it! Now you can connect to the database from QGIS (just make sure to use the correct port number!), load spatial data into it using shp2pgsl and so on. I for one am looking forward to trying out some of  the new goodies in Postgresql 9.0.

Update: There are packages available for postgis (see comments below). Also, for convenience you can set your pg port to 5433 for your session by doing:

export PGPORT=5433

The GeoNames dataset provides a list of global placenames, their location and some additional information such as population and so on. It is published under the Creative Commons Attribution 3.0 License, which allows you to use the data for your own purposes. You can download the data by country, or the entire global dataset. In this article, I will walk you though how I downloaded the entire dataset, loaded it into PostgreSQL and added a geometry column so that I could view it in QGIS. Note that you can substitute these instructions for a specific country’s data easily.

First up, lets get the data from the geonames download page!

wget -c http://download.geonames.org/export/dump/allCountries.zip

Note the download is around 196mb so if you live in an internet backwater like I do, expect it to take a little while. If the download gets disconnected, just rerun the same command again – the ‘-c’ option tells wget to continue where it left off last time.

Once the data is downloaded, unzip it:

unzip allCountries.zip

You should now have a text file called allCountries.txt weighing in at just under 900mb. Next we can load it into PostgreSQL using a variation of this article. I highly recommend the use of schemas to partition your database into logical units. In the code listings that follow, it is assumed you have a schema called ‘world’. If you need to create it, simply do:

create schema world;

From the psql prompt. Since I am only interested in the geoname table at the moment I simply do this in my database.

create table world.geoname (
         geonameid       int,
         name            varchar(200),
         asciiname        varchar(200),
         alternatenames  varchar(8000),
         latitude        float,
         longitude       float,
         fclass  char(1),
         fcode   varchar(10),
         country varchar(2),
         cc2 varchar(60),
         admin1  varchar(20),
         admin2  varchar(80),
         admin3  varchar(20),
         admin4  varchar(20),
         population      bigint,
         elevation       int,
         gtopo30         int,
         timezone varchar(40),
         moddate         date
 );

You will notice that I extended the alternatenames field size from the original tutorial’s 4000 characters to 8000 characters in order to accommodate some longer entries that were causing my imports to fail with 4000 chars.

Next we can import the data (also from the psql prompt):

copy world.geoname (geonameid,name,asciiname,alternatenames,
latitude,longitude,fclass,fcode,country,cc2,
admin1,admin2,admin3,admin4,population,elevation,gtopo30,
timezone,moddate) from '/home/web/allCountries.txt' null as '';

Once again this is similar to the import line used by the original article I used, except I have used a full path to my allCountries.txt file. The import may take a little while depending on the speed of your computer.

When it is completed, you should have a bunch of data (~7.5 million records) in your table:

gis=# select count(*) from world.geoname ;
  count
---------
 7664712
(1 row)

It is going to be useful to have a unique id for every record – QGIS for one would like to have it, so lets add one (in addition to the geonameid field):

alter table world.geoname add column id serial not null;
CREATE UNIQUE INDEX idx_geoname_id ON world.geoname (id);

Because I know I will be using some other fields as basis for subset queries etc., I added some indexes on those too:

CREATE INDEX idx_geoname_population ON world.geoname (population);
CREATE INDEX idx_geoname_fcode ON world.geoname(fcode);

Ok thats all great, but there is no geometry column yet so we can’t view this in QGIS easily. So lets GIS enable the table! First we add a new geometry column:

alter table world.geoname add column the_geom geometry;

Now we populate the geometry column:

update world.geoname set the_geom = st_makepoint(longitude,latitude);

Next we add a constraint to ensure that the column always contains a point record or is null

alter table world.geoname add constraint enforce_geotype_the_geom
CHECK (geometrytype(the_geom) = 'POINT'::text OR the_geom IS NULL);

Finally lets index the table on the geometry field:

CREATE INDEX idx_geoname_the_geom ON world.geoname USING gist(the_geom);

Ok next we can connect to the database using QGIS and view our data! Note that you may want to filter the data or set up scale dependent visibility so that QGIS doesn’t try to render all 7.5 million points when you zoom out.

I added a query filter like this to the layer properties -> general tab -> subset:

"population" >= 10000 AND fcode != 'PPLA' and fcode != 'PCLI' AND fcode != 'ADM1'

You should experiment and read the geonames documentation in order to define a filter that is useful for your purposes.

Geonames layer loaded in QGIS

The Geonames dataset is a wonderful asset to the GIS community, particularly to those of us who value Free Software and Free Data.

Update 6 Aprl 2011:

I encountered one issue with the above procedure: the SRID for the imported points is -1 when loaded instead of 4326. This cause problems e.g. in mapserver you get an error like this:

ERROR:  Operation on two geometries with different SRIDs

To fix the issue I ran the following update query to assign all points a proper SRID:

update world.geoname set the_geom = ST_SETSRID(the_geom,4326);

Note that it takes quite a while to run.  When it completes, you can verify that the points are all in the correct CRS by running this little query:

gis=# select distinct ST_SRID(the_geom) from world.geoname;

Which should produce something like this:

 st_srid
---------
 4326
(1 row)

For good measure, I also added the geoname table to my geometry columns table:

insert into geometry_columns (f_table_catalog,f_table_schema,f_table_name,f_geometry_column,
coord_dimension,srid,type) values ('','world','geoname','the_geom',2,4326,'POINT');

Lastly I also gave select permissions to my readonly user (which I use when publishing in a web environment):

grant usage on schema world to readonly;
grant select on world to readonly;

I have also fixed the formatting of some of the SQL snippets above for those who struggled to read it within the width of this page.