Coordinate system conversion using cs2cs from PROJ

I hate geodesy. It's confusing, multifarious, and ever-changing (and I call myself a scientist!). But we all need to know where we, and more importantly our data, are. I have need of a single tool which painlessly does all my conversions, preferably without me having to think too much about it.

You could use an online tool, but when you're scripting (and I'm ALWAYS scripting) you want something a little less manual. Enter cs2cs from the PROJ.4 initiative. It will convert pretty much anything to anything else. Perfect.

In Ubuntu:
sudo apt-get install proj-bin

In Fedora I searched yum for proj-4

If the man page is a little too much too soon to take in, this wonderful page helps you decide which flags and parameters to use, with a rudimentary understanding of what you want

For example, if I want to convert coordinates in Arizona Central State Plane to WGS84 Latitude/Longitude, I use this:

cs2cs -f "%.6f" +proj=tmerc +lat_0=31 +lon_0=-111.9166666666667 +k=0.9999 +x_0=213360 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs infile > outfile

view raw cs2cs 1.sh hosted with ❤ by GitHub

(the State Plane Coordinate System is new to me, but my view of the whole field of geodesy is 'let's make things even more complicated, sorry, accurate!')

Let's decompose that a bit:

-f "%.6f"

view raw proj1.sh hosted with ❤ by GitHub

tells it I want decimal degrees with 6 values after the decimal point (you could leave this out if you want deg, min, sec)

+proj=tmerc

view raw proj2.sh hosted with ❤ by GitHub

says use the Transverse Mercator projection

+lat_0=31 +lon_0=-111.9166666666667 +k=0.9999 +x_0=213360 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs

view raw proj3.sh hosted with ❤ by GitHub

are all parameters related to the input coordinates in Arizona Central State Plane (the 'to_meter' is a conversion from feet, the default unit, to metres)

+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs

view raw proj4.sh hosted with ❤ by GitHub

are the parameters related to the output coordinates in Latitude/Longitude

infile is a two column list of points e.g.
2.19644131000e+005 6.11961823000e+005
2.17676234764e+005 6.11243565478e+005
2.19763457634e+005 6.11234534908e+005
2.19786524782e+005 6.11923555789e+005
2.19762476867e+005 6.11378246389e+005

outfile will be created with the results:
-111.846501    36.517758 0.000000
-111.868478    36.511296 0.000000
-111.845175    36.511202 0.000000
-111.844912    36.517412 0.000000
-111.845185    36.512498 0.000000

By default it will always give you the third dimension (height) even if you didn't ask for it, like above, so I tend to trim it by asking awk to give me just the first 2 columns separated by a tab:

awk '{print $(NF-2),"\t",$(NF-1)}' outfile > outfile2

view raw awk trim.awk hosted with ❤ by GitHub

-111.846501      36.517758
-111.868478      36.511296
-111.845175      36.511202
-111.844912      36.517412
-111.845185      36.512498