use encoding_rs::SHIFT_JIS;
use jgdtrans::{Format, ParData, Parameter, Point, Transformer};
use std::collections::HashMap;
use std::fs;
use std::path::Path;
use thiserror::Error;
#[derive(Error, Debug)]
pub enum TransformError {
#[error("パラメータファイル読み込みエラー: {0}")]
FileRead(#[from] std::io::Error),
#[error("パラメータパースエラー: {0}")]
Parse(String),
#[error("変換エラー: {0}")]
Transform(String),
#[error("パラメータディレクトリが存在しません: {0}")]
DirectoryNotFound(String),
}
#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord)]
pub struct SurveyStartDate {
pub year: i32,
pub month: u32,
pub day: u32,
}
impl SurveyStartDate {
pub fn new(year: i32, month: u32, day: u32) -> Self {
Self { year, month, day }
}
pub fn from_yyyymmdd(s: &str) -> Option<Self> {
if s.len() < 8 {
return None;
}
let year = s[0..4].parse().ok()?;
let month = s[4..6].parse().ok()?;
let day = s[6..8].parse().ok()?;
Some(Self { year, month, day })
}
pub fn from_parts(year: Option<&str>, month: Option<&str>, day: Option<&str>) -> Option<Self> {
Some(Self {
year: year?.parse().ok()?,
month: month?.parse().ok()?,
day: day?.parse().ok()?,
})
}
}
struct EarthquakeInfo {
file_pattern: &'static str,
date: SurveyStartDate,
}
const EARTHQUAKES: &[EarthquakeInfo] = &[
EarthquakeInfo {
file_pattern: "tokachi2003",
date: SurveyStartDate { year: 2003, month: 9, day: 26 },
},
EarthquakeInfo {
file_pattern: "fukuoka2005",
date: SurveyStartDate { year: 2005, month: 3, day: 20 },
},
EarthquakeInfo {
file_pattern: "noto2007",
date: SurveyStartDate { year: 2007, month: 3, day: 25 },
},
EarthquakeInfo {
file_pattern: "chuetsuoki2007",
date: SurveyStartDate { year: 2007, month: 7, day: 16 },
},
EarthquakeInfo {
file_pattern: "iwatemiyagi2008",
date: SurveyStartDate { year: 2008, month: 6, day: 14 },
},
EarthquakeInfo {
file_pattern: "touhokutaiheiyouoki2011",
date: SurveyStartDate { year: 2011, month: 3, day: 11 },
},
EarthquakeInfo {
file_pattern: "kumamoto2016",
date: SurveyStartDate { year: 2016, month: 4, day: 14 },
},
];
fn should_apply_earthquake_correction(file_name: &str, survey_date: Option<SurveyStartDate>) -> bool {
let earthquake = EARTHQUAKES.iter().find(|eq| file_name.contains(eq.file_pattern));
match (earthquake, survey_date) {
(Some(eq), Some(date)) => {
date < eq.date
}
(Some(_), None) => false,
(None, _) => true,
}
}
pub fn read_param_file<P: AsRef<Path>>(path: P) -> Result<String, std::io::Error> {
let bytes = fs::read(&path)?;
if let Ok(s) = String::from_utf8(bytes.clone()) {
return Ok(s);
}
let (decoded, _, _) = SHIFT_JIS.decode(&bytes);
Ok(decoded.into_owned())
}
pub struct JgdTransformer {
tky_tf: Transformer<ParData>,
unified_tf: Transformer<ParData>,
jgd2011_to_wgs84: proj::Proj,
}
impl JgdTransformer {
pub fn with_survey_date<P: AsRef<Path>>(
tky_path: P,
patch_dir: P,
survey_date: SurveyStartDate,
) -> Result<Self, TransformError> {
let tky_content = read_param_file(&tky_path)?;
let tky_tf = Transformer::from_str(&tky_content, Format::TKY2JGD)
.map_err(|e| TransformError::Parse(format!("TKY2JGD parse error: {}", e)))?;
let unified_tf = Self::create_unified_transformer(patch_dir, Some(survey_date))?;
let jgd2011_to_wgs84 = proj::Proj::new_known_crs("EPSG:6668", "EPSG:4326", None)
.map_err(|e| TransformError::Transform(format!("proj init error: {}", e)))?;
Ok(Self { tky_tf, unified_tf, jgd2011_to_wgs84 })
}
pub fn new<P: AsRef<Path>>(tky_path: P, patch_dir: P) -> Result<Self, TransformError> {
let tky_content = read_param_file(&tky_path)?;
let tky_tf = Transformer::from_str(&tky_content, Format::TKY2JGD)
.map_err(|e| TransformError::Parse(format!("TKY2JGD parse error: {}", e)))?;
let unified_tf = Self::create_unified_transformer(patch_dir, None)?;
let jgd2011_to_wgs84 = proj::Proj::new_known_crs("EPSG:6668", "EPSG:4326", None)
.map_err(|e| TransformError::Transform(format!("proj init error: {}", e)))?;
Ok(Self { tky_tf, unified_tf, jgd2011_to_wgs84 })
}
pub fn new_tky_only<P: AsRef<Path>>(tky_path: P) -> Result<Self, TransformError> {
let tky_content = read_param_file(&tky_path)?;
let tky_tf = Transformer::from_str(&tky_content, Format::TKY2JGD)
.map_err(|e| TransformError::Parse(format!("TKY2JGD parse error: {}", e)))?;
let empty_params: HashMap<u32, Parameter> = HashMap::new();
let empty_data =
ParData::with_description(Format::PatchJGD_HV, empty_params, "Empty".to_string());
let unified_tf = Transformer::new(empty_data);
let jgd2011_to_wgs84 = proj::Proj::new_known_crs("EPSG:6668", "EPSG:4326", None)
.map_err(|e| TransformError::Transform(format!("proj init error: {}", e)))?;
Ok(Self { tky_tf, unified_tf, jgd2011_to_wgs84 })
}
fn create_unified_transformer<P: AsRef<Path>>(
dir_path: P,
survey_date: Option<SurveyStartDate>,
) -> Result<Transformer<ParData>, TransformError> {
let dir = dir_path.as_ref();
if !dir.exists() {
let empty_params: HashMap<u32, Parameter> = HashMap::new();
let empty_data = ParData::with_description(
Format::PatchJGD_HV,
empty_params,
"No patch files".to_string(),
);
return Ok(Transformer::new(empty_data));
}
let mut unified_params: HashMap<u32, Parameter> = HashMap::new();
let mut description = String::from("Unified Parameter Set:\n");
if let Some(date) = survey_date {
description.push_str(&format!(
"Survey date: {}/{}/{}\n",
date.year, date.month, date.day
));
} else {
description.push_str("Survey date: unknown (earthquake corrections disabled)\n");
}
let mut loaded_count = 0;
let horizontal_dir = dir.join("horizontal");
let elevation_dir = dir.join("elevation");
if horizontal_dir.exists() {
description.push_str("\n[Horizontal corrections]\n");
loaded_count += Self::load_param_files(
&horizontal_dir,
Format::PatchJGD,
survey_date,
&mut unified_params,
&mut description,
)?;
}
if elevation_dir.exists() {
description.push_str("\n[Elevation corrections]\n");
loaded_count += Self::load_param_files_elevation(
&elevation_dir,
survey_date,
&mut unified_params,
&mut description,
)?;
}
if loaded_count == 0 {
description = "No patch files loaded".to_string();
}
let unified_data =
ParData::with_description(Format::PatchJGD_HV, unified_params, description);
Ok(Transformer::new(unified_data))
}
fn load_param_files(
dir: &Path,
format: Format,
survey_date: Option<SurveyStartDate>,
unified_params: &mut HashMap<u32, Parameter>,
description: &mut String,
) -> Result<usize, TransformError> {
let mut loaded_count = 0;
for entry in fs::read_dir(dir)? {
let entry = entry?;
let path = entry.path();
if path.extension().and_then(|s| s.to_str()) != Some("par") {
continue;
}
let file_name = match path.file_name().and_then(|s| s.to_str()) {
Some(name) => name.to_string(),
None => continue,
};
if !should_apply_earthquake_correction(&file_name, survey_date) {
description.push_str(&format!(" - {} (skipped: survey after earthquake)\n", file_name));
continue;
}
let content = read_param_file(&path)?;
let par_data = ParData::from_str(&content, format)
.map_err(|e| TransformError::Parse(format!("{}: {}", file_name, e)))?;
for (mesh_code, param) in par_data.parameter {
unified_params
.entry(mesh_code)
.and_modify(|existing| {
existing.latitude += param.latitude;
existing.longitude += param.longitude;
existing.altitude += param.altitude;
})
.or_insert(param);
}
description.push_str(&format!(" - {} (applied)\n", file_name));
loaded_count += 1;
}
Ok(loaded_count)
}
fn load_param_files_elevation(
dir: &Path,
survey_date: Option<SurveyStartDate>,
unified_params: &mut HashMap<u32, Parameter>,
description: &mut String,
) -> Result<usize, TransformError> {
let mut loaded_count = 0;
for entry in fs::read_dir(dir)? {
let entry = entry?;
let path = entry.path();
if path.extension().and_then(|s| s.to_str()) != Some("par") {
continue;
}
let file_name = match path.file_name().and_then(|s| s.to_str()) {
Some(name) => name.to_string(),
None => continue,
};
if !should_apply_earthquake_correction(&file_name, survey_date) {
description.push_str(&format!(" - {} (skipped: survey after earthquake)\n", file_name));
continue;
}
let format = if file_name.contains("hyoko") {
Format::HyokoRev } else {
Format::PatchJGD_H };
let content = read_param_file(&path)?;
let par_data = ParData::from_str(&content, format)
.map_err(|e| TransformError::Parse(format!("{}: {}", file_name, e)))?;
for (mesh_code, param) in par_data.parameter {
unified_params
.entry(mesh_code)
.and_modify(|existing| {
existing.latitude += param.latitude;
existing.longitude += param.longitude;
existing.altitude += param.altitude;
})
.or_insert(param);
}
description.push_str(&format!(" - {} (applied)\n", file_name));
loaded_count += 1;
}
Ok(loaded_count)
}
pub fn transform_horizontal(
&self,
lat: f64,
lon: f64,
system_code: &str,
) -> Option<(f64, f64)> {
let result = self.transform_point(lat, lon, 0.0, system_code).ok()?;
Some((result.latitude, result.longitude))
}
pub fn transform_full(
&self,
lat: f64,
lon: f64,
alt: f64,
system_code: &str,
) -> Result<(f64, f64, f64), TransformError> {
let result = self.transform_point(lat, lon, alt, system_code)?;
Ok((result.latitude, result.longitude, result.altitude))
}
pub fn to_wgs84(
&self,
lat: f64,
lon: f64,
system_code: &str,
) -> Option<(f64, f64)> {
let (lat_jgd2011, lng_jgd2011) = self.transform_horizontal(lat, lon, system_code)?;
self.jgd2011_to_wgs84.convert((lng_jgd2011, lat_jgd2011)).ok()
}
fn transform_point(
&self,
lat: f64,
lon: f64,
alt: f64,
system_code: &str,
) -> Result<Point, TransformError> {
let input_point = Point::new_unchecked(lat, lon, alt);
let code = system_code.trim().trim_start_matches('0');
match code {
"" | "0" => {
let corr_tky = self
.tky_tf
.forward_corr(&input_point)
.map_err(|e| TransformError::Transform(format!("TKY2JGD: {}", e)))?;
let p_2000 = &input_point + corr_tky;
match self.unified_tf.forward_corr(&p_2000) {
Ok(corr_unified) => Ok(&p_2000 + corr_unified),
Err(_) => Ok(p_2000), }
}
"1" => match self.unified_tf.forward_corr(&input_point) {
Ok(corr_unified) => Ok(&input_point + corr_unified),
Err(_) => Ok(input_point), },
"2" => Ok(input_point),
_ => Ok(input_point),
}
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_code_normalization() {
let code = "00".trim().trim_start_matches('0');
assert_eq!(code, "");
let code = "01".trim().trim_start_matches('0');
assert_eq!(code, "1");
let code = "02".trim().trim_start_matches('0');
assert_eq!(code, "2");
}
#[test]
fn test_tky2jgd_transform() {
let workspace_root = std::path::Path::new(env!("CARGO_MANIFEST_DIR")).parent().unwrap();
let tky_path = workspace_root.join("params/TKY2JGD.par");
if !tky_path.exists() {
println!("Skipping test: TKY2JGD.par not found at {:?}", tky_path);
return;
}
let transformer = JgdTransformer::new_tky_only(tky_path)
.expect("Failed to create transformer");
let (lat, lon) = transformer
.transform_horizontal(35.666667, 139.766667, "00")
.expect("Transform failed");
assert!((lat - 35.666667).abs() < 0.01, "Latitude shift too large");
assert!((lon - 139.766667).abs() < 0.01, "Longitude shift too large");
assert!((lat - 35.666667).abs() > 0.0001, "Latitude should have changed");
assert!((lon - 139.766667).abs() > 0.0001, "Longitude should have changed");
println!("Tokyo Datum: (35.666667, 139.766667)");
println!("JGD2000: ({}, {})", lat, lon);
}
#[test]
fn test_jgd2011_no_transform() {
let workspace_root = std::path::Path::new(env!("CARGO_MANIFEST_DIR")).parent().unwrap();
let tky_path = workspace_root.join("params/TKY2JGD.par");
if !tky_path.exists() {
return;
}
let transformer = JgdTransformer::new_tky_only(&tky_path)
.expect("Failed to create transformer");
let (lat, lon) = transformer
.transform_horizontal(35.666667, 139.766667, "02")
.expect("Transform failed");
assert!((lat - 35.666667).abs() < 0.000001);
assert!((lon - 139.766667).abs() < 0.000001);
}
#[test]
fn test_full_transform_with_survey_date_before_earthquake() {
let workspace_root = std::path::Path::new(env!("CARGO_MANIFEST_DIR")).parent().unwrap();
let tky_path = workspace_root.join("params/TKY2JGD.par");
let patch_dir = workspace_root.join("params/patch_files");
if !tky_path.exists() || !patch_dir.exists() {
println!("Skipping test: params not found");
return;
}
let survey_date = SurveyStartDate::new(2010, 5, 15);
let transformer = JgdTransformer::with_survey_date(&tky_path, &patch_dir, survey_date)
.expect("Failed to create transformer");
let result = transformer
.transform_full(38.2601, 140.8824, 15.0, "01")
.expect("Transform failed");
println!("Survey date: 2010/5/15 (before Tohoku earthquake)");
println!("JGD2000: (38.2601, 140.8824, 15.0m)");
println!("JGD2011: ({}, {}, {}m)", result.0, result.1, result.2);
}
#[test]
fn test_full_transform_with_survey_date_after_earthquake() {
let workspace_root = std::path::Path::new(env!("CARGO_MANIFEST_DIR")).parent().unwrap();
let tky_path = workspace_root.join("params/TKY2JGD.par");
let patch_dir = workspace_root.join("params/patch_files");
if !tky_path.exists() || !patch_dir.exists() {
println!("Skipping test: params not found");
return;
}
let survey_date = SurveyStartDate::new(2015, 6, 20);
let transformer = JgdTransformer::with_survey_date(&tky_path, &patch_dir, survey_date)
.expect("Failed to create transformer");
let result = transformer
.transform_full(38.2601, 140.8824, 15.0, "01")
.expect("Transform failed");
println!("Survey date: 2015/6/20 (after Tohoku earthquake)");
println!("JGD2000: (38.2601, 140.8824, 15.0m)");
println!("JGD2011: ({}, {}, {}m)", result.0, result.1, result.2);
}
#[test]
fn test_survey_date_parsing() {
let date = SurveyStartDate::from_yyyymmdd("20101015");
assert!(date.is_some());
let date = date.unwrap();
assert_eq!(date.year, 2010);
assert_eq!(date.month, 10);
assert_eq!(date.day, 15);
let date = SurveyStartDate::from_parts(Some("2011"), Some("3"), Some("10"));
assert!(date.is_some());
let date = date.unwrap();
assert_eq!(date.year, 2011);
assert_eq!(date.month, 3);
assert_eq!(date.day, 10);
}
#[test]
fn test_earthquake_correction_should_apply() {
let before = SurveyStartDate::new(2010, 1, 1);
assert!(should_apply_earthquake_correction("touhokutaiheiyouoki2011.par", Some(before)));
let after = SurveyStartDate::new(2012, 1, 1);
assert!(!should_apply_earthquake_correction("touhokutaiheiyouoki2011.par", Some(after)));
assert!(should_apply_earthquake_correction("hyokorev2014_geoid2011_h.par", Some(before)));
assert!(should_apply_earthquake_correction("hyokorev2014_geoid2011_h.par", Some(after)));
assert!(should_apply_earthquake_correction("hyokorev2014_geoid2011_h.par", None));
assert!(!should_apply_earthquake_correction("touhokutaiheiyouoki2011.par", None));
}
}