#!/usr/bin/env python """ MODULE: m.proj AUTHOR(S): M. Hamish Bowman, Dept. Marine Science, Otago University, 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). #% keyword: miscellaneous #% keyword: projection #% keyword: transformation #%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: '-' for 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.utils import separator, parse_key_val, encode, decode from grass.script import core as gcore 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, ' ') line = encode(line) 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 gcore.find_program('cs2cs'): gcore.fatal(_( "cs2cs program not found, install PROJ.4 first: \ http://proj.maptools.org")) # check for overenthusiasm if proj_in and ll_in: gcore.fatal(_("Choose only one input parameter method")) if proj_out and ll_out: gcore.fatal(_("Choose only one output parameter method")) if ll_in and ll_out: gcore.fatal(_("Choose only one auto-projection parameter method")) if output and not gcore.overwrite() and os.path.exists(output): gcore.fatal(_("Output file already exists")) if not coords and not input: gcore.fatal(_("One of and must be given")) if coords and input: gcore.fatal(_( "Options and 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 = separator(ifs) ofs = separator(ofs) # set up projection params s = gcore.read_command("g.proj", flags='j') kv = parse_key_val(s) if "XY location" in kv['+proj'] and (ll_in or ll_out): gcore.fatal(_("Unable to project to or from a XY location")) in_proj = None if ll_in: in_proj = "+proj=longlat +datum=WGS84" gcore.verbose( "Assuming LL WGS84 as input, current projection as output ") if ll_out: in_proj = gcore.read_command('g.proj', flags='jf') if proj_in: if '+' in proj_in: in_proj = proj_in else: gcore.fatal(_("Invalid PROJ.4 input specification")) if not in_proj: gcore.verbose("Assuming current location as input") in_proj = gcore.read_command('g.proj', flags='jf') in_proj = in_proj.strip() gcore.verbose("Input parameters: '%s'" % in_proj) out_proj = None if ll_out: out_proj = "+proj=longlat +datum=WGS84" gcore.verbose( "Assuming current projection as input, LL WGS84 as output ") if ll_in: out_proj = gcore.read_command('g.proj', flags='jf') if proj_out: if '+' in proj_out: out_proj = proj_out else: gcore.fatal(_("Invalid PROJ.4 output specification")) if not out_proj: gcore.fatal(_("Missing output projection parameters ")) out_proj = out_proj.strip() gcore.verbose("Output parameters: '%s'" % out_proj) # set up input file if coords: x, y = coords.split(',') tmpfile = gcore.tempfile() fd = open(tmpfile, "w") fd.write("%s%s%s\n" % (x, ifs, y)) fd.close() inf = open(tmpfile) else: if input == '-': infile = None inf = sys.stdin else: infile = input if not os.path.exists(infile): gcore.fatal(_("Unable to read input data")) inf = open(infile) gcore.debug("input file=[%s]" % infile) # set up output file if not output: outfile = None outf = sys.stdout else: outfile = output outf = open(outfile, 'w') gcore.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 = gcore.Popen(cmd, stdin=gcore.PIPE, stdout=gcore.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 = decode(line).split(' ', 1) x, y = xy.split('\t') except ValueError: gcore.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 = decode(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: gcore.warning(_( "Projection transform probably failed, please investigate")) if __name__ == "__main__": options, flags = gcore.parser() main()