niszetの日記

アナログCMOS系雑用エンジニアが頑張る備忘録系日記

(R)震源レコードフォーマットを読み込む

固定文字長フォーマットの読み込み。

各日のページからデータを読み込むことが出来たので不要になったが、これが出来ないと思っていたので過去の震源データを下記から入手して確認しようとしていた。

www.data.jma.go.jp

欲しい年を選択すると、例えば2016年ならh2016.zipみたいなファイルがダウンロードできる。これを展開すれば2016なるファイルがある。拡張子はないが、これは文字長固定で区切られたファイルで、csvなどのカンマや空白で区切られたファイルではない。

このファイルフォーマットの説明はこちら。

www.data.jma.go.jp

これを読み混むことにする。これはreadr::read_lines()stringr::str_sub()を組み合わせれば出来る。

  ff <- readr::read_lines("h2016")

fx <- function(f){
  d <- data.frame(
  record = stringr::str_sub(f, 1, 1), # レコード種別ヘッダ
  year = stringr::str_sub(f, 2, 5),
  month = stringr::str_sub(f, 6, 7),
  day = stringr::str_sub(f, 8, 9),
  hour = stringr::str_sub(f, 10, 11),
  min = stringr::str_sub(f, 12, 13),
  sec = stringr::str_sub(f, 14, 17),
  err_sec = stringr::str_sub(f, 18, 21), # 標準誤差(秒)
  lat_deg = stringr::str_sub(f, 22, 24), # 緯度(度)
  lat_min = stringr::str_sub(f, 25, 28), # 緯度(分)
  err_lat_min = stringr::str_sub(f, 29, 32), # 標準誤差(分)
  long_deg = stringr::str_sub(f, 33, 36), # 経度(度)
  long_min = stringr::str_sub(f, 37, 40), # 経度(分)
  err_long_min = stringr::str_sub(f, 41, 44), # 標準誤差(分)
  depth = stringr::str_sub(f, 45, 49),
  err_depth = stringr::str_sub(f, 50, 52),
  magnitude1 =stringr::str_sub(f, 53, 54),
  magnitude1_type = stringr::str_sub(f, 55, 55),
  magnitude2 = stringr::str_sub(f, 56, 57),
  magnitude2_type = stringr::str_sub(f, 58, 58),
  used_soujihyo = stringr::str_sub(f, 59, 59),
  shingen_hyoka = stringr::str_sub(f, 60, 60),
  shingen_hojo = stringr::str_sub(f, 61, 61),
  max_shindo = stringr::str_sub(f, 62, 62),
  higai = stringr::str_sub(f, 63, 63),
  tsunami = stringr::str_sub(f, 64, 64),
  daichiku = stringr::str_sub(f, 65, 65),
  shochiku = stringr::str_sub(f, 66, 68),
  shino_chimei = stringr::str_sub(f, 69, 92),
  kansoku_tensu = stringr::str_sub(f, 93, 95),
  shingen_kettei_flag = stringr::str_sub(f, 96, 96)
  )
  d
}

fx(ff) -> df

これで良い。ただし得られるものは全て文字なので、適当に型変換する。depthは" 0 4" みたいなデータがあり、これを適切に変換する方法がわからないのでNAが返るようにしている(暗黙にNAになるので他のデータのエラーを見逃す可能性があるため注意。エラー処理を入れることをお勧めします。マグニチュードは例にCまでしか書いてなかったのでこのようにした。D以降があれば定義追加すればよいでしょう。 " " をNAとするかそのままにするかは好みかもしれません。NAでない場合の方が取り扱いやすいこともあるので。 たまに表示幅あわせで空白が入ってるケースがあるので注意。明示的にas.integerせずにwrite_csvしてread_csvしても、"00"は文字列の扱いなので注意、など細かい点に注意が必要です(下記では明示的に変換をしている)。 また、記号ではよくわからないもの(震度など)は別途変換関数を作ってあげると良さそうです(やってません)

fx2 <- function(df){
  df$year <- as.integer(df$year)

  df$month <- as.integer(df$month)
  df$day <- as.integer(df$day)
  df$hour <- as.integer(df$hour)
  df$min <- as.integer(df$min)
  df$sec <- as.integer(df$sec)/100
  df$err_sec <- as.integer(df$err_sec)/100
  df$lat_deg <- as.integer(str_remove(df$lat_deg, " ")) # remove space between "-" and number.
  df$lat_min <- as.integer(df$lat_min)/100
  df$err_lat_min <- as.integer(df$err_lat_min)/100
  df$long_deg <- as.integer(str_remove(df$long_deg, " ")) # remove space between "-" and number.
  df$long_min <- as.integer(df$long_min)/100
  df$err_long_min <- as.integer(df$err_long_min)/100

  df$depth = as.integer(depth)/100 # "注意。  0 4" みたいな謎フォーマットはNAになる。
  df$err_depth <- as.integer(df$err_depth)/100
  df$magnitude1 <- sup_magnitude(df$magnitude1)
  df$magnitude1_type = dplyr::na_if(df$magnitude1_type , " ")
  df$magnitude2 <- sup_magnitude(df$magnitude2)
  df$magnitude2_type = dplyr::na_if(df$magnitude2_type , " ")
  df$used_soujihyo = dplyr::na_if(df$used_soujihyo , " ")
  df$shingen_hyoka = dplyr::na_if(df$shingen_hyoka , " ")
  df$shingen_hojo = dplyr::na_if(df$shingen_hojo , " ")
  df$max_shindo = dplyr::na_if(df$max_shindo , " ")
  df$higai = dplyr::na_if(df$higai , " ")
  df$tsunami = dplyr::na_if(df$tsunami , " ")

  df$daichiku <- as.integer(df$daichiku)
  df$shochiku <- as.integer(df$shochiku)
  df$kansoku_tensu <- as.integer(df$kansoku_tensu)

  df
}

# 補助関数
sup_magnitude <- function(mag){
  mag <- stringr::str_replace(mag, "A", "-1")
  mag <- stringr::str_replace(mag, "B", "-2")
  mag <- stringr::str_replace(mag, "C", "-3")
  mag <- as.integer(mag)/10
  mag
}

fx2(df) -> df