Rで GPSのNMEA形式のデータであそぶ(その1)

[その2('$GPRMC'編)](http://qiita.com/mt08/items/9164fa3dedb984754cf6)かきました。

## 手順
1. NMEAログ取得
 * 1. GPSロガーから。
 * 2. GPSジュールを、パソコンにシリアル接続、TeraTermなどで、テキスト保存。
 * 3. 'GPSロガー　サンプルデータ　ダウンロード' などで、グーぐる先生にきいてみる
2. Rに食わせる。下の`geo_position_from_nmealog`を使用。最後に改行がないと怒られたりする。うまく対応する。やってることは、$GPGGAをフィルタして、latとlngをうまく変換して返している。
 * 例: `> rt = geo_position_from_nmealog("C:/gps/gps1234.log")`  

3.　軽く plotしてみる
　* 実行例1: `> plot(rt$Latitude, rt$Longitude)`
　* 実行例2: `> plot(rt$Longitude-rt$Longitude[1], rt$Latitude-rt$Latitude[1])` // 最初の位置からの相対座標
 ![plot.png](https://qiita-image-store.s3.amazonaws.com/0/75594/34be6f10-381f-bf5f-b1d2-a32f112511e3.png)
4. `library(rgl)`というパッケージで、plot3dしてみる。(マウスでぐりぐりできる)
　* 実行例3: `> plot3d(rt$Longitude-rt$Longitude[1], rt$Latitude-rt$Latitude[1], rt$Altitude)`
 ![plot3d.png](https://qiita-image-store.s3.amazonaws.com/0/75594/a088c9a3-52b6-fd90-46d1-2d1e779cca99.png)

***
```r:geo_position_from_nmealog
geo_position_from_nmealog <- function(log_filename, is0000 = TRUE) {
     # ファイルを読む
     lines = readLines(log_filename)
     # フィルタ: $GPGGAのラインを取得
	 rexp_str = ifelse(is0000 == TRUE, "^\\$GPGGA,.*0000\\*[0-9A-F][0-9A-F]$", "^\\$GPGGA,.*\\*[0-9A-F][0-9A-F]$")
     lines_GPGGA = lines[grep(rexp_str, lines)]
     # tableへ
     tbl_GPGGA = read.table( text =lines_GPGGA, sep=",", header=FALSE)
     # カラムに名前をつける
     colnames(tbl_GPGGA) <- (c("SentenceId","TimeUTC", "Latitude", "LatitudeHemisphere", "Longitude", "LongitudeHemisphere", "PositionFixQuality", "NumOfSatelites", "HDOP", "Altitude", "AltitudeUnit", "GeoidalHeight","GeoidalHeightUnit", "TimeSeinceLastDGPSUpdate", "_CheckSum"))
     # フィルタ: 取得できたデータのみ
     tbl_GPGGA =tbl_GPGGA [!(tbl_GPGGA $PositionFixQuality == 0), ]

     # 緯度Latitudeと経度Longitudeを変換
     rt <- data.frame(
          Latitude  = ifelse(tbl_GPGGA$LatitudeHemisphere  == "N",  1, -1) * (as.integer(tbl_GPGGA$Latitude /100) +  (tbl_GPGGA$Latitude /100 - as.integer(tbl_GPGGA$Latitude /100)) * 100/60),
          Longitude = ifelse(tbl_GPGGA$LongitudeHemisphere == "E",  1, -1) * (as.integer(tbl_GPGGA$Longitude/100) +  (tbl_GPGGA$Longitude/100 - as.integer(tbl_GPGGA$Longitude/100)) * 100/60), 
          Altitude = tbl_GPGGA$Altitude
     )
     return(rt)
}
```


## Googleマップに投影 (`RgoogleMaps`を使用)

```r:test.R
install.packages("RgoogleMaps") # 入れてなければインストール
library(RgoogleMaps)

setwd("C:/gps/") #作業フォルダ指定
rt = geo_position_from_nmealog("gps1234.log")

lat <- as.vector(rt[,1])
lon <- as.vector(rt[,2])

# センター位置。　maxとminの間
center = c ((max(lat)+min(lat))/2, (max(lon)+min(lon))/2)
# zoom値計算: サンプルコードから... 自分で、zoom <- 13とかでも、おｋ
zoom <- min(MaxZoom(range(lat), range(lon)))

# グーぐる先生から地図画像を頂戴する。 日本語でおk
MyMap <- GetMap(center=center, zoom=zoom, size=c(640, 640), destfile = "./dest.png")

#色指定。
col <- rep(3, length(lat))

# プロット
tmp <- PlotOnStaticMap(
	MyMap,
	lat = lat,
	lon = lon,
	destfile = "./dest.png",
	cex=1.0,
	pch=20,
	col=col,
	add=FALSE,
	NEWMAP = FALSE
);

# PNG 保存
dev.print(bmp, width = 640,height=640, file = "./output.png")

```
## 例
[こちら](http://www.oiccam.com/reno/gps/tahoe/nmea/2005-09-25-north-tahoe_nmea.txt)を無断でお借りしました。
`> rt = geo_position_from_nmealog("http://www.oiccam.com/reno/gps/tahoe/nmea/2005-09-25-north-tahoe_nmea.txt", is0000=FALSE)
`
![output.png](https://qiita-image-store.s3.amazonaws.com/0/75594/0ed649e9-89b5-95fe-5944-a6d28cf7e239.png)