Coordinate Transforms

From Marine Science
Jump to: navigation, search

Main Page > How To Do Stuff > Coordinate Transforms

Contents

Coordinate transforms

Geographic coordinate systems

These are degree based spherical (aka polar) coordinates. The distance and bearing between two points need special care as Great-Circle calculations must be used. If you want a highly inaccurate cheat, you can remember that there are 1852m in a nautical mile, and 60 nautical miles in a degree of latitude (111km). The width of a degree longitude varies as the cosine of the latitude. So at the equator it is 111km wide and by the poles it is nil.

Latitude / Longitude WGS84

  • EPSG code 4326

(for info about EPSG codes, see below)

GEOGCS["wgs84",
   DATUM["WGS_1984",
       SPHEROID["WGS_1984",6378137,298.257223563],
       TOWGS84[0.000,0.000,0.000]],
   PRIMEM["Greenwich",0],
   UNIT["degree",0.0174532925199433]]

The WGS84 chart datum

Navigation charts and GPSs will typically churn out lat/lon coordinates. Now-a-days (since 1984) this will be used together with the WGS84 datum. This datum defines the 0,0,0 x,y,z orgin of the Earth at its centre of mass and the difference between a perfect sphere and the lumpy ellipsoid we actually live on. In general the centre-spot and gravitational lumpiness wasn't very well known before we sent up satellites, so older dautms will vary widely while newer versions will only vary be a few cm. (besides the 1984 one some common modern datums are 1996 and 2000 versions)

Older NZ charts will generally use the NZ 1949 geoid datum, which is about 100m offset from WGS84.

To avoid confusion make sure your handheld GPS is set to use WGS84. As tempting as it may be, don't ever mess with a ship's navigation GPS without the skipper looking on, they will have set it up to work a certain way and will be expecting that to be the case. Woe be the student who invokes the captain's ire.

Lat/Long or Long/Lat?

The common use is to put latitude first, but keep in mind that when plotting or converting coordinates this puts the northing first, i.e. "y x".

Degrees, minutes, and seconds

While it may seem obvious, the conversion between degrees, minutes, and seconds always seems to lead to mass confusion, often because what was used wasn't written down. (e.g. instead of punctuation just a string of numbers was written down) So always write the little =o= degree sign, the ' minute sign, and " second sign, and a clear decimal point or you'll regret it. It helps to also create a waypoint in the GPS whenever you have to write one down. You can download it later which is less work than transcribing and less error prone after a long day of staring at a spreadsheet.

The common ways of presenting lat/lon data are:

  • Decimal degrees. Like 166.6546323,-45.2587535. In this form South is replaced by a negative sign (East is positive) and to preserve precision if possible try to keep at least 6-7 numbers after the decimal place. This is the most computer friendly but least human friendly version.
  • Degrees, Minutes. Like 166d 39.2779' E, 45d 15.5252' S. This is the typical form which a handheld GPS will use. Divide minutes by 60 to get decimal degrees. Be careful in a spreadsheet not to add + minutes part to a negative latitude.
  • Degrees, Minutes, Seconds. Like 166d 39' 16.676" E, 45d 15' 31.513" S. This is the typical form which a ship's GPS will use. As with the Deg Min form, divide seconds by 60 to get decimal minutes.

Cartesian coordinate systems

These are typically meters based so distance and area calculations are much much simpler to work with. Unlike lat/lon usually they are limited in their range, the further you go away from the centre the more distortion is introduced. Coordinates are usually referred to as "eastings" and "northings" and their values will typically be in the millions.

New Zealand Map Grid (NZMG)

  • Used by tramping topo maps and older GIS data.
  • Based on the line-of-sight between mountain top trig station network, so accuracy isn't really meaningful under a meter or so. Depending on the transform method you use the error could be as much as 5m out. Generally your handheld GPS can't do much better, so unless you care about millimeter precision don't worry about it.
  • This typically will use the NZ 1949 geoid datum. There are 3 and 7 parameter transforms available for the conversion to the WGS84 datum, but if you can help it you will want to use the LINZ NTv2 distortion grid as it will give a much better result.
  • EPSG code 27200

PROJCS["New Zealand Map Grid",
   GEOGCS["international",
       DATUM["New_Zealand_Geodetic_Datum_1949",
           SPHEROID["International_1924",6378388,297]],
       PRIMEM["Greenwich",0],
       UNIT["degree",0.0174532925199433]],
   PROJECTION["New_Zealand_Map_Grid"],
   PARAMETER["latitude_of_origin",-41],
   PARAMETER["central_meridian",173],
   PARAMETER["false_easting",2510000],
   PARAMETER["false_northing",6023150],
   UNIT["meter",1]]

New Zealand Transverse Mercator (NZTM)

  • Used by newer GIS data. Centered on Nelson.
  • This typically will use the NZ 2000 geoid datum. (NZGD2k is almost identical to WGS84)
  • EPSG code 2193

PROJCS["Transverse Mercator",
   GEOGCS["grs80",
       DATUM["New_Zealand_Geodetic_Datum_2000",
           SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101],
           TOWGS84[0,0,0,0,0,0,0]],
       PRIMEM["Greenwich",0],
       UNIT["degree",0.0174532925199433]],
   PROJECTION["Transverse_Mercator"],
   PARAMETER["latitude_of_origin",0],
   PARAMETER["central_meridian",173],
   PARAMETER["scale_factor",0.9996],
   PARAMETER["false_easting",1600000],
   PARAMETER["false_northing",10000000],
   UNIT["metre",1]]

North Taieri Circuit 2000

  • Used by newer GIS data produced for the Dunedin area.
  • Like NZTM, N. Taieri Circuit is a transverse mercator projection, but it is centered in the Dunedin area to keep distortion to a minimum. Use this if working locally, millimeter accuracy is important, and you have access to fancy RTK GPS equipment from the surveying department.
  • EPSG code 2131
  • Keep an eye out for the 1949 version of this. Luckily the 2000 version was set up to be a few million away from the 1949 version so it is easy to spot if you have been handed the wrong flavour, as there is no overlap.

PROJCS["Transverse Mercator",
   GEOGCS["grs80",
       DATUM["New_Zealand_Geodetic_Datum_2000",
           SPHEROID["Geodetic_Reference_System_1980",6378137,298.257222101],
           TOWGS84[0,0,0,0,0,0,0]],
       PRIMEM["Greenwich",0],
       UNIT["degree",0.0174532925199433]],
   PROJECTION["Transverse_Mercator"],
   PARAMETER["latitude_of_origin",-45.86138888888889],
   PARAMETER["central_meridian",170.2825],
   PARAMETER["scale_factor",0.99996],
   PARAMETER["false_easting",400000],
   PARAMETER["false_northing",800000],
   UNIT["metre",1]]

Further reading

You can find an extensive collection of documents and help on the internet.

Transform software

Each map projection will have its own set of transform parameters. Often your GPS data will come with a .prj file with "Well Known Text" (WKT) describing the projection. Just like in biology there are often many names for the same thing, so the European Petroleum Survey Group (EPSG) decided to give each unique projection its own ID code to help ease the pain of translation. If your software supports it EPSG codes are the easiest way of identifying which coordinate system you are actually talking about.

Online

  • LINZ's website provides an online converter tool
    for up to 5 points.
    for Lat/Lon WGS84 choose "World Geodetic System 1984"

Multiplatform

(Mac, PC, Linux, etc.)

PROJ.4

  • PROJ.4 is the industry standard. It can be a bit tricky to get going, but it does the business. With the right recipe it will work well, hence this help page.
  • FWTools is a handy package which includes PROJ.4 but also a nice GIS data viewer called OpenEV and a massive GIS data format converter called GDAL, plus other goodies.


Using PROJ.4

The cs2cs program will convert a data file between two coordinate systems.

  • Proj will accept decimal degree or DDDdMM'SS.SSS" format.
  • Remember to present your data as "x y", ie "long lat"
From the command prompt

example:

  • Lon/Lat WGS84 to North Taieri Circuit 2000
   cs2cs -f "%.3f" +init=epsg:4326 +to +init=epsg:2131 < coord_LL.txt > coord_NTai2k.txt
  • Lat/Lon WGS84 to North Taieri Circuit 2000 (y,x)
   cs2cs -r -f "%.3f" +init=epsg:4326 +to +init=epsg:2131 < coord_LL.txt > coord_NTai2k.txt
  • New Zealand Map Grid to NZ Transverse Mercator 2000 (using default datum transform; 3-5m accuracy)
   cs2cs -f "%.3f" +init=epsg:27200 +to +init=epsg:2193 < coord_NZMG.txt > coord_NZTM.txt
  • New Zealand Map Grid to NZ Transverse Mercator 2000 (using best datum transform via the distortion grid file, download it from LINZ website, 10-30cm accuracy)
   cs2cs -f "%.3f" +proj=nzmg +lat_0=-41.0000000000 +lon_0=173.0000000000 \
      +x_0=2510000.0000000000 +y_0=6023150.0000000000 +a=6378388 +rf=297 \
      +no_defs +nadgrids=/path/to/nzgd2kgrid0005.gsb +to_meter=1.0 \
      +to  +init=epsg:2193 < coord_NZMG.txt > coord_NZTM.txt


  • NZMG to NZTM, single point
   echo "2312901  5486097" | cs2cs -f "%.3f" +init=epsg:27200 +to +init=epsg:2193
Using FWTools
  • Install FWTools(version).exe from the Bambi repository
  • Click on the FWTools Shell icon, run cs2cs from the command line as above.
From Matlab

If you prefer to use something less costly than Matlab, Octave is Free and runs in a manner "not dissimilar".

Windows users:

  • Add setfw.bat in the FWTools dir to autoexec.bat or whatever XP uses for startup commands
  • Alternatively use Matlab's setenv() command to recreate FWTools\setfw.bat and FWTools\bin\setfwenv.bat

Everyone:

  • Put coordinates in a 2 column array, with x first. (ie longitude, latitude) Make sure latitude is negative for the southern hemisphere!
  • Save coordinates to transform with save coords_LL.txt coord_array -ASCII
  • Then use !cs2cs= or =[exitcode, result] = system('cs2cs +parameters < inputfile.txt'). An exitcode of 0 is success, 1 means an error.
  • inspect the new result output, you may have to parse it with sscanf(). If saving the results to a file with cs2cs, you can use Matlab's load output.txt command to split the results into columns automatically.
From Python

Python is an easy to use and learn, yet powerful, scripting language which is quickly becoming the de facto standard for scientific programming.

  • TODO

DOS

  • LINZ puts out a DOS program called Concord. It works.
  • I think there's one call GeoCalc for MS Windows. From memory it's a bit of a pain in the neck too and the conversions aren't as exact.

Mapping software

Quantum GIS

Control of map projections in QGIS is somewhat free-form. It will let you load up the data first and worry about projections later.

  • Pressing the "P" key or Settings → Project Properties → Projection from the menu will bring up the Projection menu

TODO: how to reproject

  • See chapter 7 of the QGIS user manual

GRASS GIS

Control of map projections in GRASS can be a bit more strict than users of other GIS mapping software may be used to. GRASS will try to stop you from doing the wrong thing projection-wise. So you may see more error messages up front but your data is less likely to get silently broken somewhere along the way. This is simply a different compromise based on a different set of values.

  • Create your project location using the data source's projection settings. New locations can also be created by EPSG code or by manual selection from a list.
  • Load additional data into a location of its native projection with the v.in.ogr and r.in.gdal modules. These will load most known vector and raster data formats. These modules can also create a new project location for you if the data source's projection does not match the current location's settings.
  • .csv and ASCII text data can be imported with the v.in.ascii module.
  • GRASS's degree minute second format can use decimal degrees or DDD:MM.MMMMh or DDD:MM:SS.SSSh, where h represents the N,S,E,W hemisphere letter.
  • View your current map projection settings with ' g.proj -p '
  • Pull data in to your target location with v.proj (for vector data) and r.proj (for raster data) modules.
  • Do conversions for a series of data points with the m.proj module. The -i and -o flags are simple shortcuts for converting from and to WGS84 lat/lon from the current map projection. This module is really just a wrapper for PROJ.4's cs2cs which fills in all the complicated projection terms for you, making it much easier to use.

ArcGIS

Control of map projections in Arc is somewhat free-form. It will let you load up the data first and worry about projections later.

  • Arc can reproject on the fly, but it seems that it will not apply the datum shift even if it knows they don't match. What this means is that going to and from NZ Map Grid will introduce a >100m error! (usually points are translated north-sound) The solution to this is to use concord or cs2cs to do the reprojection + datum transform, then load those into Arc without the reprojection on the fly.
  • If you have loaded points from a .csv file and you know the projection they are in, but Arc doesn't, you can set it as follows: Right click on the layer, select Layer Properties -> Source tab -> Set Data Source... button -> Projected -> National Grids -> NZ -> then NZ Map Grid or NZ Transverse Mercator 2000. For lat/lon data choose geographic coordinate systems choose Geographic Coord system instead of projected. Note a projected map system like NZMG will still use a Geographic man datum, so you'll still see that word floating around in the projection description.
  • Good luck.


-- HamishBowman - 10 Mar 2008