GPSデータの解析・処理

Google Map 上に描く

fn <- "new"
ne <- decodeGPS(fn)
genMap(fn, ne, size = c(640, 640))

使う関数は以下の2つ

genMap <- function(fn, ne, col = "red", lwd = 3, maptype = "roadmap",
                   size = NULL, zoom = NULL, format = "png") {
    library(RgoogleMaps)
    N <- ne$N
    E <- ne$E
    lat <- range(N)
    lon <- range(E)
    center <- c(mean(lat), mean(lon))
    if (is.null(size)) {
        size <- c(diff(lon), diff(lat))
        size <- round(640 * size / max(size))
    }
    if (is.null(zoom)) {
        zoom <- min(MaxZoom(lat, lon), size = size)
    }
    destfile <- paste(fn, "-orig.png", sep = "")
    map <- GetMap(center = center, zoom = zoom, destfile = destfile,
                  maptype = maptype, size = size)
    size <- dim(map$myTile)
    if (format == "png") {
        png(paste(fn, "-res.png", sep=""), width=size[2],
            height=size[1])
    }
    else {
        pdf(paste(fn, "-res.pdf", sep = ""), width = size[2] / 72 * 3,
                  height = size[1] / 72 *  3)
    }
    PlotOnStaticMap(map, ne$N, ne$E, destfile = destfile,
                    size = c(size[2], size[1]), FUN = lines,
                    col = col, lwd = lwd)
    dev.off()
}

decodeGPS <- function(fn) {
    base.dir <- "/foo"
    file.name <- paste(base.dir, "/", fn, ".log", sep = "")
    a <- readLines(file.name)
    a <- a[-1]
    n <- length(a)/3
    quality <- hh <- mm <- ss <- integer(n)
    km <- m <- N <- E <- numeric(n)
    cnt <- 0
    for (i in 1:n * 3 - 2) {
        d <- unlist(strsplit(a[i], ","))
        stopifnot(d[1] == "$GPGGA")
        cnt <- cnt + 1
        hh[cnt] <- as.integer(substr(d[2], 1, 2))
        mm[cnt] <- as.integer(substr(d[2], 3, 4))
        ss[cnt] <- as.integer(substr(d[2], 5, 10))
        N[cnt] <- as.double(substr(d[3], 1, 2)) +
                  as.double(substr(d[3],  3, 9)) / 60
        E[cnt] <- as.double(substr(d[5], 1, 3)) +
                  as.double(substr(d[5],  4, 10)) / 60
        quality[cnt] <- as.integer(d[7])
        m[cnt] <- as.double(d[10])
        d <- unlist(strsplit(a[i + 1], ","))
        stopifnot(d[1] == "$GPRMC")
        d <- unlist(strsplit(a[i + 2], ","))
        stopifnot(d[1] == "$GPVTG")
        km[cnt] <- as.double(d[8])
    }
    m[quality == 0] <- NA
    km[quality == 0] <- NA
    d <- data.frame(hh, mm, ss, E, N, m, km, quality)
    layout(matrix(c(1, 1, 2, 3), 2))
    par(mar = c(3, 3, 0.5, 0.5), mgp = c(1.8, 0.6, 0))
    plot(E, N, type = "l", asp = 1)
    points(E[km == 0], N[km == 0], col = 2, pch = 19, cex = 0.2)
    plot(m, type = "l", ylab = "標高(m)")
    plot(km, type = "l", ylab = "対地速度(km)")
    layout(1)
    return(list(N=N, E=E))
}