File: m.proj.py

package info (click to toggle)
grass 8.4.2-1
  • links: PTS, VCS
  • area: main
  • in suites: forky, sid
  • size: 277,040 kB
  • sloc: ansic: 460,798; python: 227,732; cpp: 42,026; sh: 11,262; makefile: 7,007; xml: 3,637; sql: 968; lex: 520; javascript: 484; yacc: 450; asm: 387; perl: 157; sed: 25; objc: 6; ruby: 4
file content (301 lines) | stat: -rwxr-xr-x 8,243 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
#!/usr/bin/env python3

"""
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-2019 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
# %rules
# % required: coordinates, input
# % exclusive: coordinates, input
# % exclusive: proj_in, -i
# % exclusive: proj_out, -o
# % exclusive: -i, -o
# %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 first: \
            https://proj.org"
            )
        )

    # 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:
        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()