Related Plugins and Tags

QGIS Planet

Installing PostGIS 2.0 on ubuntu

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 begin, you should uninstall your existing postgis packages:

sudo dpkg --purge postgis postgresql-9.1-postgis

Then add a new repository and install PostGIS from there (based on this post):

sudo apt-add-repository ppa:sharpie/for-science  # To get GEOS 3.3.2
sudo apt-add-repository ppa:sharpie/postgis-nightly
sudo apt-get update
sudo apt-get install postgresql-9.1-postgis

Next we should create a new template database (optional but recommended).

createdb -E UTF8 template_postgis2
createlang -d template_postgis2 plpgsql
psql -d postgres -c "UPDATE pg_database SET datistemplate='true' WHERE datname='template_postgis2'"
psql -d template_postgis2 -f /usr/share/postgresql/9.1/contrib/postgis-2.0/postgis.sql
psql -d template_postgis2 -f /usr/share/postgresql/9.1/contrib/postgis-2.0/spatial_ref_sys.sql
psql -d template_postgis2 -f /usr/share/postgresql/9.1/contrib/postgis-2.0/rtpostgis.sql
psql -d template_postgis2 -c "GRANT ALL ON geometry_columns TO PUBLIC;"
psql -d template_postgis2 -c "GRANT ALL ON geography_columns TO PUBLIC;"
psql -d template_postgis2 -c "GRANT ALL ON spatial_ref_sys TO PUBLIC;"
createdb training -T template_postgis2

Ok now we can load a raster (see sample data download below):

raster2pgsql -s 4326 srtm_4326.tif | psql training
shp2pgsql -s 4326 -d -g geom -I places.shp places| psql training

Good - now our spatial database is ready to use - and has raster support! Here is a nice example of what you can do. The query looks up the altitude from the SRTM raster for each place listed using the ST_Value function:

select ST_Value(rast, geom, true) from places, srtm_4326;

It should produce something like this:

Doing a 'point on raster' query on our raster in QGIS

Further reading: A really good place to start is the Boston GIS cheatsheets - I am really looking forward to exploring all the new features that are available in PostGIS 2.0, thanks to all those involved in building it!

Sample data for the example listed

Another bash one-liner - load all natural earth layers into postgis in one go

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 'full' one) is provided as shapefiles. Of course as a PostGIS fan I'd rather have these in my geospatial database, so I wrote a one liner in bash to load all the datasets.

First download the full dataset from here, and then extract the contents into a working directory somewhere. Then simply run this line, adjusting your destination database and schema as needed:

for FILE in `find . -name *.shp`; do \
  BASE=`basename $FILE .shp`; \
  /usr/lib/postgresql/9.1/bin/shp2pgsql -s 4326 -I $FILE world.$BASE \
  | psql gis; done

Ok that looks like four lines above, but only because I added some line breaks so things would be readable on this blog.

Note that you should replace 'world' with the name of the schema you are importing the data into, and 'gis' with the name of the database you are importing your data into. Thanks for the awesome dataset Natural Earth Team!

Listing the number of records in all postgresql tables

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.

Fun with GeoNames

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.

A week in Tanzania

I spent most of last week in Dar es Salaam, Tanzania. A lovely tropical country in the heart of Africa. I was there as part of a project I am working to create tools for Biodiversity Informatics practitioners. Of course the tools are based on Free Software:Quantum GIS and openModeller.

The attendees at the workshop were entertained by my talk about what FOSS is and why it is important, an introduction to QGIS slideshow (superbly presented by Marco Hugentobler), and ending with a tour of openModellerDesktop. We also did some live demonstrations of QGIS and openModeller, before going on to discuss details about how these tools can be used to support their Biodiversity Informatics workflows.

The meeting was funded by the Global Biodiversity Information Facility (GBIF) with Juan Bello as their representitive, and hosted by the Tanzanian Commission for Science and Technology (COSTECH).

In case you are unfamiliar with the aims of GBIF, they are facilitating the digitisation (or digitization for our american readers) of the worlds biodiversity records - herbarium records, museum collections and so on. COSTECH provides the local infrastructure and staff for the 'TanBif' node in Tanzania.

The meeting also included 'in-country' experts in the fields of GIS, Meteorology, Ecology, IT and so on. I think for all of the attendees, the concept of FOSS was a real eye-opener. African economies can't compare with those in Europe and the USA and the capital outlay for proprietary software that presents an irritation in the Western world is a major burden in the third world. So just knowing that they could dive in and use QGIS was a great revelation.

We finished our workshop a little early on the Friday so Marco and I offered to go along to the COSTECH offices and geo-enable their PostgreSQL species occurrence database and install QGIS on their desktop PC's running Windows XP. In the space of a couple of hours we were done - the major part of which was spent showing the TanBif staff members how to bring up the PostGIS layer in QGIS, perform simple queries and make maps. Having spent days in the past trying to get proprietary software like Oracle and Arc*** configured, optimised, licensed and generally usable, I was struck by just how easy and quick it is to get someone up and running with a robust enterprise ready PostGIS geospatial datastore and a user friendly Free Software desktop GIS like QGIS.

Thanks to the friendly Tanzanian folks for their hospitality - I look forward to my next visit! Here are some piccies from the trip...

image0

Juan Bello telling us about the cool things you can do with a good Biodiversity Information repository.

image1

The workshop attendees (Marco and Juan out of shot)

image2

Marco showing Godfrey how to use QGIS to bring up their PostGIS Biodiversity dataset.

image3

Godfrey proudly showing off his first map (made with QGIS)!

image4

Marco killing a mosquito - he became something of an expert!

Overpainting with Mapnik

The problem

I've been having a little poke around with Mapnik today (awesome software!). One of the things on my todo list has been to sort out rendering issues with roads we have been having. Our last iteration described roads something like this:

A style...

<Style name="Freeway30th_style">
    <Rule>
        <LineSymbolizer>
            <CssParameter name="stroke">rgb(169,170,153)</CssParameter>
            <CssParameter name="stroke-width">12.26</CssParameter>
            <CssParameter name="stroke-linejoin">bevel</CssParameter>
            <CssParameter name="stroke-linecap">round</CssParameter>
            <CssParameter name="stroke-opacity">1</CssParameter>
        </LineSymbolizer>
        <LineSymbolizer>
            <CssParameter name="stroke">rgb(255,172,88)</CssParameter>
            <CssParameter name="stroke-width">12.16</CssParameter>
            <CssParameter name="stroke-linejoin">miter</CssParameter>
            <CssParameter name="stroke-linecap">round</CssParameter>
        </LineSymbolizer>
    </Rule>
</Style>

...and this layer definition...

<Layer name="Freeway30th" srs="+init=epsg:&srid;" maxzoom="39105.90277777778">
  <StyleName>Freeway30th_style</StyleName>
  <Datasource>
      <Parameter name="dbname">&dbname;</Parameter>
      <Parameter name="estimate_extent">0</Parameter>
      <Parameter name="extent">&extent;</Parameter>
      <Parameter name="geometry_field">&geometry_field;</Parameter>
      <Parameter name="host">&host;</Parameter>
      <Parameter name="password">&password;</Parameter>
      <Parameter name="port">&port;</Parameter>
      <Parameter name="srid">&srid;</Parameter>
      <Parameter name="table">(SELECT * FROM "l_roads" WHERE "type" = \
      'Freeway' ORDER BY LENGTH(&geometry_field;) DESC) as "l_roads"</Parameter>
      <Parameter name="type">&datasourcetype;</Parameter>
      <Parameter name="user">&password;</Parameter>
  </Datasource>
</Layer>

With the idea being to render freeways with a gray outline and orange center. Unfortunately, it doesnt produce good results:

image0

The problem being those little line ends you see making gray splodges at the end of each segment.

The solution

Michael Migurski's blog discusses this issue a little in this article but doesnt directly explain how to achieve the desired effect. So here is what you do:

First the styles are split into two...

<Style name="Freeway30th_style-bottom">
   <Rule>
       <LineSymbolizer>
           <CssParameter name="stroke">rgb(169,170,153)</CssParameter>
           <CssParameter name="stroke-width">12.26</CssParameter>
           <CssParameter name="stroke-linejoin">bevel</CssParameter>
           <CssParameter name="stroke-linecap">round</CssParameter>
           <CssParameter name="stroke-opacity">1</CssParameter>
       </LineSymbolizer>
   </Rule>
 </Style>
 <Style name="Freeway30th_style-top">
   <Rule>
     <LineSymbolizer>
       <CssParameter name="stroke">rgb(255,172,88)</CssParameter>
       <CssParameter name="stroke-width">12.16</CssParameter>
       <CssParameter name="stroke-linejoin">miter</CssParameter>
       <CssParameter name="stroke-linecap">round</CssParameter>
     </LineSymbolizer>
   </Rule>
 </Style>

and then the layer is now rendered as two layers, the bottom layer first, then the top:

<Layer name="Freeway30th-bottom" srs="+init=epsg:&srid;" maxzoom="39105.90277777778">
    <StyleName>Freeway30th_style-bottom</StyleName>
    <Datasource>
        <Parameter name="dbname">&dbname;</Parameter>
        <Parameter name="estimate_extent">0</Parameter>
        <Parameter name="extent">&extent;</Parameter>
        <Parameter name="geometry_field">&geometry_field;</Parameter>
        <Parameter name="host">&host;</Parameter>
        <Parameter name="password">&password;</Parameter>
        <Parameter name="port">&port;</Parameter>
        <Parameter name="srid">&srid;</Parameter>
        <Parameter name="table">(SELECT * FROM "l_roads" WHERE "type" = \
        'Freeway' ORDER BY LENGTH(&geometry_field;) DESC) as "l_roads"</Parameter>
        <Parameter name="type">&datasourcetype;</Parameter>
        <Parameter name="user">&password;</Parameter>
    </Datasource>
</Layer>
<Layer name="Freeway30th-top" srs="+init=epsg:&srid;" maxzoom="39105.90277777778">
    <StyleName>Freeway30th_style-top</StyleName>
    <Datasource>
        <Parameter name="dbname">&dbname;</Parameter>
        <Parameter name="estimate_extent">0</Parameter>
        <Parameter name="extent">&extent;</Parameter>
        <Parameter name="geometry_field">&geometry_field;</Parameter>
        <Parameter name="host">&host;</Parameter>
        <Parameter name="password">&password;</Parameter>
        <Parameter name="port">&port;</Parameter>
        <Parameter name="srid">&srid;</Parameter>
        <Parameter name="table">(SELECT * FROM "l_roads" WHERE "type" = \
        'Freeway' ORDER BY LENGTH(&geometry_field;) DESC) as "l_roads"</Parameter>
        <Parameter name="type">&datasourcetype;</Parameter>
        <Parameter name="user">&password;</Parameter>
     </Datasource>
</Layer>

image1

A much cleaner rendering!

Note

This approach consumes more cpu time and hits your database harder than the 'messier' approach shown first.

Also you can see in the example above, I have adopted Michaels approach of rendering long lines first.

Have fun with your mapnik maps!

Introduction to PostGIS

Horst and I are spending the week up in Johannesburg at the Satellite Applications Center in Hartebeeshoek. We are doing yet another week long training course (I hope I'm not working the poor guy too hard :-P ). This time we are doing:

- Two days QGIS (with a little GRASS)
- One day PostGIS
- Two days geospatial programming with Bash, Python and QGIS

Tomorrow we start with the PostGIS component. Horst and I have been compiling some course notes for the PostGIS module which we are making available to the world as per usual. The pdf still has some rendering issues - we are aware of that. The document tries to walk the reader through the basics of using SQL and then some basic activities with PostGIS and working with geometries.

I hope some of you out there find it useful - let us know if you do! Also if you have any improvements to make, we'd love to hear from you.

Here is a quick pic or two from the course:

image0
image1
image2
image3

Back to Top

Sustaining Members