Rで GPSのNMEA形式のデータであそぶ(その2) '$GPRMC'編 : GPSの速度と方角の情報から矢印を描く

$GPRMCではじまる行では、緯度(Latitude)、経度(Longitude)と速度(knot)、進行方向(TrackAngle)が取得できる。ので、それらを元に、矢印を表示、描いてみよう。

ちょっとその前に、まずは、

GPSのデータチェック

grepのフィルタだけでは、抜けてくるのもあるので、チェックサムの計算して、はじいたほうがいいので、コードを書きました。
Rのお作法を知らないので、うごきゃーいいとなってるため、コメントなどで、ご指導いただけると助かります。
やってることは、左端1文字目は"$" 右端３文字は"XX"となるとして、2文字目から(length(line)-3)文字目まで、ASCIIコードの値でXORしていき、結果を、""後の2文字のCheckSumの値と比較する、です。
bitXorのために、"bitops"パッケージを入れてください。
> install.packages("bitops")

GPSの行データチェック

library(bitops)

is_valid_gps_nmea_checksum <-function(gpsstr) {
     len <- nchar(gpsstr)
     result = strtoi(substr(gpsstr, len-1,len), 16L)

     len <- len - 3
     tbl<-rep(0, len)

     for (k in 1:len) {
       ch = substr(gpsstr, k,k)
       tbl[k] = strtoi(charToRaw(ch), 16L)
     }

     chksum <- tbl[2]
     for (k in 3:len) {
        chksum <- bitXor(chksum, tbl[k])
     }
     return (chksum == result )
}

get_gps_valid_lines = function(lines) {
     valid <- 0
     i <- 0
     for (k in 1:length(lines)) { 
          if(is_valid_gps_nmea_checksum(lines[k])) {
               i<-i+1
               valid[i] = k
          }
     }
     return(valid)
}

データの準備

今回も、前回使用させていただいたデータをお借りします。
(http://www.oiccam.com/reno/gps/tahoe/nmea/2005-09-25-north-tahoe_nmea.txt)
なんども使う場合は、ローカルにダウンロードしといたほうがよろしいかと。

データrt0の準備

# 今回も`RgoogleMaps`にお世話になります。
library(RgoogleMaps)
library(bitops)

# 作業フォルダへ移動
setwd("C:/R_de_gpslog_asobu/")

# 1. データ読み込み
lines = readLines("2005-09-25-north-tahoe_nmea.txt")

# 2. "GPRMC"でフィルタ. で、その数を表示
lines_GPRMC <- lines[grep("GPRMC", lines)];  length(lines_GPRMC)

# 3. CheckSumでフィルタ. で、その数を表示(減ったかな？)
lines_GPRMC <- lines_GPRMC[get_gps_valid_lines(lines_GPRMC)] ; length(lines_GPRMC)
 
# 4. テーブルに変換
tbl_GPRMC = read.table( text =lines_GPRMC , sep=",", header=FALSE)

# 5. 行名をつける
colnames(tbl_GPRMC) <- (c("SentenceId", "TimeUTC", "Status",  "Latitude", "LatitudeHemisphere", "Longitude", "LongitudeHemisphere", "SpeedInKnots", "TrackAngle", "Date", "MagneticVariation", "MV", "_CheckSum"))

# 6. 捕捉されてるデータのみ
tbl_GPRMC = tbl_GPRMC[tbl_GPRMC$Status == "A", ]

# 7. tbl_GPRMCから、rt0 (時間、時速km、 時速mile、緯度、経度、方角)のデータを作成

rt0 <- data.frame(
	time = as.integer(tbl_GPRMC$Date/10000)*86400 + as.integer(tbl_GPRMC$TimeUTC / 10000) *3600 + (as.integer(tbl_GPRMC$TimeUTC / 100) - as.integer(tbl_GPRMC$TimeUTC / 10000)*100) * 60 +  as.integer(tbl_GPRMC$TimeUTC) - as.integer(tbl_GPRMC$TimeUTC / 100)*100 ,
	speed_kmph = tbl_GPRMC$SpeedInKnots *  1.85200 , 
	speed_miph =  tbl_GPRMC$SpeedInKnots * 1.15077945,
	Latitude  = ifelse(tbl_GPRMC$LatitudeHemisphere  == "N",  1, -1) * (as.integer(tbl_GPRMC$Latitude /100) +  (tbl_GPRMC$Latitude /100 - as.integer(tbl_GPRMC$Latitude /100)) * 100/60),
	Longitude = ifelse(tbl_GPRMC$LongitudeHemisphere == "E",  1, -1) * (as.integer(tbl_GPRMC$Longitude/100) +  (tbl_GPRMC$Longitude/100 - as.integer(tbl_GPRMC$Longitude/100)) * 100/60),
	TrackAngle = tbl_GPRMC$TrackAngle 
)

# 時間の単位秒にして、日付*86400 + 時*3600 + 分*60 + 秒
# 速度は、グーぐる先生によると
# 1 knot = 1.15077945 mile per hour`
# 1 knot = 1.85200 kilometer per hour`
# 緯度経度変換は前回のをこぴぺ

RgoogleMapsでプロット

前回の復習(PlotOnStaticMap)

前回と同様のプロット

PlotOnStaticMap使用

rt <- rt0
lat <- as.vector(rt$Latitude)
lon <- as.vector(rt$Longitude)
center = c ((max(lat)+min(lat))/2, (max(lon)+min(lon))/2)
zoom <- min(MaxZoom(range(lat), range(lon)))

MyMap <- GetMap(center=center, zoom=zoom, size=c(640, 640), destfile = "./dest.png", maptype = 'terrain',  hl = "ja")

tmp <- PlotOnStaticMap(MyMap,
                       lat = lat,
                       lon = lon,
                       destfile = "./dest.png",
                       cex=1.0,
                       pch=20,
                       col='red',
                       add=FALSE,
                       NEWMAP = FALSE
                      )

dev.print(png, file = "./result1.png",width = 640, height=640)

矢印でPlotArrowsOnStaticMap (その1)

PlotArrowsOnStaticMap で、Arrowが描ける。
パラメータで (lat0, lon0) から、(lat1, lon1)に向かって矢印がかけるので、lat1, lon1を速度と方角からうまく計算する。

PlotArrowsOnStaticMapその1

rt <- rt0

#
lat <- as.vector(rt$Latitude)
lon <- as.vector(rt$Longitude)
center = c ((max(lat)+min(lat))/2, (max(lon)+min(lon))/2)
zoom <- min(MaxZoom(range(lat), range(lon)))

# lat1, lon1の計算
arrow_size = 0.007 #矢の長さ
arrow_length = 0.04 #矢の傘の大きさ
lat1 <- cos(rt$TrackAngle *pi/180)*arrow_size *rt$speed_miph/mean(rt$speed_miph)
lon1 <- sin(rt$TrackAngle *pi/180)*arrow_size *rt$speed_miph/mean(rt$speed_miph)
lat1 <- lat + lat1
lon1 <- lon + lon1


# maptype: roadmap, satellite, hybrid, terrain
MyMap <- GetMap(center=center, zoom=zoom, size=c(640, 640), destfile = "./dest.png", maptype = 'terrain',  hl = "ja")

tmp <-PlotArrowsOnStaticMap(MyMap,
                            lat0 = lat,
                            lon0 = lon,
                            lat1 = lat1,
                            lon1 = lon1,
                            cex=1.0,
                            pch=20,
                            col='red',
                            add=FALSE,
                            length = arrow_length 
                           )

dev.print(png, file = "./result2.png",width = 640, height=640)

！！！
データが多すぎる...

矢印でPlotArrowsOnStaticMap (その2)

rtを作る際、データを間引きする。

データ間引き

step <- 15 # 15毎にひとつ。
rt <- rt0[0:(length(rt0[,1])/step) *step +1, ]

# 以下、その1と同じ

# やってることは、
# 0 から、データ数/step までのVector作成して、stepかけて、+1
# 例:データ数=5236, step=15
# 5236/15 = 349.0667..
# 0:349 :=  0 1 2 3 4 ... 348 349
# (0:349) * 15 :=  0 15 30 45  ... 5220 5235  (15毎)
# (0:349) * 15 +1 := 1 16 31 46 ... 5221 5236 
# 　↑で列フィルタ　data[(0:349) * 15 +1, ]　:= data[1, ], data[16,] ... data[5236, ]