stud

study spacejunk
Log | Files | Refs | LICENSE

measure.py (8550B) - Raw


      1 #!/usr/bin/env python3
      2 """
      3 Execute like this:
      4     ./measure.py  | column -t -s $'\t'
      5 """
      6 
      7 from collections import namedtuple
      8 from decimal import Decimal as Dec
      9 from math import sin, cos, pi
     10 from shapely.geometry import LineString, asPolygon, Point as sPoint
     11 import numpy as np
     12 
     13 
     14 class Deg(namedtuple('Deg', ['deg', 'mm', 'ss'])):
     15     def __str__(self):
     16         return "%03d-%02d-%04.1f" % (self.deg, self.mm, self.ss)
     17 
     18 def hms(deg):
     19     assert isinstance(deg, Dec)
     20     pdeg, pmm = divmod(deg, 1)
     21     pmm = pmm * Dec(60)
     22     pmm, pss = divmod(pmm, 1)
     23     pss = pss * Dec(60)
     24     return Deg(pdeg, pmm, pss)
     25 
     26 def normalize(ang):
     27     while ang > 180:
     28         ang -= 360
     29     while ang <= -180:
     30         ang += 360
     31     return ang
     32 
     33 def guess(inp):
     34     if isinstance(inp, str) and '-' in inp:
     35         deg, mm, ss = inp.split('-')
     36         ddeg, dmm, dss = Dec(deg), Dec(mm), Dec(ss)
     37         return ddeg + dmm/60 + dss/3600
     38     else:
     39         return Dec(instr)
     40 
     41 class Point(namedtuple('Point', ['acadx', 'acady'])):
     42     @property
     43     def lksx(self):
     44         return self.acady
     45     @property
     46     def lksy(self):
     47         return self.acadx
     48 
     49 class Vertex:
     50     def __init__(self, point, length, angle, dirang=Dec(), coords = Point(Dec(), Dec())):
     51         self.point = point
     52         self.len = length
     53         self.ang = angle
     54         self.dirang = dirang
     55         self.coords = coords
     56         self.dx, self.dy = Dec(), Dec()
     57 
     58     @property
     59     def xy(self):
     60         """xy returns a tuple of lksx and lksy coordinates"""
     61         return np.array([float(self.coords.lksx), float(self.coords.lksy)])
     62 
     63 def heptagon(d1):
     64     angles = np.linspace(0, 2*pi, num=8)
     65     R = float(D1)/2/sin(pi/7)
     66     heptagon_xy = (np.array([np.cos(angles), np.sin(angles)])*R).T
     67     return asPolygon(heptagon_xy)
     68 
     69 juosta = namedtuple('juosta', ['plotis', 'kryptis', 'dashes', 'spalva'])
     70 kelias = namedtuple('kelias', ['virsunes', 'plotis', 'kat', 'dashes', 'spalva', 'juostos'])
     71 
     72 # Kategorijos
     73 KAT0, KAT1, KAT2, KAT3, KAT4 = range(5,0,-1)
     74 
     75 A = Dec('6.094')
     76 B = Dec('-2.923')
     77 C = Dec('-13.462')
     78 N = Dec('9.512')
     79 
     80 L1 = Dec('16.321')
     81 # === Kelias A-05 ===
     82 L2 = Dec('9.109')
     83 L3 = Dec('4.819')
     84 # === Kelias A-08 ===
     85 L4 = Dec('2.675')
     86 L5 = Dec('2.059')
     87 L6 = Dec('1.262')
     88 L7 = Dec('4.170')
     89 L8 = Dec('6.005')
     90 L9 = Dec('6.453')
     91 # === Griovys G-11 ===
     92 L10 = Dec('4.882')
     93 L11 = Dec('3.305')
     94 L12 = Dec('2.210')
     95 L13 = Dec('4.381')
     96 
     97 A03_plotis = Dec('17.401') + A
     98 A05_plotis = Dec('13.705') + B
     99 A08_plotis = Dec('29.006') + C
    100 G11_plotis = Dec('14.776') + N
    101 
    102 # Directional coords + angle
    103 X11 = Dec('6091968.055')
    104 Y11 = Dec('485944.146')
    105 A11_2 = guess('70-16-17')
    106 
    107 vertices = [
    108   #  point   len             angle               dirangle  coords
    109   Vertex(11, Dec('164.126'), guess('103-03-03'), A11_2,    Point(X11, Y11)),
    110   Vertex(2,  Dec('149.851'), guess('218-27-42')),
    111   Vertex(19, Dec('82.384' ), guess('211-44-30')),
    112   Vertex(3,  Dec('259.022'), guess('67-26-49' )),
    113   Vertex(24, Dec('319.331'), guess('67-33-06' )),
    114   Vertex(12, Dec('74.764' ), guess('279-03-59')),
    115   Vertex(13, Dec('81.640' ), guess('278-54-55')),
    116   Vertex(14, Dec('31.888' ), guess('119-27-45')),
    117   Vertex(15, Dec('84.073' ), guess('160-50-28')),
    118   Vertex(16, Dec('70.072' ), guess('207-42-31')),
    119   Vertex(17, Dec('73.378' ), guess('206-18-01')),
    120   Vertex(10, Dec('66.625' ), guess('90-55-10' )),
    121   Vertex(18, Dec('97.003' ), guess('100-18-10')),
    122   Vertex(9,  Dec('121.003'), guess('148-30-56')),
    123   Vertex(8,  Dec('131.915'), guess('285-20-57')),
    124   Vertex(23, Dec('102.086'), guess('29-44-22' )),
    125   Vertex(22, Dec('158.324'), guess('276-33-49')),
    126   Vertex(7,  Dec('72.157' ), guess('82-07-47' )),
    127   Vertex(6,  Dec('107.938'), guess('104-15-46')),
    128   Vertex(21, Dec('104.082'), guess('234-17-37')),
    129   Vertex(5,  Dec('154.332'), guess('283-30-57')),
    130   Vertex(20, Dec('68.972' ), guess('152-15-58')),
    131   Vertex(1,  Dec('151.531'), guess('101-20-01')),
    132   Vertex(4,  Dec('179.336'), guess('150-15-41')),
    133 ]
    134 
    135 for i, v in enumerate(vertices[1:]):
    136     prev = vertices[i]
    137     v.dirang = prev.dirang + 180 - v.ang
    138     v.dx = Dec(float(prev.len) * cos(float(prev.dirang) * pi/180))
    139     v.dy = Dec(float(prev.len) * sin(float(prev.dirang) * pi/180))
    140     v.coords = Point(prev.coords.acadx + v.dx, prev.coords.acady + v.dy)
    141 
    142 angle_sum = Dec(0)
    143 for v in vertices:
    144     angle_sum += v.ang
    145 theoretical_angle_sum = Dec(int((len(vertices)-2)*180))
    146 
    147 # 9-kampio krastine D1
    148 D1 = Dec('174.667') + C
    149 # Daugiakampio pasukimo kampas (K1)
    150 K1 = Dec('13.147') + B
    151 # Atstumas iki tikrosios uzliejimo zonos (A1) (0.001 tikslumu)
    152 A1 = Dec('67.536') + B
    153 
    154 circle_radius = float(D1)/2/sin(pi/7)-float(A1)
    155 heptagon_area = heptagon(float(D1)).area
    156 circle_area = sPoint(0,0).buffer(circle_radius).area
    157 
    158 # Points is vertice map by id
    159 Points = {}
    160 for v in vertices:
    161     Points[v.point] = v
    162 
    163 CONTINUOUS = (1,0)
    164 DASHDOTX2 = (10,3,2,3)
    165 DASHED = (100,20)
    166 
    167 keliai = {
    168     'A-08': kelias(
    169         virsunes=[1,2,3],
    170         plotis=A08_plotis,
    171         kat=KAT1,
    172         dashes=DASHDOTX2,
    173         spalva='xkcd:red',
    174         juostos=(
    175             juosta(L6+L5+L4, 'right', DASHED,     'xkcd:lightgreen'),
    176             juosta(L6+L5,    'right', DASHED,     'xkcd:lightgreen'),
    177             juosta(L6,       'right', CONTINUOUS, 'xkcd:black'),
    178             juosta(L7,       'left',  CONTINUOUS, 'xkcd:black'),
    179             juosta(L7+L8,    'left',  DASHED,     'xkcd:lightgreen'),
    180             juosta(L7+L8+L9, 'left',  DASHED,     'xkcd:lightgreen'),
    181         ),
    182     ),
    183     'A-05': kelias(
    184         virsunes=[4,5,6,7,8,9,10],
    185         plotis=A05_plotis,
    186         kat=KAT2,
    187         dashes=DASHDOTX2,
    188         spalva='xkcd:red',
    189         juostos=(
    190             juosta(L3, 'right', CONTINUOUS, 'xkcd:brown'),
    191             juosta(L2, 'left',  CONTINUOUS, 'xkcd:brown'),
    192         ),
    193     ),
    194     'A-03': kelias(
    195         virsunes=[11,12,13,14,15,16,17,18],
    196         plotis=A03_plotis,
    197         kat=KAT3,
    198         dashes=CONTINUOUS,
    199         spalva='xkcd:magenta',
    200         juostos=(
    201             juosta(L1, 'right', DASHED, 'xkcd:magenta'),
    202             juosta(0,   'left', DASHED, 'xkcd:white'),
    203         ),
    204     ),
    205     'G-11': kelias(
    206         virsunes=[19,20,21,22,23,24],
    207         plotis=G11_plotis,
    208         kat=KAT4,
    209         dashes=CONTINUOUS,
    210         spalva='xkcd:red',
    211         juostos=(
    212             juosta(L10+L11, 'right', CONTINUOUS, 'xkcd:blue'),
    213             juosta(L11,     'right', CONTINUOUS, 'xkcd:lightblue'),
    214             juosta(L12,     'left',  CONTINUOUS, 'xkcd:lightblue'),
    215             juosta(L12+L13, 'left',  CONTINUOUS, 'xkcd:blue'),
    216         ),
    217     ),
    218 }
    219 
    220 keliu_ilgiai = {}
    221 for id, kelias in keliai.items():
    222     keliu_ilgiai[id] = LineString([Points[i].xy for i in kelias.virsunes]).length
    223 
    224 if __name__ == '__main__':
    225     print("tšk. nr.\tišmatuotas kampas\tdirekcinis kampas\tilgis\tdx\tdy\tx\ty")
    226     for i, v in enumerate(vertices):
    227         print("\t".join([
    228             "%d" % v.point,
    229             "%s" % str(hms(v.ang)),
    230             "%s" % str(hms(v.dirang)),
    231             "%.3f" % v.len,
    232             "%.3f" % v.dx,
    233             "%.3f" % v.dy,
    234             "%.3f" % v.coords.acadx,
    235             "%.3f" % v.coords.acady,
    236         ]))
    237 
    238         #acad coords for drawing
    239         """
    240         nxt = vertices[0 if i == len(vertices) - 1 else i+1]
    241         pts = "%d-%d" % (v.point, nxt.point)
    242         draw = "@%.3f<%.4f" % (v.len, normalize(90 - v.dirang))
    243         print("%5s: %19s  acadcoords:(%.3f,%.3f)" % \
    244                 (pts, draw, v.coords.acadx, v.coords.acady))
    245         """
    246 
    247     # debugging & while drawing
    248     ("""
    249     Kelio A-03 plotis = 17.401 + A = %.3f""" % A03_plotis + """
    250     Kelio A-05 plotis = 13.705 + B = %.3f""" % A05_plotis + """
    251     Kelio A-08 plotis = 29.006 + C = %.3f""" % A08_plotis + """
    252     Griovio G-11 plotis = 14.776 + N = %.3f""" % G11_plotis + """
    253 
    254     Prognozuojamo uzliejimo zona, tai taisyklingas 9-kampis
    255     9-kampio krastine D1 = %.3f""" % D1 + """
    256     Daugiakampio pasukimo kampas (K1) (0.0001 laipsnio tikslumu)
    257     K1 = %.4f""" % K1 + """
    258 
    259     Tikroji uzliejimo zona, tai taisyklingas apskritimas, kurio centras TURI SUTAPTI su daugiakampio centru.
    260     Atstumas iki tikrosios uzliejimo zonos (A1) (0.001 tikslumu)
    261     A1 = %.3f""" % A1 + """
    262 
    263     A-05:
    264     x(l) = %.3f""" % (A05_plotis*L3/(L2+L3)) + """
    265     x(r) = %.3f""" % (A05_plotis*L2/(L2+L3)) + """
    266 
    267     A-08:
    268     x(l) = %.3f""" % (A08_plotis*(L7+L8+L9)/(L7+L8+L9+L6+L5+L4)) + """
    269     x(r) = %.3f""" % (A08_plotis*(L6+L5+L4)/(L7+L8+L9+L6+L5+L4)) + """
    270 
    271     G-11:
    272     x(l) = %.3f""" % (G11_plotis*(L12+L13)/(L10+L11+L12+L13)) + """
    273     x(r) = %.3f""" % (G11_plotis*(L10+L11)/(L10+L11+L12+L13)) + """
    274     """)