티스토리 툴바

BLOG main image
분류 전체보기 (463)
Random (335)
컴퓨터 (37)
Favorites (20)
사진 (9)
아포리즘 (8)
자전거 (33)
여행 (21)

Twitter Updates

    follow me on Twitter
    «   2012/02   »
          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      
    162,105 Visitors up to today!
    Today 10 hit, Yesterday 22 hit
    rss
    2008/04/05 20:49

    지난번에 말씀드린 Holux M-241…

    시간 간격(1, 5, 10, 15, 30, 60, 120초) 혹은 지정 거리(50, 100, 150, 300, 500, 1000m)별로 위치를 저장하고 총 13만개정도의 위치를 저장할 수 있으니 36시간 400초(13만초) – 180일 13시간 20분(15600000초) 정도의 시간, 혹은 6500Km – 130000Km의 거리를 저장할 수 있는 셈입니다.

    BT747 프로그램을 사용하여 맥에서 로그 데이터를 받아 구글 어스에서 받는 것까지는 성공했지만 막상 전체 거리, 속도 등의 자료를 볼 수 없어서 이를 계산하는 프로그램을 짜 보았습니다. 프로그램은 XML의 형식을 가지는 KML 파일을 읽어들여 저장되어 있는 위도, 경도와 시간을 사용하여 이동거리 및 속도를 기록하고 이것을 기초로 간단한 그래프를 보여줍니다.

    프로그램은 GPS 자료의 위도, 경도 차이를 기준으로 거리를 계산합니다. 위도, 경도를 계산하는 방법은 Haversine 방법과 Vincenty 방법 이 알려져 있는데 Vincenty 방법은 지구를 타원으로 측정하여 1m 정도까지 정확하게 계산할 수 있는 모양입니다. 집에서 직장까지의 거리를 Vincenty 방법과 Haversine 방법으로 계산한 결과는 각각 9.212Km, 9.206Km로 대략 60m 정도 차이가 나지만 어차피 GPS가 크게 정확하지 않으리라 생각하고 Haversine 방법으로 거리를 계산했습니다.

    다음은 지난 금요일 자전거를 타고 출퇴근하면서 10초 간격으로 GPS를 기록한 것을 그래프로 나타낸 것입니다. 퇴근때는 GPS를 켜자마자 바로 달려서 초반 부분이 기록되어 있지 않습니다.

    사용자 삽입 이미지
    먼저 출근, 거리는 10.4Km를 달렸고 시간은 31분 정도가 걸렸으면 평속 20Km정도 입니다. 10초 이상 신호등에 걸린 부분에서는 속도가 0이 됩니다. 그리고 후반부에 속도가 상대적으로 감소되는데 이는 산중턱에 있는 직장 근처에서 시작되는 오르막길 때문입니다.

    사용자 삽입 이미지
    다음은 퇴근입니다. 코스가 약간 차이나지만 거리가 9Km로 초반부 GPS 신호를 잡지 못했기 때문으로 생각됩니다. 전박적으로 오전에 비해 속도가 빠른 편이고 경성대학교 부근에서 1번 신호에 잡힌것 말고는 계속 달렸음을 알수 있습니다.

    전체 프로그램의 소스입니다. 그래프는 rubyCocoa를 이용했기 때문에 Leopard이상의 OSX에서 실행해야 합니다. kml파일 목록을 인자로 실행시키면 기본적으로 GPS 기록이 1시간 이상 차이날 때마다 새로운 그래프를 만들어줍니다. KML 파일 이외의 형식도 마지막 부분의 소스를 조금만 수정하면 사용할 수 있도록 만들어져 있습니다.

    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
    
    require 'rexml/document'
    require 'time'
    require 'osx/cocoa'
    # require 'rubygems'
    # require 'units'
    
    # great resources about distance calculations on web
    # http://williams.best.vwh.net/avform.htm
    # http://www.movable-type.co.uk/scripts/latlong.html
    # http://www.movable-type.co.uk/scripts/latlong-vincenty.html
    # http://ajax.suaccess.org/rubyisms-in-rails/converting-between-degrees-and-radians/
    
    include Math
    include REXML
    
    SUMMARY_KO = { :dist=>"전체 거리 : %.2f Km", :time=>"시간 : %d분 %d초", :max_v=>"최고 속도 : %.2f Km/시", :avg_v=>"평균 속도 : %.2f Km/시" }
    SUMMARY_EN = { :dist=>"Total Distance : %.2f Km", :time=>"Time : %d:%d", :max_v=>"Max Velocity : %.2f Km/hr", :avg_v=>"Avg Velocity : %.2f Km/hr" }
    SUMMARY = SUMMARY_KO
    
    class Numeric
      def to_rad
        self*Math::PI/180
      end
      # add_unit_conversions(:angle => { :radians => 1, :degrees => Math::PI/180 })
      # add_unit_aliases(:angle => { :degrees => [:degree], :radians => [:radian] })
    end
    
    def distVincenty(lat1, lon1, lat2, lon2) 
      a, b = 6378137.0, 6356752.3142
      f = 1/298.257223563;                              # WGS-84 ellipsiod
      l = (lon2-lon1).to_rad
      u1 = atan((1-f) * tan(lat1.to_rad))
      u2 = atan((1-f) * tan(lat2.to_rad))
      sinU1, cosU1 = sin(u1), cos(u1)
      sinU2, cosU2 = sin(u2), cos(u2)
      
      lambda, lambdaP = l, 2*Math::PI
      iterLimit = 19
      while ((lambda-lambdaP).abs > 1e-12 and iterLimit>0) do
        sinLambda, cosLambda = sin(lambda), cos(lambda)
        sinSigma = sqrt((cosU2*sinLambda) ** 2 + (cosU1*sinU2-sinU1*cosU2*cosLambda) ** 2)
        return 0 if (sinSigma==0)                       # co-incident points
        cosSigma = sinU1*sinU2 + cosU1*cosU2*cosLambda
        sigma = atan2(sinSigma, cosSigma)
        sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma
        cosSqAlpha = 1 - sinAlpha ** 2
        cos2SigmaM = cosSigma - 2*sinU1*sinU2/cosSqAlpha
        # if (isNaN(cos2SigmaM)) cos2SigmaM = 0          # equatorial line: cosSqAlpha=0 (§6)
        c = f/16*cosSqAlpha*(4+f*(4-3*cosSqAlpha))
        lambdaP = lambda
        lambda = l + (1-c) * f * sinAlpha * \
          (sigma + c*sinSigma*(cos2SigmaM+c*cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)))
        iterLimit -= 1
      end
    
      return nil if (iterLimit==0)                      # formula failed to converge
    
      uSq = cosSqAlpha * (a*a - b*b) / (b*b)
      biga = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)))
      bigb = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)))
      deltaSigma = bigb*sinSigma*(cos2SigmaM+bigb/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM) - \
        bigb/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)))
      b*biga*(sigma-deltaSigma)
    end
    
    def distHaversine(lat1, lon1, lat2, lon2)
      # R = 6371
      dLat = (lat2-lat1).to_rad
      dLon = (lon2-lon1).to_rad
      a = sin(dLat/2) * sin(dLat/2) + cos(lat1.to_rad) * cos(lat2.to_rad) * sin(dLon/2)**2
      c = 2 * atan2(sqrt(a), sqrt(1-a)) 
      6371 * c * 1000;
    end
    
    def velocity(dist, t_int)
      (dist / 1000) / t_int * 3600
    end
    
    def gps_calc(doc, path, loc_f, h_f, t_f, cut_interval = 60 * 5)
      result = []
      rhash = { :data => [] }
      
      pre_lat, pre_lon, pre_h, pre_t = 0, 0, 0, 0
      start_time, total_dist, total_points, time_int, max_v = nil, 0, 0, 0, 0
    
      doc.elements.each(path) do |p|
        lat, lon = loc_f.call(p)
        h = h_f.call(p)
        t = t_f.call(p)
        
        start_time = t if start_time.nil?
        
        if pre_lat != 0 then
          time_int = t - pre_t
          if time_int <= cut_interval then
            dist = distHaversine(lat, lon, pre_lat, pre_lon)
            total_dist += dist
            total_points += 1
            v = velocity(dist, time_int)
            max_v = v if v > max_v
            rhash[:data] << [t - start_time, dist, v]
          end
        end
    
        if time_int > cut_interval then
          avg_v = velocity(total_dist, pre_t - start_time)
          rhash[:start_time] = start_time
          rhash[:total_points]= total_points
          rhash[:total_dist] = total_dist / 1000
          rhash[:time_int] = pre_t - start_time
          rhash[:avg_v] = avg_v
          rhash[:max_v] = max_v 
          result << rhash
          rhash = { :data => [] }      
          total_dist, total_points, max_v, time_int = 0, 0, 0, 0
          pre_lat, pre_lon, start_time = 0, 0, nil
        else
          pre_lat, pre_lon, pre_h, pre_t = lat, lon, h, t
        end    
      end
    
      if pre_lat != 0 then
        avg_v = velocity(total_dist, pre_t - start_time)
        rhash[:start_time] = start_time
        rhash[:total_points]= total_points
        rhash[:total_dist] = total_dist / 1000
        rhash[:time_int] = pre_t - start_time
        rhash[:avg_v] = avg_v
        rhash[:max_v] = max_v 
        result << rhash
      end
      
      result    
    end
    
    def mark_x(gps_data, value, para, emphasis=false)
      position = value * para[:graph_width] / gps_data[:total_dist] + para[:margin]
    
      vstring = value.integer? ? value.to_s : "%.2f" % value
      label = OSX::NSString.alloc.initWithString(vstring)
      font_dict = emphasis ? para[:em_font_dict] : para[:font_dict]
      size = label.sizeWithAttributes(font_dict)
      label.drawAtPoint_withAttributes([position-size.width/2, para[:zero].y-size.height-para[:font_margin]], font_dict)
      OSX::NSRectFill([position, para[:zero].y-para[:font_margin], 1, para[:font_margin]])
    end
    
    def mark_y(gps_data, value, para, emphasis=false)
      position = value * para[:graph_height] / (gps_data[:max_v] * 1) + para[:margin]
      
      vstring = value.integer? ? value.to_s : "%.2f" % value
      label = OSX::NSString.alloc.initWithString(vstring)
      font_dict = emphasis ? para[:em_font_dict] : para[:font_dict]
      size = label.sizeWithAttributes(font_dict)
      label.drawAtPoint_withAttributes([para[:margin]-size.width-para[:font_margin], position-size.height/2], font_dict)
      OSX::NSRectFill([para[:zero].x-para[:font_margin], position, para[:font_margin], 1])
    end
    
    def make_graph(gps_data, para={})
      default = { :width=>700, :height=>500, :mark_int=>60*10, :format=>OSX::NSPNGFileType, :margin=>50, :font_size=>12, :font_margin=>3 }
      dist_scales = [[1, 2, 5, 10, 20, 50, 100], [5, 10, 20, 50, 100, 200, 1000]]
      
      default.update(para)
    
      canvas = OSX::NSBitmapImageRep.alloc.initWithBitmapDataPlanes_pixelsWide_pixelsHigh_bitsPerSample_samplesPerPixel_hasAlpha_isPlanar_colorSpaceName_bytesPerRow_bitsPerPixel(nil, default[:width], default[:height], 8, 4, true, false, OSX::NSDeviceRGBColorSpace, 0, 0)
      context = OSX::NSGraphicsContext.graphicsContextWithBitmapImageRep(canvas)
      OSX::NSGraphicsContext.setCurrentContext(context)
    
      # font
      white = OSX::NSColor.whiteColor
      white.set
      yellow = OSX::NSColor.yellowColor
      
      font = OSX::NSFont.fontWithName_size('Helvetica', default[:font_size])
      font_dict = OSX::NSMutableDictionary.alloc.init
      font_dict.setObject_forKey(font, OSX::NSFontAttributeName)
      font_dict.setObject_forKey(white, OSX::NSForegroundColorAttributeName)
      default[:font_dict] = font_dict
    
      em_font = OSX::NSFont.boldSystemFontOfSize(default[:font_size]+2)
      em_font_dict = OSX::NSMutableDictionary.alloc.init
      em_font_dict.setObject_forKey(em_font, OSX::NSFontAttributeName)
      em_font_dict.setObject_forKey(OSX::NSColor.yellowColor, OSX::NSForegroundColorAttributeName)
      default[:em_font_dict] = em_font_dict  
      
      # background gradient
      gradient = OSX::NSGradient.alloc.initWithStartingColor_endingColor(OSX::NSColor.blueColor, OSX::NSColor.blackColor)
      gradient.drawInRect_angle([0, 0, default[:width], default[:height]], 90)
      
      # lines   
      default[:zero] = OSX::NSMakePoint(default[:margin], default[:margin])
      default[:x_end] = OSX::NSMakePoint(default[:width]-default[:margin], default[:margin])
      default[:y_end] = OSX::NSMakePoint(default[:margin], default[:height]-default[:margin])
      path = OSX::NSBezierPath.bezierPath
      path.moveToPoint(default[:x_end])
      path.lineToPoint(default[:zero])
      path.lineToPoint(default[:y_end])
      path.stroke
      
      # labels
      # dist_label = OSX::NSString.alloc.initWithString('Km')
      # dist_size = dist_label.sizeWithAttributes(font_dict)
      # dist_label.drawAtPoint_withAttributes([default[:x_end].x+default[:font_margin], default[:x_end].y-dist_size.height-default[:font_margin]], font_dict)
      # velo_label = OSX::NSString.alloc.initWithString('Km/Hr')
      # velo_size = velo_label.sizeWithAttributes(font_dict)
      # velo_label.drawAtPoint_withAttributes([default[:y_end].x-velo_size.width-default[:font_margin], default[:y_end].y+default[:font_margin]], font_dict)
    
      # 
      default[:graph_width] = default[:width] - default[:margin] * 2
      default[:graph_height] = default[:height] - default[:margin] * 2
      
      # determine x scale units
      scale_index = 0
      (dist_scales[1].size-1).downto(0) { |i| scale_index=i; break if gps_data[:total_dist] > dist_scales[1][i] }
      scale = dist_scales[0][scale_index]
      
      # draw x legends
      # mark_x(gps_data, gps_data[:total_dist], default, true)
      dist = scale
      while (dist < gps_data[:total_dist]) do 
        mark_x(gps_data, dist, default)
        dist += scale
      end
      
      # mark_y(gps_data, gps_data[:max_v], default)
      velo = 10
      while (velo < gps_data[:max_v]) do
        mark_y(gps_data, velo, default)
        velo += 10
      end
      
      # draw graph
      graph = nil
      total = 0
      gps_data[:data].each do |time_int, dist, velo|
        total += dist
        point_x = (total / 1000) * default[:graph_width] / gps_data[:total_dist] + default[:margin]
        point_y = velo * default[:graph_height] / (gps_data[:max_v] * 1) + default[:margin]
        
        if graph.nil? then
          graph = OSX::NSBezierPath.bezierPath
          graph.moveToPoint(OSX::NSMakePoint(point_x, point_y))
        else
          graph.lineToPoint(OSX::NSMakePoint(point_x, point_y))
        end
        
        if (default[:mark_int] != 0) and (time_int % default[:mark_int] == 0) then
          mark_label = OSX::NSString.alloc.initWithString("%d:%02d" % [time_int / 60, time_int % 60])
          mark_size = mark_label.sizeWithAttributes(default[:font_dict])
          mark_label.drawAtPoint_withAttributes([point_x-mark_size.width/2, point_y+default[:font_margin]], default[:font_dict])
          mark = OSX::NSBezierPath.bezierPathWithOvalInRect([point_x-2, point_y-2, 4, 4])
          yellow.set
          mark.stroke
          white.set      
        end
      end
      graph.stroke
      
      # summary
      sum_y = default[:height] - default[:margin] / 2
      dist_label = OSX::NSString.alloc.initWithString(SUMMARY[:dist] % gps_data[:total_dist])
      dist_size = dist_label.sizeWithAttributes(default[:em_font_dict])
      dist_label.drawAtPoint_withAttributes([default[:width]-default[:margin]/2-dist_size.width, sum_y-dist_size.height], default[:em_font_dict])
      sum_y -= dist_size.height + default[:font_margin]
      time_label = OSX::NSString.alloc.initWithString(SUMMARY[:time] % [gps_data[:time_int]/60, gps_data[:time_int]%60])
      time_size = time_label.sizeWithAttributes(default[:em_font_dict])
      time_label.drawAtPoint_withAttributes([default[:width]-default[:margin]/2-time_size.width, sum_y-time_size.height], default[:em_font_dict])
      sum_y -= time_size.height + default[:font_margin]
      maxv_label = OSX::NSString.alloc.initWithString(SUMMARY[:max_v] % gps_data[:max_v])
      maxv_size = maxv_label.sizeWithAttributes(default[:em_font_dict])
      maxv_label.drawAtPoint_withAttributes([default[:width]-default[:margin]/2-maxv_size.width, sum_y-maxv_size.height], default[:em_font_dict])
      sum_y -= maxv_size.height + default[:font_margin]
      avgv_label = OSX::NSString.alloc.initWithString(SUMMARY[:avg_v] % gps_data[:avg_v])
      avgv_size = avgv_label.sizeWithAttributes(default[:em_font_dict])
      avgv_label.drawAtPoint_withAttributes([default[:width]-default[:margin]/2-avgv_size.width, sum_y-avgv_size.height], default[:em_font_dict])
    
      
      canvas.representationUsingType_properties(default[:format], nil).rubyString  
    end
    
    if ARGV.size > 0 then
      ARGV.each do |a|
        doc = Document.new(open(a))
        
        loc_f = proc { |e| e.elements['Point'].elements['coordinates'].text.split(',').collect { |s| s.to_f}[0..2] }
        h_f = proc { |e| e.elements['Point'].elements['coordinates'].text.split(',').collect { |s| s.to_f}[2] }
        t_f = proc { |e| Time.parse(e.elements['TimeStamp'].elements['when'].text).localtime }
        
        r = gps_calc(doc, 'kml/Document/Folder/Folder/Placemark', loc_f, h_f, t_f, 60*60)
        r.each do |result|
          open("graph_#{result[:start_time].strftime('%Y%m%d%H%M')}.png", 'wb') { |f| f.write(make_graph(result))} 
        end
      end
    end