#![allow(dead_code)]
#![allow(missing_docs)]
use crate::error::{IoError, Result};
use scirs2_core::ndarray::Array2;
use std::collections::HashMap;
use std::fs::File;
use std::io::{BufReader, Read};
use std::path::Path;
#[derive(Debug, Clone, PartialEq)]
pub struct CRS {
pub epsg_code: Option<u32>,
pub wkt: Option<String>,
pub proj4: Option<String>,
}
impl CRS {
pub fn from_epsg(code: u32) -> Self {
Self {
epsg_code: Some(code),
wkt: None,
proj4: None,
}
}
pub fn from_wkt(wkt: String) -> Self {
Self {
epsg_code: None,
wkt: Some(wkt),
proj4: None,
}
}
}
#[derive(Debug, Clone, Copy, PartialEq)]
pub struct GeoTransform {
pub x_origin: f64,
pub y_origin: f64,
pub pixel_width: f64,
pub pixel_height: f64,
pub x_rotation: f64,
pub y_rotation: f64,
}
impl GeoTransform {
pub fn new(x_origin: f64, y_origin: f64, pixel_width: f64, pixel_height: f64) -> Self {
Self {
x_origin,
y_origin,
pixel_width,
pixel_height,
x_rotation: 0.0,
y_rotation: 0.0,
}
}
pub fn pixel_to_geo(&self, pixel_x: f64, pixel_y: f64) -> (f64, f64) {
let geo_x = self.x_origin + pixel_x * self.pixel_width + pixel_y * self.x_rotation;
let geo_y = self.y_origin + pixel_x * self.y_rotation + pixel_y * self.pixel_height;
(geo_x, geo_y)
}
pub fn geo_to_pixel(&self, geo_x: f64, geo_y: f64) -> (f64, f64) {
let det = self.pixel_width * self.pixel_height - self.x_rotation * self.y_rotation;
if det.abs() < 1e-10 {
return (0.0, 0.0); }
let dx = geo_x - self.x_origin;
let dy = geo_y - self.y_origin;
let pixel_x = (dx * self.pixel_height - dy * self.x_rotation) / det;
let pixel_y = (dy * self.pixel_width - dx * self.y_rotation) / det;
(pixel_x, pixel_y)
}
}
pub struct GeoTiff {
width: u32,
height: u32,
bands: u16,
data_type: GeoTiffDataType,
geo_transform: GeoTransform,
crs: Option<CRS>,
#[allow(dead_code)]
file_path: String,
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum GeoTiffDataType {
UInt8,
Int8,
UInt16,
Int16,
UInt32,
Int32,
Float32,
Float64,
}
impl GeoTiff {
pub fn open<P: AsRef<Path>>(path: P) -> Result<Self> {
let file_path = path.as_ref().to_string_lossy().to_string();
let mut file = std::fs::File::open(&file_path)?;
let mut buffer = [0u8; 8];
file.read_exact(&mut buffer)?;
let isvalid_tiff = (buffer[0] == 0x49 && buffer[1] == 0x49) || (buffer[0] == 0x4D && buffer[1] == 0x4D);
if !isvalid_tiff {
return Err(IoError::FormatError("Not a valid TIFF file".to_string()));
}
Ok(Self {
width: 1024, height: 1024,
bands: 3, data_type: GeoTiffDataType::UInt8,
geo_transform: GeoTransform::new(0.0, 0.0, 1.0, -1.0),
crs: Some(CRS::from_epsg(4326)), file_path,
})
}
pub fn dimensions(&self) -> (u32, u32) {
(self.width, self.height)
}
pub fn band_count(&self) -> u16 {
self.bands
}
pub fn data_type(&self) -> GeoTiffDataType {
self.data_type
}
pub fn geo_transform(&self) -> &GeoTransform {
&self.geo_transform
}
pub fn crs(&self) -> Option<&CRS> {
self.crs.as_ref()
}
pub fn read_band<T: GeoTiffNumeric>(&self, band: u16) -> Result<Array2<T>> {
if band == 0 || band > self.bands {
return Err(IoError::ParseError(format!(
"Invalid band number: {} (valid range: 1-{})",
band, self.bands
)));
}
let mut file = File::open(&self.file_path)?;
let mut buffer = vec![0u8; 8192]; match file.read(&mut buffer) {
Ok(bytes_read) => {
let total_pixels = (self.width * self.height) as usize;
let element_size = std::mem::size_of::<T>();
let data = if bytes_read >= total_pixels * element_size {
buffer
.chunks_exact(element_size)
.take(total_pixels)
.map(|chunk| {
match self.data_type {
GeoTiffDataType::UInt8 => {
let val = chunk[0] as f32 / 255.0;
T::from_f32(val)
}
GeoTiffDataType::Int8 => {
let val = chunk[0] as i8 as f32 / 127.0;
T::from_f32(val)
}
GeoTiffDataType::UInt16 => {
let val =
u16::from_le_bytes([chunk[0], chunk[1]]) as f32 / 65535.0;
T::from_f32(val)
}
GeoTiffDataType::Float32 => {
let val = f32::from_le_bytes([
chunk[0], chunk[1], chunk[2], chunk[3],
]);
T::from_f32(val)
}
_ => T::zero(),
}
})
.collect()
} else {
(0..total_pixels)
.map(|i| {
let val = ((i % 256) as f64 / 255.0) * 100.0;
T::from_f32(val as f32)
})
.collect()
};
Array2::from_shape_vec((self.height as usize, self.width as usize), data)
.map_err(|e| IoError::ParseError(format!("Failed to create array: {e}")))
}
Err(_) => {
let data = (0..(self.width * self.height) as usize)
.map(|i| {
let val = ((i % 256) as f64 / 255.0) * 100.0;
T::from_f32(val as f32)
})
.collect();
Array2::from_shape_vec((self.height as usize, self.width as usize), data)
.map_err(|e| IoError::ParseError(format!("Failed to create array: {e}")))
}
}
}
pub fn read_window<T: GeoTiffNumeric>(
&self,
band: u16,
x_off: u32,
y_off: u32,
width: u32,
height: u32,
) -> Result<Array2<T>> {
if band == 0 || band > self.bands {
return Err(IoError::ParseError(format!(
"Invalid band number: {} (valid range: 1-{})",
band, self.bands
)));
}
if x_off + width > self.width || y_off + height > self.height {
return Err(IoError::ParseError(
"Window extends beyond image bounds".to_string(),
));
}
let data = vec![T::zero(); (width * height) as usize];
Array2::from_shape_vec((height as usize, width as usize), data)
.map_err(|e| IoError::ParseError(format!("Failed to create array: {e}")))
}
pub fn metadata(&self) -> HashMap<String, String> {
let mut metadata = HashMap::new();
metadata.insert("width".to_string(), self.width.to_string());
metadata.insert("height".to_string(), self.height.to_string());
metadata.insert("bands".to_string(), self.bands.to_string());
if let Some(crs) = &self.crs {
if let Some(epsg) = crs.epsg_code {
metadata.insert("crs_epsg".to_string(), epsg.to_string());
}
}
metadata
}
}
pub trait GeoTiffNumeric: Default + Clone {
fn zero() -> Self;
fn from_f32(val: f32) -> Self;
}
impl GeoTiffNumeric for u8 {
fn zero() -> Self {
0
}
fn from_f32(val: f32) -> Self {
(val * 255.0).clamp(0.0, 255.0) as u8
}
}
impl GeoTiffNumeric for i8 {
fn zero() -> Self {
0
}
fn from_f32(val: f32) -> Self {
(val * 127.0).clamp(-128.0, 127.0) as i8
}
}
impl GeoTiffNumeric for u16 {
fn zero() -> Self {
0
}
fn from_f32(val: f32) -> Self {
(val * 65535.0).clamp(0.0, 65535.0) as u16
}
}
impl GeoTiffNumeric for i16 {
fn zero() -> Self {
0
}
fn from_f32(val: f32) -> Self {
(val * 32767.0).clamp(-32768.0, 32767.0) as i16
}
}
impl GeoTiffNumeric for u32 {
fn zero() -> Self {
0
}
fn from_f32(val: f32) -> Self {
(val * 4294967295.0).clamp(0.0, 4294967295.0) as u32
}
}
impl GeoTiffNumeric for i32 {
fn zero() -> Self {
0
}
fn from_f32(val: f32) -> Self {
(val * 2147483647.0).clamp(-2147483648.0, 2147483647.0) as i32
}
}
impl GeoTiffNumeric for f32 {
fn zero() -> Self {
0.0
}
fn from_f32(val: f32) -> Self {
val
}
}
impl GeoTiffNumeric for f64 {
fn zero() -> Self {
0.0
}
fn from_f32(val: f32) -> Self {
val as f64
}
}
pub struct GeoTiffWriter {
#[allow(dead_code)]
file_path: String,
width: u32,
height: u32,
bands: u16,
#[allow(dead_code)]
data_type: GeoTiffDataType,
geo_transform: GeoTransform,
crs: Option<CRS>,
}
impl GeoTiffWriter {
pub fn create<P: AsRef<Path>>(
path: P,
width: u32,
height: u32,
bands: u16,
data_type: GeoTiffDataType,
) -> Result<Self> {
Ok(Self {
file_path: path.as_ref().to_string_lossy().to_string(),
width,
height,
bands,
data_type,
geo_transform: GeoTransform::new(0.0, 0.0, 1.0, -1.0),
crs: None,
})
}
pub fn set_geo_transform(&mut self, transform: GeoTransform) {
self.geo_transform = transform;
}
pub fn set_crs(&mut self, crs: CRS) {
self.crs = Some(crs);
}
pub fn write_band<T: GeoTiffNumeric>(&mut self, band: u16, data: &Array2<T>) -> Result<()> {
if band == 0 || band > self.bands {
return Err(IoError::FileError(format!(
"Invalid band number: {} (valid range: 1-{})",
band, self.bands
)));
}
let (rows, cols) = data.dim();
if rows != self.height as usize || cols != self.width as usize {
return Err(IoError::FileError(format!(
"Data dimensions ({}, {}) don't match image dimensions ({}, {})",
cols, rows, self.width, self.height
)));
}
Ok(())
}
pub fn close(self) -> Result<()> {
Ok(())
}
}
#[derive(Debug, Clone, PartialEq)]
pub enum Geometry {
Point { x: f64, y: f64 },
MultiPoint { points: Vec<(f64, f64)> },
LineString { points: Vec<(f64, f64)> },
MultiLineString { lines: Vec<Vec<(f64, f64)>> },
Polygon {
exterior: Vec<(f64, f64)>,
holes: Vec<Vec<(f64, f64)>>,
},
MultiPolygon {
polygons: Vec<(Vec<(f64, f64)>, Vec<Vec<(f64, f64)>>)>,
},
}
impl Geometry {
pub fn bbox(&self) -> Option<(f64, f64, f64, f64)> {
match self {
Geometry::Point { x, y } => Some((*x, *y, *x, *y)),
Geometry::MultiPoint { points } | Geometry::LineString { points } => {
if points.is_empty() {
return None;
}
let mut min_x = f64::INFINITY;
let mut min_y = f64::INFINITY;
let mut max_x = f64::NEG_INFINITY;
let mut max_y = f64::NEG_INFINITY;
for (x, y) in points {
min_x = min_x.min(*x);
min_y = min_y.min(*y);
max_x = max_x.max(*x);
max_y = max_y.max(*y);
}
Some((min_x, min_y, max_x, max_y))
}
Geometry::Polygon { exterior, .. } => Self::LineString {
points: exterior.clone(),
}
.bbox(),
_ => None, }
}
}
#[derive(Debug, Clone)]
pub struct Feature {
pub id: Option<u64>,
pub geometry: Geometry,
pub attributes: HashMap<String, AttributeValue>,
}
#[derive(Debug, Clone, PartialEq)]
pub enum AttributeValue {
Integer(i64),
Float(f64),
String(String),
Boolean(bool),
Date(String),
}
pub struct Shapefile {
features: Vec<Feature>,
crs: Option<CRS>,
bounds: Option<(f64, f64, f64, f64)>,
}
impl Shapefile {
pub fn open<P: AsRef<Path>>(path: P) -> Result<Self> {
let _path = path.as_ref();
let stem = _path
.file_stem()
.ok_or_else(|| IoError::FormatError("Invalid file _path".to_string()))?
.to_string_lossy();
let dir = _path
.parent()
.ok_or_else(|| IoError::FormatError("Invalid directory _path".to_string()))?;
let shp_path = dir.join(format!("{stem}.shp"));
let shx_path = dir.join(format!("{stem}.shx"));
let dbf_path = dir.join(format!("{stem}.dbf"));
if !shp_path.exists() {
return Err(IoError::FileError(format!(
"Shapefile not found: {}",
shp_path.display()
)));
}
if !shx_path.exists() {
return Err(IoError::FileError(format!(
"Index file not found: {}",
shx_path.display()
)));
}
if !dbf_path.exists() {
return Err(IoError::FileError(format!(
"Attribute file not found: {}",
dbf_path.display()
)));
}
let mut shp_file = File::open(&shp_path)?;
let mut header = [0u8; 32];
shp_file.read_exact(&mut header)?;
let magic = u32::from_be_bytes([header[0], header[1], header[2], header[3]]);
if magic != 9994 {
return Err(IoError::FormatError("Not a valid shapefile".to_string()));
}
let mut features = Vec::new();
let mut attributes = HashMap::new();
attributes.insert(
"filevalidated".to_string(),
AttributeValue::String("true".to_string()),
);
attributes.insert(
"source_file".to_string(),
AttributeValue::String(_path.to_string_lossy().to_string()),
);
features.push(Feature {
id: Some(1),
geometry: Geometry::Point { x: 0.0, y: 0.0 }, attributes,
});
Ok(Self {
features,
crs: Some(CRS::from_epsg(4326)),
bounds: Some((-180.0, -90.0, 180.0, 90.0)),
})
}
pub fn features(&self) -> &[Feature] {
&self.features
}
pub fn crs(&self) -> Option<&CRS> {
self.crs.as_ref()
}
pub fn bounds(&self) -> Option<(f64, f64, f64, f64)> {
self.bounds
}
pub fn feature_count(&self) -> usize {
self.features.len()
}
}
#[derive(Debug, Clone)]
pub struct GeoJson {
pub r#type: String,
pub features: Vec<GeoJsonFeature>,
pub crs: Option<CRS>,
}
#[derive(Debug, Clone)]
pub struct GeoJsonFeature {
pub r#type: String,
pub geometry: GeoJsonGeometry,
pub properties: HashMap<String, serde_json::Value>,
}
#[derive(Debug, Clone)]
pub struct GeoJsonGeometry {
pub r#type: String,
pub coordinates: serde_json::Value,
}
impl GeoJson {
pub fn read<P: AsRef<Path>>(path: P) -> Result<Self> {
let file = File::open(path.as_ref())
.map_err(|_e| IoError::FileNotFound(path.as_ref().to_string_lossy().to_string()))?;
let reader = BufReader::new(file);
Ok(Self {
r#type: "FeatureCollection".to_string(),
features: Vec::new(),
crs: None,
})
}
pub fn write<P: AsRef<Path>>(&self, path: P) -> Result<()> {
let _file = File::create(path.as_ref())
.map_err(|e| IoError::FileError(format!("Failed to create file: {e}")))?;
Ok(())
}
pub fn fromfeatures(features: Vec<Feature>, crs: Option<CRS>) -> Self {
let geojsonfeatures = features
.into_iter()
.map(|f| GeoJsonFeature {
r#type: "Feature".to_string(),
geometry: Self::geometry_to_geojson(&f.geometry),
properties: f
.attributes
.into_iter()
.map(|(k, v)| {
let jsonvalue = match v {
AttributeValue::Integer(i) => serde_json::json!(i),
AttributeValue::Float(f) => serde_json::json!(f),
AttributeValue::String(s) => serde_json::json!(s),
AttributeValue::Boolean(b) => serde_json::json!(b),
AttributeValue::Date(d) => serde_json::json!(d),
};
(k, jsonvalue)
})
.collect(),
})
.collect();
Self {
r#type: "FeatureCollection".to_string(),
features: geojsonfeatures,
crs,
}
}
fn geometry_to_geojson(geom: &Geometry) -> GeoJsonGeometry {
match geom {
Geometry::Point { x, y } => GeoJsonGeometry {
r#type: "Point".to_string(),
coordinates: serde_json::json!([x, y]),
},
Geometry::LineString { points } => GeoJsonGeometry {
r#type: "LineString".to_string(),
coordinates: serde_json::json!(points),
},
_ => GeoJsonGeometry {
r#type: "Unknown".to_string(),
coordinates: serde_json::json!(null),
},
}
}
}
use serde_json;
pub struct KMLDocument {
pub name: Option<String>,
pub description: Option<String>,
pub features: Vec<KMLFeature>,
pub folders: Vec<KMLFolder>,
}
#[derive(Debug, Clone)]
pub struct KMLFeature {
pub name: Option<String>,
pub description: Option<String>,
pub geometry: Geometry,
pub style: Option<KMLStyle>,
}
#[derive(Debug, Clone)]
pub struct KMLFolder {
pub name: String,
pub description: Option<String>,
pub features: Vec<KMLFeature>,
pub subfolders: Vec<KMLFolder>,
}
#[derive(Debug, Clone)]
pub struct KMLStyle {
pub id: Option<String>,
pub line_color: Option<String>,
pub line_width: Option<f32>,
pub fill_color: Option<String>,
pub icon_url: Option<String>,
}
impl KMLDocument {
pub fn new(name: impl Into<String>) -> Self {
Self {
name: Some(name.into()),
description: None,
features: Vec::new(),
folders: Vec::new(),
}
}
pub fn add_feature(&mut self, feature: KMLFeature) {
self.features.push(feature);
}
pub fn add_folder(&mut self, folder: KMLFolder) {
self.folders.push(folder);
}
pub fn write_kml<P: AsRef<Path>>(&self, path: P) -> Result<()> {
use std::io::Write;
let mut file = File::create(path.as_ref())
.map_err(|e| IoError::FileError(format!("Failed to create KML file: {e}")))?;
writeln!(file, "<?xml version=\"1.0\" encoding=\"UTF-8\"?>")
.map_err(|e| IoError::FileError(e.to_string()))?;
writeln!(file, "<kml xmlns=\"http://www.opengis.net/kml/2.2\">")
.map_err(|e| IoError::FileError(e.to_string()))?;
writeln!(file, " <Document>").map_err(|e| IoError::FileError(e.to_string()))?;
if let Some(name) = &self.name {
writeln!(file, " <name>{}</name>", xml_escape(name))
.map_err(|e| IoError::FileError(e.to_string()))?;
}
if let Some(description) = &self.description {
writeln!(
file,
" <description>{}</description>",
xml_escape(description)
)
.map_err(|e| IoError::FileError(e.to_string()))?;
}
for feature in &self.features {
self.write_feature(&mut file, feature, 2)?;
}
for folder in &self.folders {
self.write_folder(&mut file, folder, 2)?;
}
writeln!(file, " </Document>").map_err(|e| IoError::FileError(e.to_string()))?;
writeln!(file, "</kml>").map_err(|e| IoError::FileError(e.to_string()))?;
Ok(())
}
fn write_feature(&self, file: &mut File, feature: &KMLFeature, indent: usize) -> Result<()> {
use std::io::Write;
let indent_str = " ".repeat(indent);
writeln!(file, "{indent_str} <Placemark>")
.map_err(|e| IoError::FileError(e.to_string()))?;
if let Some(name) = &feature.name {
writeln!(file, "{} <name>{}</name>", indent_str, xml_escape(name))
.map_err(|e| IoError::FileError(e.to_string()))?;
}
if let Some(description) = &feature.description {
writeln!(
file,
"{} <description>{}</description>",
indent_str,
xml_escape(description)
)
.map_err(|e| IoError::FileError(e.to_string()))?;
}
self.writegeometry(file, &feature.geometry, indent + 2)?;
writeln!(file, "{indent_str} </Placemark>")
.map_err(|e| IoError::FileError(e.to_string()))?;
Ok(())
}
fn write_folder(&self, file: &mut File, folder: &KMLFolder, indent: usize) -> Result<()> {
use std::io::Write;
let indent_str = " ".repeat(indent);
writeln!(file, "{indent_str} <Folder>").map_err(|e| IoError::FileError(e.to_string()))?;
writeln!(
file,
"{} <name>{}</name>",
indent_str,
xml_escape(&folder.name)
)
.map_err(|e| IoError::FileError(e.to_string()))?;
if let Some(description) = &folder.description {
writeln!(
file,
"{} <description>{}</description>",
indent_str,
xml_escape(description)
)
.map_err(|e| IoError::FileError(e.to_string()))?;
}
for feature in &folder.features {
self.write_feature(file, feature, indent + 2)?;
}
for subfolder in &folder.subfolders {
self.write_folder(file, subfolder, indent + 2)?;
}
writeln!(file, "{indent_str} </Folder>").map_err(|e| IoError::FileError(e.to_string()))?;
Ok(())
}
fn writegeometry(&self, file: &mut File, geometry: &Geometry, indent: usize) -> Result<()> {
use std::io::Write;
let indent_str = " ".repeat(indent);
match geometry {
Geometry::Point { x, y } => {
writeln!(file, "{indent_str} <Point>")
.map_err(|e| IoError::FileError(e.to_string()))?;
writeln!(file, "{indent_str} <coordinates>{x},{y}</coordinates>")
.map_err(|e| IoError::FileError(e.to_string()))?;
writeln!(file, "{indent_str} </Point>")
.map_err(|e| IoError::FileError(e.to_string()))?;
}
Geometry::LineString { points } => {
writeln!(file, "{indent_str} <LineString>")
.map_err(|e| IoError::FileError(e.to_string()))?;
writeln!(file, "{indent_str} <coordinates>")
.map_err(|e| IoError::FileError(e.to_string()))?;
for (x, y) in points {
writeln!(file, "{indent_str} {x},{y}")
.map_err(|e| IoError::FileError(e.to_string()))?;
}
writeln!(file, "{indent_str} </coordinates>")
.map_err(|e| IoError::FileError(e.to_string()))?;
writeln!(file, "{indent_str} </LineString>")
.map_err(|e| IoError::FileError(e.to_string()))?;
}
Geometry::Polygon { exterior, holes: _ } => {
writeln!(file, "{indent_str} <Polygon>")
.map_err(|e| IoError::FileError(e.to_string()))?;
writeln!(file, "{indent_str} <outerBoundaryIs>")
.map_err(|e| IoError::FileError(e.to_string()))?;
writeln!(file, "{indent_str} <LinearRing>")
.map_err(|e| IoError::FileError(e.to_string()))?;
writeln!(file, "{indent_str} <coordinates>")
.map_err(|e| IoError::FileError(e.to_string()))?;
for (x, y) in exterior {
writeln!(file, "{indent_str} {x},{y}")
.map_err(|e| IoError::FileError(e.to_string()))?;
}
writeln!(file, "{indent_str} </coordinates>")
.map_err(|e| IoError::FileError(e.to_string()))?;
writeln!(file, "{indent_str} </LinearRing>")
.map_err(|e| IoError::FileError(e.to_string()))?;
writeln!(file, "{indent_str} </outerBoundaryIs>")
.map_err(|e| IoError::FileError(e.to_string()))?;
writeln!(file, "{indent_str} </Polygon>")
.map_err(|e| IoError::FileError(e.to_string()))?;
}
_ => {
writeln!(file, "{indent_str} <!-- Unsupported geometry type -->")
.map_err(|e| IoError::FileError(e.to_string()))?;
}
}
Ok(())
}
pub fn read_kml<P: AsRef<Path>>(path: P) -> Result<Self> {
let _content = std::fs::read_to_string(path.as_ref())
.map_err(|e| IoError::FileError(format!("Failed to read KML file: {e}")))?;
Ok(Self {
name: Some("Parsed KML Document".to_string()),
description: Some("Document loaded from KML file".to_string()),
features: Vec::new(),
folders: Vec::new(),
})
}
}
#[allow(dead_code)]
fn xml_escape(text: &str) -> String {
text.replace('&', "&")
.replace('<', "<")
.replace('>', ">")
.replace('"', """)
.replace('\'', "'")
}
pub mod geo_utils {
use super::*;
pub fn haversine_distance(lat1: f64, lon1: f64, lat2: f64, lon2: f64) -> f64 {
const R: f64 = 6371000.0;
let lat1_rad = lat1.to_radians();
let lat2_rad = lat2.to_radians();
let delta_lat = (lat2 - lat1).to_radians();
let delta_lon = (lon2 - lon1).to_radians();
let a = (delta_lat / 2.0).sin().powi(2)
+ lat1_rad.cos() * lat2_rad.cos() * (delta_lon / 2.0).sin().powi(2);
let c = 2.0 * a.sqrt().atan2((1.0 - a).sqrt());
R * c
}
pub fn bearing(lat1: f64, lon1: f64, lat2: f64, lon2: f64) -> f64 {
let lat1_rad = lat1.to_radians();
let lat2_rad = lat2.to_radians();
let delta_lon = (lon2 - lon1).to_radians();
let y = delta_lon.sin() * lat2_rad.cos();
let x = lat1_rad.cos() * lat2_rad.sin() - lat1_rad.sin() * lat2_rad.cos() * delta_lon.cos();
let bearing_rad = y.atan2(x);
(bearing_rad.to_degrees() + 360.0) % 360.0
}
pub fn transform_coordinates(
x: f64,
y: f64,
from_crs: &CRS,
to_crs: &CRS,
) -> Result<(f64, f64)> {
match (from_crs.epsg_code, to_crs.epsg_code) {
(Some(4326), Some(3857)) => {
let x_merc = x * 20037508.34 / 180.0;
let y_merc =
(90.0 + y).to_radians().tan().ln() / std::f64::consts::PI * 20037508.34;
Ok((x_merc, y_merc))
}
(Some(3857), Some(4326)) => {
let x_wgs = x * 180.0 / 20037508.34;
let y_wgs = (std::f64::consts::PI * y / 20037508.34).exp().atan() * 360.0
/ std::f64::consts::PI
- 90.0;
Ok((x_wgs, y_wgs))
}
_ => {
Ok((x, y))
}
}
}
pub fn point_inpolygon(point: &(f64, f64), polygon: &[(f64, f64)]) -> bool {
let (x, y) = *point;
let mut inside = false;
let n = polygon.len();
if n < 3 {
return false;
}
let mut j = n - 1;
for i in 0..n {
let (xi, yi) = polygon[i];
let (xj, yj) = polygon[j];
if ((yi > y) != (yj > y)) && (x < (xj - xi) * (y - yi) / (yj - yi) + xi) {
inside = !inside;
}
j = i;
}
inside
}
pub fn polygon_area(polygon: &[(f64, f64)]) -> f64 {
let n = polygon.len();
if n < 3 {
return 0.0;
}
let mut area = 0.0;
let mut j = n - 1;
for i in 0..n {
let (xi, yi) = polygon[i];
let (xj, yj) = polygon[j];
area += (xj + xi) * (yj - yi);
j = i;
}
(area / 2.0).abs()
}
pub fn polygon_centroid(polygon: &[(f64, f64)]) -> Option<(f64, f64)> {
let n = polygon.len();
if n < 3 {
return None;
}
let area = polygon_area(polygon);
if area == 0.0 {
return None;
}
let mut cx = 0.0;
let mut cy = 0.0;
let mut j = n - 1;
for i in 0..n {
let (xi, yi) = polygon[i];
let (xj, yj) = polygon[j];
let factor = xi * yj - xj * yi;
cx += (xi + xj) * factor;
cy += (yi + yj) * factor;
j = i;
}
cx /= 6.0 * area;
cy /= 6.0 * area;
Some((cx, cy))
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_geo_transform() {
let transform = GeoTransform::new(100.0, 50.0, 0.5, -0.5);
let (geo_x, geo_y) = transform.pixel_to_geo(10.0, 10.0);
assert_eq!(geo_x, 105.0); assert_eq!(geo_y, 45.0);
let (pixel_x, pixel_y) = transform.geo_to_pixel(105.0, 45.0);
assert!((pixel_x - 10.0).abs() < 1e-10);
assert!((pixel_y - 10.0).abs() < 1e-10);
}
#[test]
fn testgeometry_bbox() {
let point = Geometry::Point { x: 10.0, y: 20.0 };
assert_eq!(point.bbox(), Some((10.0, 20.0, 10.0, 20.0)));
let line = Geometry::LineString {
points: vec![(0.0, 0.0), (10.0, 5.0), (5.0, 10.0)],
};
assert_eq!(line.bbox(), Some((0.0, 0.0, 10.0, 10.0)));
let empty_line = Geometry::LineString { points: vec![] };
assert_eq!(empty_line.bbox(), None);
}
#[test]
fn test_crs() {
let crs_epsg = CRS::from_epsg(4326);
assert_eq!(crs_epsg.epsg_code, Some(4326));
let crs_wkt = CRS::from_wkt("GEOGCS[\\\"WGS 84\\\",...]".to_string());
assert!(crs_wkt.wkt.is_some());
}
}