#  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_latlong(584385) -> z
#  paste("latitide  = ", z$lat)
#  paste("longitude = ", z$lon)
#
# 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)
}