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

#!/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 <coordinates> and <input> must be given"))
    if coords and input:
        gcore.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 = 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()