Mercurial Mercurial > wsgi > hgweb.wsgi > reduce / comparison
summary | shortlog | changelog | graph | tags | bookmarks | branches | files | changeset | file | latest | revisions | annotate | diff | comparison | raw | help
Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
gpx_plot.py
changeset 10 89108adbc468
parent 4 1b96bb9de1f3
child 12 1d6e37b3ebf2
equal deleted inserted replaced
9:b379310ee21c 10:89108adbc468 
     1 #!/usr/bin/env python
 
     1 #!/usr/bin/env python
 
     2 # -*- coding: utf8 -*-
 
     2 # -*- coding: utf8 -*-
 
     3 
     3 
     4 '''
 
     4 '''
 
       
     5 gpx_plot.py calls gnuplot to draw the gpx-tracks given as arguments.
 
       
     6 Copyright 2017 willem179
 
       
     7 Extracted from the code of "gpx_reduce.py"
 
       
     8 Copyright (C) 2011,2012,2013,2015,2016,2017 travelling_salesman on OpenStreetMap
 
       
     9 
     5 This program is free software: you can redistribute it and/or modify
 
    10 This program is free software: you can redistribute it and/or modify
 
     6 it under the terms of the GNU General Public License as published by
 
    11 it under the terms of the GNU General Public License as published by
 
     7 the Free Software Foundation, either version 3 of the License, or
 
    12 the Free Software Foundation, either version 3 of the License, or
 
     8 (at your option) any later version.
 
    13 (at your option) any later version.
 
     9 
    14 
    17 '''
 
    22 '''
 
    18 
    23 
    19 import datetime
 
    24 import datetime
 
    20 import sys
 
    25 import sys
 
    21 import time
 
    26 import time
 
    22 import numpy as np
 
       
    23 import numpy.linalg as la
 
       
    24 from math import *
 
    27 from math import *
 
    25 import xml.etree.ElementTree as etree
 
    28 import xml.etree.ElementTree as etree
 
    26 from optparse import OptionParser
 
    29 from optparse import OptionParser
 
    27 
    30 
       
    31 # the path to the gnuplot binary
 
       
    32 gnuPlotCmd = 'gnuplot'
 
    28 
    33 
    29 parser = OptionParser('usage: %prog [options] input-file.gpx')
 
    34 parser = OptionParser('usage: %prog [options] input-file.gpx')
 
    30 (options, args) = parser.parse_args()
 
    35 (options, args) = parser.parse_args()
 
    31 
    36 
    32 
    37 
    40 a = 6378137.0
 
    45 a = 6378137.0
 
    41 b = 6356752.314245179
 
    46 b = 6356752.314245179
 
    42 
    47 
    43 timeformat = '%Y-%m-%dT%H:%M:%SZ'
 
    48 timeformat = '%Y-%m-%dT%H:%M:%SZ'
 
    44 
    49 
       
    50 # the linear algebra with lists
 
       
    51 norm = lambda p: sqrt (sum (a * a for a in p))
 
       
    52 
    45 def rotate(x, y, phi):
 
    53 def rotate(x, y, phi):
 
    46     return x*cos(phi) - y*sin(phi), x*sin(phi) + y*cos(phi)
 
    54     return x*cos(phi) - y*sin(phi), x*sin(phi) + y*cos(phi)
 
    47 
    55 
    48 
    56 
    49 def project_to_meters(lat, lon, latm, lonm):
 
    57 def project_to_meters(lat, lon, latm, lonm):
 
    50     # azimuthal map projection centered at average track coordinate
 
    58     # azimuthal map projection centered at average track coordinate
 
    51     lon -= lonm
 
    59     lon -= lonm
 
    52     xyz = latlonele_to_xyz(lat, lon, 0.0)
 
    60     xyz = latlonele_to_xyz(lat, lon, 0.0)
 
    53     zy = rotate(xyz[2], xyz[1], radians(90 - latm))
 
    61     zy = rotate(xyz[2], xyz[1], radians(90 - latm))
 
    54     lat2 = atan2(zy[0], la.norm([zy[1], xyz[0]]))
 
    62     lat2 = atan2(zy[0], norm([zy[1], xyz[0]]))
 
    55     lon2 = atan2(xyz[0], -zy[1])
 
    63     lon2 = atan2(xyz[0], -zy[1])
 
    56     x_meters = rE * sin(lon2) * (pi / 2.0 - lat2)
 
    64     x_meters = rE * sin(lon2) * (pi / 2.0 - lat2)
 
    57     y_meters = -rE * cos(lon2) * (pi / 2.0 - lat2)
 
    65     y_meters = -rE * cos(lon2) * (pi / 2.0 - lat2)
 
    58     return x_meters, y_meters
 
    66     return x_meters, y_meters
 
    59 
    67 
    60 
    68 
    61 def latlonele_to_xyz(lat, lon, ele):
 
    69 def latlonele_to_xyz(lat, lon, ele):
 
    62     s = sin(radians(lat))
 
    70     s = sin(radians(lat))
 
    63     c = cos(radians(lat))
 
    71     c = cos(radians(lat))
 
    64     r = ele + a * b / la.norm([s*a, c*b])
 
    72     r = ele + a * b / norm([s*a, c*b])
 
    65     lon = radians(lon)
 
    73     lon = radians(lon)
 
    66     return r * c * sin(lon), r * c * (-cos(lon)), r * s
 
    74     return r * c * sin(lon), r * c * (-cos(lon)), r * s
 
    67 
    75 
    68 
    76 
    69 def xyz_to_latlonele(x, y, z):
 
    77 def xyz_to_latlonele(x, y, z):
 
    70     r = la.norm([x, y, z])
 
    78     r = norm([x, y, z])
 
    71     if (r == 0):
 
    79     if (r == 0):
 
    72         return 0.0, 0.0, 0.0
 
    80         return 0.0, 0.0, 0.0
 
    73     lat = degrees(atan2(z, la.norm([x, y])))
 
    81     lat = degrees(atan2(z, norm([x, y])))
 
    74     lon = degrees(atan2(x, -y))
 
    82     lon = degrees(atan2(x, -y))
 
    75     ele = r * (1.0 - a * b / la.norm([a*z, b*x, b*y]))
 
    83     ele = r * (1.0 - a * b / norm([a*z, b*x, b*y]))
 
    76     return lat, lon, ele
 
    84     return lat, lon, ele
 
    77 
    85 
    78 
    86 
    79 ############################## main function #################################
 
    87 ############################## main function #################################
 
    80 tracks = []
 
    88 tracks = []
 
   151     dr = (dy - dx) / 2
 
   159     dr = (dy - dx) / 2
 
   152     xmax += dr
 
   160     xmax += dr
 
   153     xmin -= dr
 
   161     xmin -= dr
 
   154 
   162 
   155 from subprocess import Popen, PIPE
 
   163 from subprocess import Popen, PIPE
 
   156 gnuPlotCmd = ['gnuplot']
 
   164 plot = Popen ([gnuPlotCmd], stdin=PIPE, stdout=PIPE, stderr=PIPE)
 
   157 plot = Popen (gnuPlotCmd, stdin=PIPE, stdout=PIPE, stderr=PIPE)
 
       
   158 
   165 
   159 range = 'set xrange [%f:%f]\nset yrange [%f:%f]\n' % (xmin, xmax, ymin, ymax)
 
   166 range = 'set xrange [%f:%f]\nset yrange [%f:%f]\n' % (xmin, xmax, ymin, ymax)
 
   160 plot.stdin.write (range)
 
   167 plot.stdin.write (range)
 
   161 
   168 
   162 curves = ','.join ("'-' with linespoints ti '%s: %d punten'" % t for t in npoints)
 
   169 curves = ','.join ("'-' with linespoints ti '%s: %d punten'" % t for t in npoints)
reduce