Well, I might be a little odd, but I like the charts produced by OpenOffice 3.2 much more than those produce by Microsoft Office 2003. But I like word more than writer, so how to combine the two?

I couldn’t find any instructions online that worked for me.¬† These did though ūüôā

Copying and pasting charts from OpenOffice Calc to Microsoft Word is a recipe for disaster.  The only effective method I have found is the follow:

  1. Copy OpenOffice Calc chart
  2. Open a new document in OpenOffice Draw.  Select Edit -> Paste Special -> GDI Metafile
  3. Right-click on the image, and select Save as Picture…
  4. Save the picture as a Windows Metafile (*.wmf), not an Enhanced Metafile (*.emf)
  5. Import into Microsoft Word via Insert Picture -> From File..
  6. Sigh with relief

The resulting image should be a vector file that word can print in a nice, antialiased way.

Advertisements

In a previous installment I demonstrated how to use ogr2ogr to reproject shapefiles. It turns out that its a bit trickier when you are reprojecting from AGD66 to GDA94 (and vice versa). This is because the ellipsoids they use to model the earth are different, and it turns out the be quite difficult to calculate the exact equivalent locations between them.

Difference between AGD66 and GDA94 at the ANU

Difference between AGD66 and GDA94 at the ANU

The most accurate way to convert between the two coordinate reference systems is to use a pre-calculated distortion grid. Fortunately, some clever surveyors have made one for us, available here (Webcite archive).

To download the distortion grid and put it in the right place for us to use, you can paste the following into your bash shell:

wget http://www.icsm.gov.au/gda/gdatm/national66.zip
unzip national66.zip
mkdir -p ~/bin
mv "A66 National (13.09.01).gsb" ~/bin/a66_national.gsb
rm national66.zip

I’ve also archived the grid files using Webcite, if the ICSM link is dead replace please see the archive.

Please note that if you change the filenames to something that contains spaces or funky punctuation ogr2ogr might not find it, and will not perform the grid transformation, leaving you with incorrect coordinates. So I would recommend just running the above code verbatim.

Then, to transform a shapefile from AGD66 to GDA94 type the following:

ogr2ogr -f "ESRI Shapefile" -s_srs "+proj=longlat +ellps=aust_SA +nadgrids=~/bin/a66_national.gsb +wktext" -t_srs EPSG:4283 outputgda94.shp inputagd66.shp

You will need to change the filenames outputgda94.shp and inputagd66.shp to whatever suits you.

Yesterday’s post expounded the virtues of doing GIS analysis within SQL transactions. One annoyance of that approach is that whenever an error is made within a transaction, the whole thing must be rolled back and redone. Any further queries won’t be processed within the transaction:

ERROR: current transaction is aborted, commands ignored until end of transaction block

Fortunately, this can be avoided, at least within recent versions of psql. All that needs to be done is to set a special variable from within psql and commands can continue to be issued as per normal:

gis=# \set ON_ERROR_ROLLBACK interactive

This variable is unset everytime you close psql. If you want this to be default behaviour within psql, close your psql session and issue the follow command from the shell:

fmark@fmark-laptop:~$ echo '\set ON_ERROR_ROLLBACK interactive' >> ~/.psqlrc

Thanks to the hivemind at stackoverflow.com for this tip.

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.

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

« Previous Page