Archive for April, 2010

Yep, its one of those small annoyances – like keeping up-to-date backups – that you should always do: always use transactions when doing analysis in PostGIS.

In case you don’t know, transactions let you wrap a bunch of database modifications such as UPDATE, INSERT, CREATE, etc. in an undo-able block. This means you can recover from mistakes, like you would use undo when using Excel. When you are happy with your changes, you can then commit them to the database to store them permanently.

You start a transaction by issuing a BEGIN comand, and end it by either undoing the transaction with ROLLBACK or saving it with COMMIT. Here’s an example of how to do it, using pqsl:

Let’s say you realise you have a table two columns with latitudes that are incorrectly entered as positive numbers, although they refer to locations south of the equator:

gis=# SELECT lat_to, lat_from FROM closures ORDER BY closure_id DESC LIMIT 10;

 lat_to  | lat_from 
---------+----------
       0 | -14.9819
 13.8779 | -13.9802
 13.6661 | -14.0998
 14.6906 |        0
 13.6661 | -14.5559
       0 |        0
       0 | -13.9802
       0 |        0
 25.4141 |        0
 22.2898 | -22.9968
(10 rows)

So we, issue some updates, wrapped in a transaction:

gis=# BEGIN;
BEGIN
gis=# UPDATE closures SET lat_to = (-(abs(lat_to))) WHERE lat_to != 0;
UPDATE 1866
gis=# UPDATE closures SET lat_to = (-(abs(lat_from))) WHERE lat_from != 0;
UPDATE 1912

Now, we check our result and realise we have made a mistake, by updating the wrong field in our second query:

gis=# SELECT lat_to, lat_from FROM closures ORDER BY closure_id DESC LIMIT 10;

  lat_to  | lat_from 
----------+----------
 -14.9819 | -14.9819
 -13.9802 | -13.9802
 -14.0998 | -14.0998
 -14.6906 |        0
 -14.5559 | -14.5559
        0 |        0
 -13.9802 | -13.9802
        0 |        0
 -25.4141 |        0
 -22.9968 | -22.9968
(10 rows)

If we hadn’t used a transaction, we would now be stuffed. Fortunately, we can revert the changes:

gis=# ROLLBACK;
ROLLBACK
gis=# SELECT lat_to, lat_from FROM closures ORDER BY closure_id DESC LIMIT 10;

 lat_to  | lat_from 
---------+----------
       0 | -14.9819
 13.8779 | -13.9802
 13.6661 | -14.0998
 14.6906 |        0
 13.6661 | -14.5559
       0 |        0
       0 | -13.9802
       0 |        0
 25.4141 |        0
 22.2898 | -22.9968
(10 rows)

Now we can try again, fixing the broken query:

gis=# BEGIN;
BEGIN
gis=# UPDATE closures SET lat_to = (-(abs(lat_to))) WHERE lat_to != 0;
UPDATE 1866
gis=# UPDATE closures SET lat_from = (-(abs(lat_from))) WHERE lat_from != 0;
UPDATE 1912
gis=# SELECT lat_to, lat_from FROM closures ORDER BY closure_id DESC LIMIT 10;

  lat_to  | lat_from 
----------+----------
        0 | -14.9819
 -13.8779 | -13.9802
 -13.6661 | -14.0998
 -14.6906 |        0
 -13.6661 | -14.5559
        0 |        0
        0 | -13.9802
        0 |        0
 -25.4141 |        0
 -22.2898 | -22.9968
(10 rows)

Yay, the changes we wanted have been made, so we can now commit the transaction to save the changes permanently:

gis=# COMMIT;
COMMIT
gis=# SELECT lat_to, lat_from FROM closures ORDER BY closure_id DESC LIMIT 10;

  lat_to  | lat_from 
----------+----------
        0 | -14.9819
 -13.8779 | -13.9802
 -13.6661 | -14.0998
 -14.6906 |        0
 -13.6661 | -14.5559
        0 |        0
        0 | -13.9802
        0 |        0
 -25.4141 |        0
 -22.2898 | -22.9968
(10 rows)

UPDATE:In order to avoid rolling back transactions every time you make a typo inside them, set ON_ERROR_ROLLBACK to interactive.

Advertisements

Firstly, you might want to view the projection information in the shapefile, ROAD_CENTRELINES.shp in this case.  The perl one-liner stops ogrinfo from printing copious data about each vector feature in the shapefile:

fmark@fmark-laptop:~$ ogrinfo -ro -al ROAD_CENTERLINES.shp | perl -e 'while (<STDIN>) { if (/OGRFeature/) { exit } print; }'
INFO: Open of `ROAD_CENTERLINES.shp'
using driver `ESRI Shapefile' successful.

Layer name: ROAD_CENTERLINES
Geometry: Line String
Feature Count: 4596
Extent: (128.999678, -25.998883) - (138.000000, -11.067068)
Layer SRS WKT:
GEOGCS["Longitude / Latitude (GDA 94)",
DATUM["Unknown",
SPHEROID["GRS_80",6378137,298.2572221010595]],
PRIMEM["Greenwich",0],
UNIT["Degree",0.017453292519943295]]
ROAD_NAME: String (100.0)
ROAD_TYPE: String (15.0)
ROAD_SUFFI: String (15.0)
ROAD_LABEL: String (130.0)
Length: Real (19.11)

We can see from this that the geographic coordinate system of the shapefile is GDA94.

In order to specify the projection we want to use, we need to know its EPSG code.  What’s an EPSG code?  Its a number that uniquely defines some well-known projections (according to Wikipedia, EPSG stands for the European Petroleum Survey Group, so we have big oil to thank for our daily GIS tools, as well as big military). Some useful Australian EPSG codes are:

4326 WGS84
4283 GDA94
4202 AGD66
3577 GDA94 / Australian Albers
3112 GDA94 / Geoscience Australia Lambert
28348 GDA94 / MGA zone 48
28349 GDA94 / MGA zone 49
28350 GDA94 / MGA zone 50
28351 GDA94 / MGA zone 51
28352 GDA94 / MGA zone 52
28353 GDA94 / MGA zone 53
28354 GDA94 / MGA zone 54
28355 GDA94 / MGA zone 55
28356 GDA94 / MGA zone 56
28357 GDA94 / MGA zone 57
28358 GDA94 / MGA zone 58

More can be found at http://spatialreference.org/.

In this case, I want MGA zone 53, or EPSG:28353.  So to reproject:

ogr2ogr -t_srs 'EPSG:28353' ROAD_CENTERLINES_MGAZ53.shp ROAD_CENTERLINES.shp

Note that the output filename comes before the input filename.  Also, for some reason on ubuntu karmic there was a missing dependency, so I needed to:

sudo apt-get install libproj-dev