| Top Page | プログラミング | メッシュコード |
Updated on 2014-08-14
以下の内容のファイルへのリンク (mesh_lib.R)。 「リンク先を保存」等でダウンロードして使用。
# mesh_lib.R
#
# メッシュコード計算用関数群 (for R)
#
# 2010-07-13 TAKENAKA, Akio
# rev 2014-08-14 TAKENAKA, Akio (sapplyの用例の修正）
#
# source('mesh_lib.R') で、このファイルて定義した４つの関数を
# 読み込んでから使用する。
#
# meshcode_to_latlong, メッシュコードから緯度・経度へ
# latlong_to_meshcode, 緯度・経度からメッシュコードへ。
# dms_to_deg, 度・分・秒から度（小数点表示）へ。
# deg_to_dms, 度（小数点表示）から度・分・秒へ。
########################################
# メッシュコードを渡し、対応する緯度・経度を返す。
# 特に指定しなければメッシュの中央の座標を返す。
# メッシュコードは一次、ニ次、三次いずれも可。
#
# NW, NE, SW, SE を２つ目の引数として渡すと、それぞれ北西、北東、南西、南東の
# コーナーの座標を返す。
#
# 戻り値はリスト、要素の名前は lat と long、各要素の値は度単位の小数。
#
# Usage1:
# meshcode_to_long(584385) -> z
# paste("latitide = ", z$lat)
# paste("longitude = ", z$lat)
#
# Usage2:
# d <- read.table('...', header = T) # メッシュコードの列 (code) を持つテキストファイル
# lapply(d$code, meshcode_to_latlong) -> z # z は list の list
# d$lat <- sapply(z, "[[", "lat") # z の各要素に[["lat"]] を適用
# d$long <- sapply(z, "[[", "long") # long についても同様に。
meshcode_to_latlong <- function (code, loc = "C")
{
code <- as.character(code)
if (length(grep("^[0-9]{4}", code)) == 1) { # 一次メッシュ以上
code12 <- as.numeric(substring(code, 1, 2))
code34 <- as.numeric(substring(code, 3, 4))
lat_width <- 2 / 3;
long_width <- 1;
}
else {
return(NULL)
}
if (length(grep("^[0-9]{6}", code)) == 1) { # 二次メッシュ以上
code5 <- as.numeric(substring(code, 5, 5))
code6 <- as.numeric(substring(code, 6, 6))
lat_width <- lat_width / 8;
long_width <- long_width / 8;
}
if (length(grep("^[0-9]{8}", code)) == 1) { # 三次メッシュ
code7 <- as.numeric(substring(code, 7, 7))
code8 <- as.numeric(substring(code, 8, 8))
lat_width <- lat_width / 10;
long_width <- long_width / 10;
}
# 以下、南西コーナーの座標を求める。
lat <- code12 * 2 / 3; # 一次メッシュ
long <- code34 + 100;
if (exists("code5") && exists("code6")) { # 二次メッシュ
lat <- lat + (code5 * 2 / 3) / 8;
long <- long + code6 / 8;
}
if (exists("code7") && exists("code8")) { # 三次メッシュ
lat <- lat + (code7 * 2 / 3) / 8 / 10;
long <- long + code8 / 8 / 10;
}
if (loc == "C") { # 中央の座標
lat <- lat + lat_width / 2;
long <- long + long_width / 2;
}
if (length(grep("N", loc)) == 1) { # 北端の座標
lat <- lat + lat_width;
}
if (length(grep("E", loc) == 1)) { # 東端の座標
long <- long +long_width;
}
lat <- sprintf("%.8f", lat); # 小数点以下８桁まで。
long <- sprintf("%.8f", long);
x <- list(as.numeric(lat), as.numeric(long))
names(x) <- c("lat", "long")
return (x)
}
##################
# 緯度・経度（度単位）を受けとり、これを含むメッシュのコードを返す
#
# ３番めの引数として 1,2,3 を渡せば、一次メッシュ、二次メッシュ、
# 三次メッシュコードを計算する。指定がなければ三次メッシュ、
#
# Usage1:
# code_3 <- latlong_to_meshcode(35.3, 138.5) # 三次メッシュコード
# code_2 <- latlong_to_meshcode(35.3, 138.5, 2) # 二次メッシュコード
#
# Usage2:
# d <- read.table('lat_long.txt', header = T) # 緯度(lat)と経度(long)の列があるファイル
# mapply(latlong_to_meshcode, d$lat, d$long) -> z # z は 三次メッシュコードのベクトル
# d$code <- z # データフレームにメッシュコードの列を附加
#
# Usage2b:
# d <- read.table('lat_long.txt', header = T)
# mapply(latlong_to_meshcode, d$lat, d$long, MoreArgs = list(2)) -> z # 二次メッシュ
# d$code <- z
latlong_to_meshcode <- function(lat, long, order = 3)
{
if (length(grep("[123]", order)) == 0) {
return (NULL)
}
# Latitude
lat_in_min <- lat * 60
code12 <- as.integer(lat_in_min / 40)
lat_rest_in_min <- lat_in_min - code12 * 40
code5 = as.integer(lat_rest_in_min / 5 ); # 二次メッシュの１区画は緯度５分。
lat_rest_in_min <- lat_rest_in_min - code5 * 5
code7 <- as.integer(lat_rest_in_min / (5 / 10))
# Longitude
code34 <- as.integer(long) - 100
long_rest_in_deg <- long - as.integer(long)
code6 <- as.integer(long_rest_in_deg * 8)
long_rest_in_deg <- long_rest_in_deg - code6 / 8
code8 <- as.integer(long_rest_in_deg / (1/80) )
code <- sprintf("%02d%02d", code12, code34)
if (order >= 2) {
code <- sprintf("%s%01d%01d", code, code5, code6)
}
if (order >= 3) {
code <- sprintf("%s%01d%01d", code, code7, code8)
}
return (as.numeric(code))
}
##################
# 度、分、秒を受けとり、度の単位で小数点表示に変換して返す。
#
# Usage1:
# dms_to_deg(35, 20, 0) -> deg
#
# Usage2:
# d <- read.table('....txt', header = T) # 度(d), 分(m), 秒(s) の列があるファイル。
# mapply(dms_to_deg, d$d, d$m, d$s) -> z # １行めから順に 各列の要素を dms_to_deg へ渡す。
# # z は度の小数点表示が並んだベクトル。
# d$deg <- z # データフレームに新しい列を加える
dms_to_deg <- function(deg, min, sec)
{
in.deg <- deg + min / 60 + sec / 3600
in.deg.str <- sprintf("%.6f", in.deg) # 小数点以下6桁まで。
in.deg <- as.numeric(in.deg.str)
}
##################
# 度の単位の小数点表示を受けとり、 度、分、秒に変換して返す。
# 戻り値はリスト、要素の名前は deg, min, sec、各要素の値は数値
#
# Usage1:
# deg_to_dms(35.456) -> x
# paste(x$deg, x$min, x$sec, sep = ":") # 三つの要素を : をはさみながら連結。
# Usage2:
# d <- read.table('...', header = T) # 小数点表示の度(deg)の列があるファイル
# lapply(d$deg, deg_to_dms) -> z # z は list の list
# d$d <- sapply(z, "[[", "deg") # z の各要素に[["deg"]] を適用
# d$m <- sapply(z, "[[", "min") # "min" についても同様。
# d$s <- sapply(z, "[[", "sec") # "sec" についても同様。
deg_to_dms <- function(x)
{
deg <- as.integer(x)
rest_in_deg <- x - deg
min <- as.integer( rest_in_deg * 60 )
rest_in_min <- rest_in_deg * 60 - min
sec <- rest_in_min * 60
sec <- sprintf("%.2f", sec) # 小数点以下２桁まで。
sec <- as.numeric(sec)
# もし繰り上がってしまっていたら、その分の調整、
if (sec >= 60) { # 秒が60を越えたら
sec = sec - 60
min = min + 1
}
if (min == 60) { # 分が 60 を越えたら
min = min - 60
deg = deg + 1
}
x <- list(deg, min, sec) # 度、分、秒のベクトルを返す。
names(x) <- c("deg", "min", "sec")
return (x)
}