Geografic-Information-System/scripts/m.proj/m.proj.py

#!/usr/bin/env python

############################################################################
#
# MODULE:	m.proj
# AUTHOR:	M. Hamish Bowman, Dept. Marine Science, Otago Univeristy,
#		  New Zealand
#               Converted to Python by Glynn Clements
# PURPOSE:      cs2cs reprojection frontend for a list of coordinates.
#		Replacement for m.proj2 from GRASS 5
# COPYRIGHT:	(c) 2006-2014 Hamish Bowman, and the GRASS Development Team
#		This program is free software under the GNU General Public
#		License (>=v2). Read the file COPYING that comes with GRASS
#		for details.
#
#############################################################################

# notes:
#  - cs2cs expects "x y" data so be sure to send it "lon lat" not "lat lon"
#  - if you send cs2cs a third data column, beware it might be treated as "z"
# todo:
#  - `cut` away x,y columns into a temp file, feed to cs2cs, then `paste`
#    back to input file. see method in v.in.garmin.sh. that way additional
#    numeric and string columns would survive the trip, and 3rd column would
#    not be modified as z.

#%module
#% description: Converts coordinates from one projection to another (cs2cs frontend).
#% keywords: miscellaneous
#% keywords: projection
#%end
#%option G_OPT_M_COORDS
#% description: Input coordinates to reproject
#% guisection: Input coordinates
#%end
#%option G_OPT_F_INPUT
#% label: Name of input coordinate file
#% description: '-' to read from standard input
#% required: no
#% guisection: Input coordinates
#%end
#%option G_OPT_F_OUTPUT
#% description: Name for output coordinate file (omit to send to stdout)
#% required : no
#% guisection: Output
#%end
#%option G_OPT_F_SEP
#% label: Field separator (format: input[,output])
#% required : no
#% guisection: Input coordinates
#%end
#%option
#% key: proj_in
#% type: string
#% description: Input projection parameters (PROJ.4 style)
#% required : no
#% guisection: Projections
#%end
#%option
#% key: proj_out
#% type: string
#% description: Output projection parameters (PROJ.4 style)
#% required : no
#% guisection: Projections
#%end
#%flag
#% key: i
#% description: Use LL WGS84 as input and current location as output projection
#% guisection: Projections
#%end
#%flag
#% key: o
#% description: Use current location as input and LL WGS84 as output projection
#% guisection: Projections
#%end
#%flag
#% key: d
#% description: Output long/lat in decimal degrees, or other projections with many decimal places
#% guisection: Output
#%end
#%flag
#% key: e
#% description: Include input coordinates in output file
#% guisection: Output
#%end
#%flag
#% key: c
#% description: Include column names in output file
#% guisection: Output
#%end


import sys
import os
import threading
from grass.script import core as grass

class TrThread(threading.Thread):
    def __init__(self, ifs, inf, outf):
        threading.Thread.__init__(self)
        self.ifs = ifs
        self.inf = inf
        self.outf = outf

    def run(self):
        while True:
            line = self.inf.readline()
            if not line:
                break
            line = line.replace(self.ifs, ' ')
            self.outf.write(line)
            self.outf.flush()

        self.outf.close()

def main():
    coords = options['coordinates']
    input = options['input']
    output = options['output']
    fs = options['separator']
    proj_in = options['proj_in']
    proj_out = options['proj_out']
    ll_in = flags['i']
    ll_out = flags['o']
    decimal = flags['d']
    copy_input = flags['e']
    include_header = flags['c']

    #### check for cs2cs
    if not grass.find_program('cs2cs'):
	grass.fatal(_("cs2cs program not found, install PROJ.4 first: http://proj.maptools.org"))

    #### check for overenthusiasm
    if proj_in and ll_in:
	grass.fatal(_("Choose only one input parameter method"))

    if proj_out and ll_out:
	grass.fatal(_("Choose only one output parameter method")) 

    if ll_in and ll_out:
	grass.fatal(_("Choise only one auto-projection parameter method"))

    if output and not grass.overwrite() and os.path.exists(output):
	grass.fatal(_("Output file already exists")) 

    if not coords and not input:
        grass.fatal(_("One of <coordinates> and <input> must be given"))
    if coords and input:
        grass.fatal(_("Options <coordinates> and <input> are mutually exclusive"))

    #### parse field separator
    # FIXME: input_x,y needs to split on multiple whitespace between them
    if fs == ',':
        ifs = ofs = ','
    else:
	try:
	    ifs, ofs = fs.split(',')
	except ValueError:
	    ifs = ofs = fs

    ifs = grass.separator(ifs)
    ofs = grass.separator(ofs)

    #### set up projection params
    s = grass.read_command("g.proj", flags='j')
    kv = grass.parse_key_val(s)
    if "XY location" in kv['+proj'] and (ll_in or ll_out):
	grass.fatal(_("Unable to project to or from a XY location")) 

    in_proj = None

    if ll_in:
	in_proj = "+proj=longlat +datum=WGS84"
	grass.verbose("Assuming LL WGS84 as input, current projection as output ")

    if ll_out:
	in_proj = grass.read_command('g.proj', flags = 'jf')

    if proj_in:
	in_proj = proj_in

    if not in_proj:
	grass.verbose("Assuming current location as input")
        in_proj = grass.read_command('g.proj', flags = 'jf')
    
    in_proj = in_proj.strip()
    grass.verbose("Input parameters: '%s'" % in_proj)

    out_proj = None

    if ll_out:
	out_proj = "+proj=longlat +datum=WGS84"
	grass.verbose("Assuming current projection as input, LL WGS84 as output ")

    if ll_in:
	out_proj = grass.read_command('g.proj', flags = 'jf')

    if proj_out:
	out_proj = proj_out

    if not out_proj:
	grass.fatal(_("Missing output projection parameters "))
    out_proj = out_proj.strip()
    grass.verbose("Output parameters: '%s'" % out_proj)

    #### set up input file
    if coords:
        x, y = coords.split(',')
        tmpfile = grass.tempfile()
        fd = open(tmpfile, "w")
        fd.write("%s%s%s\n" % (x, ifs, y))
        fd.close()
        inf = file(tmpfile)
    else:
        if input == '-':
            infile = None
            inf = sys.stdin
        else:
            infile = input
            if not os.path.exists(infile):
                grass.fatal(_("Unable to read input data"))
            inf = file(infile)
            grass.debug("input file=[%s]" % infile)
    
    #### set up output file
    if not output:
	outfile = None
	outf = sys.stdout
    else:
	outfile = output
	outf = open(outfile, 'w')
	grass.debug("output file=[%s]" % outfile) 

    #### set up output style
    if not decimal:
	outfmt = ["-w5"]
    else:
	outfmt = ["-f", "%.8f"]
    if not copy_input:
	copyinp = []
    else:
	copyinp = ["-E"]

    #### do the conversion
    # Convert cs2cs DMS format to GRASS DMS format:
    #   cs2cs | sed -e 's/d/:/g' -e "s/'/:/g"  -e 's/"//g'

    cmd = ['cs2cs'] + copyinp + outfmt + in_proj.split() + ['+to'] + out_proj.split()
    p = grass.Popen(cmd, stdin = grass.PIPE, stdout = grass.PIPE)

    tr = TrThread(ifs, inf, p.stdin)
    tr.start()

    if not copy_input:
	if include_header:
	    outf.write("x%sy%sz\n" % (ofs, ofs))
	for line in p.stdout:
            try:
                xy, z = line.split(' ', 1)
                x, y = xy.split('\t')
            except ValueError:
                grass.fatal(line)
            
	    outf.write('%s%s%s%s%s\n' % \
                       (x.strip(), ofs, y.strip(), ofs, z.strip()))
    else:
	if include_header:
	    outf.write("input_x%sinput_y%sx%sy%sz\n" % (ofs, ofs, ofs, ofs))
	for line in p.stdout:
            inXYZ, x, rest = line.split('\t')
            inX, inY = inXYZ.split(' ')[:2]
	    y, z = rest.split(' ', 1)
	    outf.write('%s%s%s%s%s%s%s%s%s\n' % \
                       (inX.strip(), ofs, inY.strip(), ofs, x.strip(), \
		        ofs, y.strip(), ofs, z.strip()))

    p.wait()

    if p.returncode != 0:
	grass.warning(_("Projection transform probably failed, please investigate"))

if __name__ == "__main__":
    options, flags = grass.parser()
    main()