Merge pull request #21 from georust/feature/coordinate-transformation

Coordinate Transformations

Cargo.toml

@@ -3,10 +3,22 @@ name = "geotiff"
 description = "A GeoTIFF library for Rust"
 version = "0.0.2"
 edition = "2021"
-authors = ["Dominik Bucher <[email protected]>"]
+authors = [
+    "Dominik Bucher <[email protected]>",
+    "Gunnar Schulze <[email protected]>"
+]
 repository = "https://github.com/georust/geotiff"
 
 [dependencies]
+delaunator = { version = "1.0", optional = true }
+geo-index = { version = "0.1", optional = true }
+geo-types = { version = "0.7" }
 num_enum = "0.7"
 num-traits = "0.2"
 tiff = "0.9"
+
+[dev-dependencies]
+proj = "0.27"
+
+[features]
+tie-points = ["dep:delaunator", "dep:geo-index"]

src/coordinate_transform.rs

@@ -0,0 +1,179 @@
+#[cfg(feature = "tie-points")]
+use std::rc::Rc;
+
+#[cfg(feature = "tie-points")]
+use geo_index::rtree::OwnedRTree;
+use geo_types::Coord;
+use tiff::{TiffError, TiffFormatError, TiffResult};
+
+#[cfg(feature = "tie-points")]
+use crate::coordinate_transform::tie_points::Face;
+
+mod affine_transform;
+mod tie_point_and_pixel_scale;
+#[cfg(feature = "tie-points")]
+mod tie_points;
+
+const MODEL_TIE_POINT_TAG: &str = "ModelTiePointTag";
+const MODEL_PIXEL_SCALE_TAG: &str = "ModelPixelScaleTag";
+const MODEL_TRANSFORMATION_TAG: &str = "ModelTransformationTag";
+
+/// Defines the transformation between raster space and model space.
+///
+/// Ref: https://docs.ogc.org/is/19-008r4/19-008r4.html#_raster_to_model_coordinate_transformation_requirements
+#[derive(Debug)]
+pub enum CoordinateTransform {
+    AffineTransform {
+        transform: [f64; 6],
+        inverse_transform: [f64; 6],
+    },
+    TiePointAndPixelScale {
+        raster_point: Coord,
+        model_point: Coord,
+        pixel_scale: Coord,
+    },
+    #[cfg(feature = "tie-points")]
+    TiePoints {
+        raster_mesh: Rc<Vec<Face>>,
+        raster_index: OwnedRTree<f64>,
+        model_mesh: Rc<Vec<Face>>,
+        model_index: OwnedRTree<f64>,
+    },
+}
+
+impl CoordinateTransform {
+    pub(super) fn from_tag_data(
+        pixel_scale_data: Option<Vec<f64>>,
+        model_tie_points_data: Option<Vec<f64>>,
+        model_transformation_data: Option<Vec<f64>>,
+    ) -> TiffResult<Self> {
+        let pixel_scale = pixel_scale_data
+            .map(|data| {
+                <[f64; 3]>::try_from(data).map_err(|_| {
+                    TiffError::FormatError(TiffFormatError::Format(format!(
+                        "Number values in {MODEL_PIXEL_SCALE_TAG} must be equal to 3"
+                    )))
+                })
+            })
+            .transpose()?;
+        let tie_points = model_tie_points_data
+            .map(|data| {
+                let len = data.len();
+                if len == 0 {
+                    return Err(TiffError::FormatError(TiffFormatError::Format(format!(
+                        "Number of values in {MODEL_TIE_POINT_TAG} must be greater than 0"
+                    ))));
+                }
+
+                if len % 6 != 0 {
+                    return Err(TiffError::FormatError(TiffFormatError::Format(format!(
+                        "Number of values in {MODEL_TIE_POINT_TAG} must be divisible by 6"
+                    ))));
+                }
+
+                Ok(data)
+            })
+            .transpose()?;
+        let transformation_matrix = model_transformation_data
+            .map(|data| {
+                <[f64; 16]>::try_from(data).map_err(|_| {
+                    TiffError::FormatError(TiffFormatError::Format(format!(
+                        "Number of values in {MODEL_TRANSFORMATION_TAG} must be equal to 16"
+                    )))
+                })
+            })
+            .transpose()?;
+
+        if let Some(transformation_matrix) = transformation_matrix {
+            if pixel_scale.is_some() {
+                return Err(TiffError::FormatError(TiffFormatError::Format(
+                    format!("{MODEL_PIXEL_SCALE_TAG} must not be specified when {MODEL_TRANSFORMATION_TAG} is present"),
+                )));
+            }
+            if tie_points.is_some() {
+                return Err(TiffError::FormatError(TiffFormatError::Format(
+                    format!("{MODEL_TIE_POINT_TAG} must not be specified when {MODEL_TRANSFORMATION_TAG} is present"),
+                )));
+            }
+
+            Self::from_transformation_matrix(transformation_matrix)
+        } else {
+            let Some(tie_points) = tie_points else {
+                return Err(TiffError::FormatError(TiffFormatError::Format(
+                    format!("{MODEL_TIE_POINT_TAG} must be present when {MODEL_TRANSFORMATION_TAG} is missing"),
+                )));
+            };
+
+            if tie_points.len() == 6 {
+                let Some(pixel_scale) = pixel_scale else {
+                    return Err(TiffError::FormatError(TiffFormatError::Format(
+                        format!("{MODEL_PIXEL_SCALE_TAG} must be specified when {MODEL_TIE_POINT_TAG} contains 6 values"),
+                    )));
+                };
+
+                Self::from_tie_point_and_pixel_scale(&tie_points, &pixel_scale)
+            } else {
+                #[cfg(feature = "tie-points")]
+                {
+                    Self::from_tie_points(&tie_points)
+                }
+                #[cfg(not(feature = "tie-points"))]
+                {
+                    Err(TiffError::FormatError(TiffFormatError::Format(
+                        "Transformation by tie points is not supported".into(),
+                    )))
+                }
+            }
+        }
+    }
+
+    pub fn transform_to_model(&self, coord: &Coord) -> Coord {
+        match self {
+            CoordinateTransform::AffineTransform { transform, .. } => {
+                Self::transform_by_affine_transform(transform, coord)
+            }
+            CoordinateTransform::TiePointAndPixelScale {
+                raster_point,
+                model_point,
+                pixel_scale,
+            } => Self::transform_to_model_by_tie_point_and_pixel_scale(
+                raster_point,
+                model_point,
+                pixel_scale,
+                coord,
+            ),
+            #[cfg(feature = "tie-points")]
+            CoordinateTransform::TiePoints {
+                raster_index,
+                raster_mesh,
+                model_mesh,
+                ..
+            } => Self::transform_by_tie_points(raster_index, raster_mesh, model_mesh, coord),
+        }
+    }
+
+    pub(super) fn transform_to_raster(&self, coord: &Coord) -> Coord {
+        match self {
+            CoordinateTransform::AffineTransform {
+                inverse_transform, ..
+            } => Self::transform_by_affine_transform(inverse_transform, coord),
+            CoordinateTransform::TiePointAndPixelScale {
+                raster_point,
+                model_point,
+                pixel_scale,
+            } => Self::transform_to_raster_by_tie_point_and_pixel_scale(
+                raster_point,
+                model_point,
+                pixel_scale,
+                coord,
+            ),
+            #[cfg(feature = "tie-points")]
+            CoordinateTransform::TiePoints {
+                model_index,
+                model_mesh,
+                raster_mesh,
+                ..
+            } => Self::transform_by_tie_points(model_index, model_mesh, raster_mesh, coord),
+        }
+    }
+}

src/coordinate_transform/affine_transform.rs

@@ -0,0 +1,45 @@
+use geo_types::Coord;
+use tiff::{TiffError, TiffFormatError, TiffResult};
+
+use crate::coordinate_transform::CoordinateTransform;
+
+impl CoordinateTransform {
+    pub fn from_transformation_matrix(transformation_matrix: [f64; 16]) -> TiffResult<Self> {
+        let transform = [
+            transformation_matrix[0],
+            transformation_matrix[1],
+            transformation_matrix[3],
+            transformation_matrix[4],
+            transformation_matrix[5],
+            transformation_matrix[7],
+        ];
+
+        let det = transform[0] * transform[4] - transform[1] * transform[3];
+        if det.abs() < 0.000000000000001 {
+            return Err(TiffError::FormatError(TiffFormatError::Format(
+                "Provided transformation matrix is not invertible".into(),
+            )));
+        }
+
+        let inverse_transform = [
+            transform[4] / det,
+            -transform[1] / det,
+            (transform[1] * transform[5] - transform[2] * transform[4]) / det,
+            -transform[3] / det,
+            transform[0] / det,
+            (-transform[0] * transform[5] + transform[2] * transform[3]) / det,
+        ];
+
+        Ok(CoordinateTransform::AffineTransform {
+            transform,
+            inverse_transform,
+        })
+    }
+
+    pub(super) fn transform_by_affine_transform(transform: &[f64; 6], coord: &Coord) -> Coord {
+        Coord {
+            x: coord.x * transform[0] + coord.y * transform[1] + transform[2],
+            y: coord.x * transform[3] + coord.y * transform[4] + transform[5],
+        }
+    }
+}

src/coordinate_transform/tie_point_and_pixel_scale.rs

@@ -0,0 +1,50 @@
+use geo_types::Coord;
+use tiff::TiffResult;
+
+use crate::coordinate_transform::CoordinateTransform;
+
+impl CoordinateTransform {
+    pub(super) fn from_tie_point_and_pixel_scale(
+        tie_points: &[f64],
+        pixel_scale: &[f64],
+    ) -> TiffResult<Self> {
+        Ok(CoordinateTransform::TiePointAndPixelScale {
+            raster_point: Coord {
+                x: tie_points[0],
+                y: tie_points[1],
+            },
+            model_point: Coord {
+                x: tie_points[3],
+                y: tie_points[4],
+            },
+            pixel_scale: Coord {
+                x: pixel_scale[0],
+                y: pixel_scale[1],
+            },
+        })
+    }
+
+    pub(super) fn transform_to_model_by_tie_point_and_pixel_scale(
+        raster_point: &Coord,
+        model_point: &Coord,
+        pixel_scale: &Coord,
+        coord: &Coord,
+    ) -> Coord {
+        Coord {
+            x: (coord.x - raster_point.x) * pixel_scale.x + model_point.x,
+            y: (coord.y - raster_point.y) * -pixel_scale.y + model_point.y,
+        }
+    }
+
+    pub(super) fn transform_to_raster_by_tie_point_and_pixel_scale(
+        raster_point: &Coord,
+        model_point: &Coord,
+        pixel_scale: &Coord,
+        coord: &Coord,
+    ) -> Coord {
+        Coord {
+            x: (coord.x - model_point.x) / pixel_scale.x + raster_point.x,
+            y: (coord.y - model_point.y) / -pixel_scale.y + raster_point.y,
+        }
+    }
+}