Linfiniti Geo Blog

GIS for Open Source People

Browsing Posts in QGIS

Sometimes I make changes to my django models that causes syncdb to not run and I can’t figure out why. Here is an example of an error I might get:

python manage.py syncdb
.
. ommitted for brevity
.
    return self.get_query_set().get_or_create(**kwargs)
  File "/usr/lib/pymodules/python2.6/django/db/models/query.py",
    line 343, in get_or_create
    raise e
psycopg2.IntegrityError: duplicate key value violates unique
    constraint "auth_permission_pkey"

These messages are a bit opaque so one useful technique I found is to temporarily enable postgresql logging to a file:

sudo vim /etc/postgresql/8.3/main/postgresql.conf

Now add or uncomment the following:

log_destination = 'stderr'
logging_collector = on
log_directory = '/tmp'
log_statement = all

Next you need to restart PostgreSQL (obviously you dont want to do this on a production system without taking the needed precautions!).

sudo /etc/init.d/postgresql-8.3 restart

Next I ran my syncdb command again and at the same time watched my PostgreSQL logs being generated in /tmp (I used the unix tail command to do that):

2010-03-29 14:29:22 SAST LOG:  statement: INSERT INTO "auth_permission"
("name", "content_type_id", "codename") VALUES (E'Can add Clip', 64, E'add_clip')
2010-03-29 14:29:23 SAST ERROR:  duplicate key value violates unique
constraint "auth_permission_pkey"
2010-03-29 14:29:23 SAST STATEMENT:  INSERT INTO "auth_permission"
("name", "content_type_id", "codename") VALUES (E'Can add Clip', 64, E'add_clip')
2010-03-29 14:29:23 SAST LOG:  statement: ROLLBACK TO SAVEPOINT s1217227072_x1

Ok from that I could see that django was having trouble inserting entries into the auth_permissions table for my Clip model. The postgresql logging functions can prove indispensible when you just want to know what low level requests and responses are coming and going from your database.

Batch importing data into PostGIS is a topic I have covered previously in my article on importing CDSM (Chief Directorate: Surveys and Mapping) data into postgis. However its a popular question I get asked so today I am going to take another look at the problem. In particular the client has a directory heirachy of zipped shapefiles divided by regions in Africa. The shapefiles use a common naming convention which makes things easier. The convention is basically something like:

<number><region_name>/<region_abbreviation>_master/<region_abbreviation>_<feature_type>_<feature name>.zip

For example:

02 Mozambique Malawi/moz_master/moz_a_parks.zip

When unzipped the shapefiles are consistently named across regions e.g.:

moz_a_parks.shp
moz_a_parks.shx
moz_a_parks.dbf
moz_a_parks.sbn
moz_a_parks.sbx

From that we can define a few rules:

  • Recurse the top level directory looking for all zip files
  • Get the base name of the zip file without its extension (e.g. moz_a_parks)
  • Strip the region off the base name – this will give us the standard feature class name / table name used across all regions for that feature type.
  • Extract the zip file to /tmp
  • If the feature class / table does not exist in PostGIS, import it with the -c (create) option
  • Otherwise import it with the -a (append) option

So here is the little script I wrote that achieves this.

#!/bin/bash
DATABASE=foo
IFS="$(echo -e '\n\r')";
FILES=$(find -name "*.zip");
for FILE in $FILES
do
  BASE=$(basename $FILE .zip)
  TABLE=$( echo $BASE | sed 's/^[a-z]*_//g')
  # double check the extracted files have been cleaned up or zip will moan
  rm /tmp/${BASE}.*
  echo "$BASE to $TABLE"
  unzip -d /tmp $FILE
  MATCH=$(echo "\d" | psql ${DATABASE} | grep -o "${TABLE}")
  if [ "$MATCH" ]
  then
    # Append to the exisitng table
    echo "Appending to ${TABLE}"
    shp2pgsql -s 4326 -I -S -a -W UTF-8 "/tmp/${BASE}.shp" $TABLE | psql $DATABASE
  else
    # Create the table
    echo "Creating $TABLE"
    shp2pgsql -s 4326 -I -S -c -W UTF-8 "/tmp/${BASE}.shp" $TABLE | psql $DATABASE
  fi
  # clean up
  rm /tmp/${BASE}.*
done

If your shapefiles are already unzipped into a directory structure, you could easily modify the script above to find ‘.shp’ files rather than zip files and leave out hte unzipping , cleaning up steps.

Once you have run the scripts, open your favourite PostGIS client (I use psql) and take a look at the tables created there. If all looks good you should be able to start browsing the datasets with QGIS or publishing them with Mapserver!

Last week a QGIS hackfest was helf in Piza, Italy. Unfortunately I couldn’t attend but luckily there is a nice write up on the QGIS home page! Hopefully I will make it along to the next one!

When we released QGIS 1.4 in January, I decided to host the QGIS-1.4.0-1-No-GrassSetup.exe download file here on my server at Linfiniti so I could try to get a sense of how many people are using our software. Estimating users in a FOSS project is a hard thing to do – our software is free and freely redistributable so one copy of the software may make its way onto many desktops. Also some operating system’s packaging technology make it hard to track download numbers – we cant easily see how many downloads have come out of the ubuntu-gis PPA for example. But for the user profile that runs windows and doesn’t want GRASS functionality, I can reveal that the total downloads from my server for the months of Jan and Feb is:

22461

To which I can only say: “Wow! I didn’t know that many people use Windows!”

This has been done before by others, but its always interesting to look at image quality versus size when starting a new project. Each project has different types of data and thuse different optimal configuration. Here is a simple breakdown of size vs quality I did for a project.

First let me show you my mapserver image format definitions…:

  OUTPUTFORMAT
    NAME 'AGG_Q'
    DRIVER AGG/PNG
    IMAGEMODE RGB
    FORMATOPTION "QUANTIZE_FORCE=ON"
    FORMATOPTION "QUANTIZE_DITHER=OFF"
    FORMATOPTION "QUANTIZE_COLORS=256"
  END 

  OUTPUTFORMAT
    NAME 'AGG_JPEG'
    DRIVER AGG/JPEG
    IMAGEMODE RGB
  END 

  OUTPUTFORMAT
    NAME          "AGG_PNG"
    DRIVER        "AGG/PNG"
    IMAGEMODE     RGB
    TRANSPARENT   ON
    EXTENSION     "png"
    FORMATOPTION  "INTERLACE=OFF"
  END 

  OUTPUTFORMAT
    NAME          "AGGA_PNG"
    DRIVER        "AGG/PNG"
    IMAGEMODE     RGBA
    TRANSPARENT   ON
    EXTENSION     "png"
    FORMATOPTION  "INTERLACE=OFF"
  END 

  OUTPUTFORMAT
    NAME "GD_JPEG"
    DRIVER "GD/JPEG"
    MIMETYPE "image/jpeg"
    IMAGEMODE RGB
    EXTENSION "jpg"
  END

  OUTPUTFORMAT
    NAME "GD_PNG"
    DRIVER "GD/PNG"
    MIMETYPE "image/png"
    IMAGEMODE RGB
    EXTENSION "jpg"
  END

Now lets look at how each performed in terms of rendering quality versus file size:

File Size : 10.5kb Format: RGB JPEG using GD renderer

File Size : 19.7kb Format: RGB Png using GD renderer

File Size : 11.9kb Format: RGB Quantised Png using AGG renderer

File Size : 37.6kb Format: RGB Standard PNG using AGG renderer

File Size : 37.6kb Format: RGBA Png using AGG renderer

To be fair I haven’t added the needed ANTIALIAS clauses for the GD renderer to produce anti-aliased images for true comparison with AGG. That said even GD at it’s best isn’t a patch on AGG in my opinion. Needless to say In this project I’ve gone with the Quantised AGG outputs which are almost as small as GD JPEG images but with much nicer quality.

Note: Updated 1 March to inline images rather having to click to view the full size.

Hi my name is Robert. I am an intern at Linfiniti Consulting. I am having a great time at the company, exploring different Open Source GIS technologies. I would like to take a moment to discuss one of the projects i have been working on. I created a simple mapserver project using a mapfile that I exported from qgis. I used data that covers Mbizana municipality. MapServer is an open source development environment for building spatially-enabled internet applications. The following is a step-by-step explanation of how I did it:

1.Install a maperver export plugin on qgis (if not installed)
Click Plugins => Manage plugins
This will take you to the Qgis plugin manager window
Click on the checkbox next to Mapserver Export(Version 0.4.3)
Then click OK
2.Add some layers to your qgis
I added two layers (roads layer and mbizana_munic layer) using the data of Mbizana municipality
3.Exporting mapfile
I first saved my qgis project in /home/robert/mapserver
Then you click Plugins => Mapserver export => Mapserver export (will open the mapserver export window)
Click the “save as” and give the name to the mapfile.
I saved mine as mbizana.map in the mapserver folder (/home/robert/mapserver/).
Then click OK.
5.Edit the mapfile
In the mapfile just comment the lines that point to the symbology file and font file, to look like this:
#SYMBOLSET “./symbols/symbols.sym”
#FONTSET “./fonts/fonts.list”
6.Test your mapserver with your mapfile
I put this /localhost/cgi-bin/mapserv?map=/home/robert/mapserver/mbizana.map&mode=map

Now i can access the mapserver using Openlayers. Happy mapping !!!!!!!!!!!!

When Internet connection is a limited resource, a well-designed website doesn’t perform multiple times the same request. This little adjustment can significantly reduce the time required to load and refresh a page. First-world programmers should keep this in mind, or better come to South Africa and experience it in person…

This reminds me how life forms adapt to severe environmental conditions. But it’s a wide topic.

Let’s see how you can easily do this in Python. The snippet is a generic version of a function in views.py of the GeoDjango website we are busy developing. That function caches the result of WMS requests for layer legends in a dedicated directory, assuming that the images are not changing over time.

import urllib
import os

def retrieveDataFromUrl():
 myFileName = "file.txt"
 myLocalPath = os.path.join( os.getcwd(), myFileName )

  if not os.path.exists( myLocalPath ):
    print "Downloading data"
    myUrl = "http://whatever"
    # save it where it should have been found
    urllib.urlretrieve(myUrl, myLocalPath)
  else:
    print "Reading from local file" + myLocalPath
  # then read the file...

This code only checks if the file exists. If the file downloaded in previous run is outdated, then the newer version must be downloaded. This can be a good task for a cronjob – but it’s the topic of another post ;)

Hope this helps!

– anne

Since I got a crackberry cellphone that geotags images when I take them using its built in GPS, I have become a photo snapping lunatic. Of course one’s lunacy needs some way to manifest itself using FOSSGIS so I am going to decribe here a workflow to enable you to enjoy the thrills of geotagged imagery in QGIS even if you aren’t the proud owner of a completely proprietary crackberry phone! The process I describe below will allow you to geotag a directory full of images using a gpx track that you collect using your gps when you are out in the field. If your camera / phonecam already geotags your images, you can do step 1 and then skip straight down to step 10 below!

Step 1: Install dependencies:

sudo apt-get install gpscorrelate exiv2 python-pyexiv2

Step 2: Clean up your gpx file.

Be careful of using the QGIS editing tools to cut away any ropey looking or unneeded features – QGIS will delete the vertex attributes in the track! Rather use a text editor (like VIM!).

Step 3: Clean any existing geotags from the images:

gpscorrelate -r *.jpg;

Step 4: Find the timestamp of your first image:

exiv2 IMG_0594.jpg | grep timestamp

Image timestamp : 2009:07:31 13:02:53 IMG_0594.jpg

Step 5: Find the timestamp of the first trkpnt in the gpx file that matches the position where the photo was taken.


Step 6: Compute the time offset betweenn your camera and your gps

13 02 53 <-- photo
10 58 29 <-- gpx
 2  4 24 <-- difference

+ 2:00 hours gmt <-- timezone offset
-264 photo offset <-- difference (4min 24sec) in seconds, negative to show gps is behind the camera

Another example:

Here was a point at which we knew we had taken a picture (and we know the filename of the picture):


  685.07
  
  ACTIVE LOG
  22:10 11-Jan-10
  eTrex Venture
   0
   0

And then we got the exiv2 dump for the image:

File name       : IMG_2709.jpg
File size       : 960453 Bytes
MIME type       : image/jpeg
Image size      : 2816 x 1880
Camera make     : Canon
Camera model    : Canon EOS 1000D
Image timestamp : 2009:11:08 08:55:14
Image number    :
Exposure time   : 1/250 s
Aperture        : F8
Exposure bias   : 0 EV
Flash           : No, compulsory
Flash bias      : 0 EV
Focal length    : 70.0 mm
Subject distance: 0
ISO speed       : 100
Exposure mode   : Landscape mode
Metering mode   : Multi-segment
Macro mode      : Off
Image quality   : Normal
Exif Resolution : 2816 x 1880
White balance   : Auto
Thumbnail       : image/jpeg, 7605 Bytes
Copyright       :
Exif comment    :

From that we can work out something like this:

08 55 14 <-- photo
06 49 54 <-- gpx
 2 05 20 <-- difference
 +2hrs   <-- gmt offset
  320s   <-- second offset

Step 7: GeoTag your images using gpscorrelate:

For the first example:

gpscorrelate --timeadd +2:00 -g /tmp/knp_gpx_tim_interp.gpx --degmins -O -264 -m 1 *.jpg

Sometimes you need to make some manual adjustments to get the timestamps just right. For the images in our second example we ended up using:

gpscorrelate --timeadd +2:00 -g ../../KZN_tim.gpx --degmins -O -348 -m 30 *.jpg

Step 8: Exiv2 Image GeoTag display:

This is just to show what is going on internally - you can use the -pa option to exiv2 to see the geotags that have now been embedded in the images e.g.:

exiv2 -pa IMG_0700.jpg

produces output which includes something like this:

Exif.Image.GPSTag                            Long        1  8984
Exif.GPSInfo.GPSVersionID                    Byte        4  2.0.0.0
Exif.GPSInfo.GPSLatitudeRef                  Ascii       2  South
Exif.GPSInfo.GPSLatitude                     Rational    3  25deg 4.26000'
Exif.GPSInfo.GPSLongitudeRef                 Ascii       2  East
Exif.GPSInfo.GPSLongitude                    Rational    3  31deg 12.74000'
Exif.GPSInfo.GPSAltitudeRef                  Byte        1  Above sea level
Exif.GPSInfo.GPSAltitude                     Rational    1  675.5 m
Exif.GPSInfo.GPSTimeStamp                    SRational   3  13:25:59
Exif.GPSInfo.GPSMapDatum                     Ascii       7  WGS-84
Exif.GPSInfo.GPSDateStamp                    Ascii      11  2009:07:31

Step 9: Bonus material

gpscorrelate includes an interpolation option (enabled by default) that will derive a relative position between two vertices based on the timestamps. You can also use gpsbabel to interpolate extra vertices into your gpx file. For example I put in a vertex for every 1 second interval into the gpx file:

gpsbabel -i gpx -f knp_gpx_tim.gpx -x interpolate,time=1 -o gpx -F knp_gpx_tim_interp.gpx

Step 10: Loading in QGIS

GeoTagged image browsing with evis in QGIS

I have written a simple plugin which will create a shapefile from a directory full of geotagged images. You need to have exiv2 and python bindings for exiv2 in order for it to work, and the geotags must currently be written in degrees / decimal minutes format (which is achieved by using the --degmins option in the examples above). The plugin can be downloaded here. Point the plugin at your directory of geotagged images (it will recurse nested directories too) and afterwards a shapefile will appear in your map view. You can then query the image locations using the EVIS query tool (the EVIS plugin is included by default with QGIS).

You may also want to check out the Easy GeoTagger Project for manual options for embedding geotags into images. EasyGT also includes a vector driver which will treat a directory of geotagged images as vector layer in QGIS.

By the time a release of QGIS makes it out the door, most of us developers have long since forgotten about it and taken to the green fields of being able to add features again in QGIS trunk (i.e. ‘the fun place’). There is however a kind of sigh of relief to have made another milestone in the project and to have reduced the delta between what users want in a desktop GIS and what we provide. This (1.4) release announcement resulted in a bit of a mob rush to the QGIS website which caused quite a bit of downtime for the site over the last 24 hours. Thankfully Chris Schmidt and Frank Warmerdam (and probably others) were on hand to give the server the needed poke with a pointy stick to get it to behave better.

I decided to host the QGIS standalone installer exe on linfiniti.com for this release in order to try to get a better idea of download stats. Despite the server outages (meaning people didn’t have the link to download QGIS off my server) we have done around 900 downloads in the 24 hour period since the announcement. Considering that many people will share a download amongst friends and colleagues, especially here in South Africa where bandwidth is limited, and that there are numerous types of QGIS packages which aren’t tracked, 900 downloads probably means many times that actually installed onto people’s desktops.

QGIS 1.5 should be the final release before we start breaking API compatibility to make way for the 2.0 release. Breaking API compatibility lets use get rid of cruft from the code base and refactor the way we have designed QGIS without the overhead of having to support a large number of existing plugins and custom apps based on QGIS. I am looking forward to the rest of this year and all the new QGIS goodies the developer team will bring to the table!

I just sent out the official release announcement for QGIS 1.4 ‘Enceladus’….lets hope everyone loves it! The QGIS team has done a brilliant job for this one and I can’t wait for 1.5 because things are only going to get better! Read the QGIS Blog for details.

I put the standalone windows installer on my own server this time so I can track downloads – it will be interesting to see just how many people grab it….I’ll post some stats here later if there is anything worth reporting…