GPS 활용기

실마리 2008. 4. 5. 20:49

2008. 4. 5. 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 파일 이외의 형식도 마지막 부분의 소스를 조금만 수정하면 사용할 수 있도록 만들어져 있습니다.

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

실마리의 홈페이지

GPS 활용기

+ Recent posts

티스토리툴바