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

その2('$GPRMC'編)かきました。

手順

1. NMEAログ取得

• 1. GPSロガーから。
• 1. GPSジュールを、パソコンにシリアル接続、TeraTermなどで、テキスト保存。
• 1. 'GPSロガー　サンプルデータ　ダウンロード' などで、グーぐる先生にきいてみる

1. 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]) // 最初の位置からの相対座標

4. library(rgl)というパッケージで、plot3dしてみる。(マウスでぐりぐりできる)
　* 実行例3: > plot3d(rt$Longitude-rt$Longitude[1], rt$Latitude-rt$Latitude[1], rt$Altitude)

---

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を使用)

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")

例

こちらを無断でお借りしました。
> rt = geo_position_from_nmealog("http://www.oiccam.com/reno/gps/tahoe/nmea/2005-09-25-north-tahoe_nmea.txt", is0000=FALSE)