Skip to content

Commit

Permalink
Merge pull request #334 from geo-engine/raster_reprojection_outside
Browse files Browse the repository at this point in the history 
produce empty tiles when querying outside of area of use of projection
  • Loading branch information
@ChristianBeilschmidt
ChristianBeilschmidt authored Sep 10, 2021
2 parents 13aacd4 + b949edd commit 88c4ca3
Show file tree
Hide file tree
Showing 32 changed files with 729 additions and 287 deletions.
80 changes: 71 additions & 9 deletions datatypes/src/raster/geo_transform.rs
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,8 @@
use std::cmp::max;

use crate::primitives::{
AxisAlignedRectangle, Coordinate2D, SpatialPartition2D, SpatialResolution,
use crate::{
primitives::{AxisAlignedRectangle, Coordinate2D, SpatialPartition2D, SpatialResolution},
util::test::TestDefault,
};
use serde::{de, Deserialize, Deserializer, Serialize};

Expand All @@ -10,13 +11,16 @@ use super::{GridBoundingBox2D, GridIdx, GridIdx2D};
/// This is a typedef for the `GDAL GeoTransform`. It represents an affine transformation matrix.
pub type GdalGeoTransform = [f64; 6];

/// The `GeoTransform` is a more user friendly representation of the `GDAL GeoTransform` affine transformation matrix.
/// The `GeoTransform` specifies the relation between pixel coordinates and geographic coordinates.
/// In Geo Engine x pixel size is always postive and y pixel size is always negative. For raster tiles
/// the origin is always the upper left corner. In the global grid for the `TilingStrategy` the origin
/// is always located at (0, 0).
#[derive(Copy, Clone, PartialEq, Debug, Serialize, Deserialize)]
#[serde(rename_all = "camelCase")]
pub struct GeoTransform {
pub origin_coordinate: Coordinate2D,
pub x_pixel_size: f64,
pub y_pixel_size: f64,
x_pixel_size: f64,
y_pixel_size: f64,
}

impl GeoTransform {
Expand All @@ -32,6 +36,9 @@ impl GeoTransform {
///
#[inline]
pub fn new(origin_coordinate: Coordinate2D, x_pixel_size: f64, y_pixel_size: f64) -> Self {
debug_assert!(x_pixel_size > 0.0);
debug_assert!(y_pixel_size < 0.0);

Self {
origin_coordinate,
x_pixel_size,
Expand All @@ -55,13 +62,24 @@ impl GeoTransform {
origin_coordinate_y: f64,
y_pixel_size: f64,
) -> Self {
debug_assert!(x_pixel_size > 0.0);
debug_assert!(y_pixel_size < 0.0);

Self {
origin_coordinate: (origin_coordinate_x, origin_coordinate_y).into(),
x_pixel_size,
y_pixel_size,
}
}

pub fn x_pixel_size(&self) -> f64 {
self.x_pixel_size
}

pub fn y_pixel_size(&self) -> f64 {
self.y_pixel_size
}

/// Transforms a grid coordinate (row, column) ~ (y, x) into a SRS coordinate (x,y)
/// The resulting coordinate is the upper left coordinate of the pixel
/// See GDAL documentation for more details (including the two ignored parameters): <https://gdal.org/user/raster_data_model.html>
Expand Down Expand Up @@ -108,8 +126,10 @@ impl GeoTransform {
///
#[inline]
pub fn coordinate_to_grid_idx_2d(&self, coord: Coordinate2D) -> GridIdx2D {
let grid_x_index = ((coord.x - self.origin_coordinate.x) / self.x_pixel_size) as isize;
let grid_y_index = ((coord.y - self.origin_coordinate.y) / self.y_pixel_size) as isize;
let grid_x_index =
((coord.x - self.origin_coordinate.x) / self.x_pixel_size).floor() as isize;
let grid_y_index =
((coord.y - self.origin_coordinate.y) / self.y_pixel_size).floor() as isize;
[grid_y_index, grid_x_index].into()
}

Expand Down Expand Up @@ -171,8 +191,8 @@ impl GeoTransform {
}
}

impl Default for GeoTransform {
fn default() -> Self {
impl TestDefault for GeoTransform {
fn test_default() -> Self {
GeoTransform::new_with_coordinate_x_y(0.0, 1.0, 0.0, -1.0)
}
}
Expand Down Expand Up @@ -273,6 +293,48 @@ mod tests {
);
}

#[test]
fn geo_transform_coordinate_2d_to_global_grid_2d() {
let geo_transform = GeoTransform::new_with_coordinate_x_y(0.0, 1.0, 0.0, -1.0);
assert_eq!(
geo_transform.coordinate_to_grid_idx_2d((0.0, 0.0).into()),
GridIdx2D::new([0, 0])
);
assert_eq!(
geo_transform.coordinate_to_grid_idx_2d((0.5, 0.0).into()),
GridIdx2D::new([0, 0])
);
assert_eq!(
geo_transform.coordinate_to_grid_idx_2d((0.5, 0.5).into()),
GridIdx2D::new([-1, 0])
);
assert_eq!(
geo_transform.coordinate_to_grid_idx_2d((0.0, 0.5).into()),
GridIdx2D::new([-1, 0])
);
assert_eq!(
geo_transform.coordinate_to_grid_idx_2d((0.5, -0.5).into()),
GridIdx2D::new([0, 0])
);
assert_eq!(
geo_transform.coordinate_to_grid_idx_2d((-0.5, 0.5).into()),
GridIdx2D::new([-1, -1])
);
assert_eq!(
geo_transform.coordinate_to_grid_idx_2d((-0.5, -0.5).into()),
GridIdx2D::new([0, -1])
);

assert_eq!(
geo_transform.coordinate_to_grid_idx_2d((1.5, -0.5).into()),
GridIdx2D::new([0, 1])
);
assert_eq!(
geo_transform.coordinate_to_grid_idx_2d((-1.5, 1.5).into()),
GridIdx2D::new([-2, -2])
);
}

#[test]
fn pixel_center() {
let geo_transform = GeoTransform::new_with_coordinate_x_y(5.0, 1.0, 5.0, -1.0);
Expand Down
67 changes: 48 additions & 19 deletions datatypes/src/raster/macros_raster_tile.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -243,6 +243,7 @@ mod tests {
use crate::{
primitives::TimeInterval,
raster::{GeoTransform, Grid2D, GridIndexAccess, Pixel, RasterTile2D, TypedRasterTile2D},
util::test::TestDefault,
};
use crate::{raster::RasterDataType, util::test::catch_unwind_silent};
use std::marker::PhantomData;
Expand All @@ -254,8 +255,11 @@ mod tests {
}

let r = Grid2D::new([3, 2].into(), vec![1, 2, 3, 4, 5, 6], None).unwrap();
let t =
RasterTile2D::new_without_offset(TimeInterval::default(), GeoTransform::default(), r);
let t = RasterTile2D::new_without_offset(
TimeInterval::default(),
GeoTransform::test_default(),
r,
);
let typed_raster = TypedRasterTile2D::U32(t);

assert_eq!(
Expand All @@ -271,8 +275,11 @@ mod tests {
}

let r = Grid2D::new([3, 2].into(), vec![1, 2, 3, 4, 5, 6], None).unwrap();
let t =
RasterTile2D::new_without_offset(TimeInterval::default(), GeoTransform::default(), r);
let t = RasterTile2D::new_without_offset(
TimeInterval::default(),
GeoTransform::test_default(),
r,
);
let typed_raster = TypedRasterTile2D::U32(t);

assert_eq!(
Expand All @@ -299,14 +306,18 @@ mod tests {
];

let r = Grid2D::new([3, 2].into(), data, None).unwrap();
RasterTile2D::new_without_offset(TimeInterval::default(), GeoTransform::default(), r)
RasterTile2D::new_without_offset(
TimeInterval::default(),
GeoTransform::test_default(),
r,
)
}

assert_eq!(
generate_generic_raster_tile_2d!(RasterDataType::U8, generate()),
TypedRasterTile2D::U8(RasterTile2D::new_without_offset(
TimeInterval::default(),
GeoTransform::default(),
GeoTransform::test_default(),
Grid2D::new([3, 2].into(), vec![1, 2, 3, 4, 5, 6], None,).unwrap()
),)
);
Expand All @@ -321,13 +332,19 @@ mod tests {
}

let r = Grid2D::new([3, 2].into(), vec![1, 2, 3, 4, 5, 6], None).unwrap();
let t =
RasterTile2D::new_without_offset(TimeInterval::default(), GeoTransform::default(), r);
let t = RasterTile2D::new_without_offset(
TimeInterval::default(),
GeoTransform::test_default(),
r,
);
let typed_raster_a = TypedRasterTile2D::U32(t);

let r = Grid2D::new([3, 2].into(), vec![1, 2, 3, 4, 5, 6], None).unwrap();
let t =
RasterTile2D::new_without_offset(TimeInterval::default(), GeoTransform::default(), r);
let t = RasterTile2D::new_without_offset(
TimeInterval::default(),
GeoTransform::test_default(),
r,
);
let typed_raster_b = TypedRasterTile2D::U16(t);

assert_eq!(
Expand All @@ -348,13 +365,19 @@ mod tests {
}

let r = Grid2D::new([3, 2].into(), vec![1, 2, 3, 4, 5, 6], None).unwrap();
let t =
RasterTile2D::new_without_offset(TimeInterval::default(), GeoTransform::default(), r);
let t = RasterTile2D::new_without_offset(
TimeInterval::default(),
GeoTransform::test_default(),
r,
);
let typed_raster_a = TypedRasterTile2D::U32(t);

let r = Grid2D::new([3, 2].into(), vec![1, 2, 3, 4, 5, 6], None).unwrap();
let t =
RasterTile2D::new_without_offset(TimeInterval::default(), GeoTransform::default(), r);
let t = RasterTile2D::new_without_offset(
TimeInterval::default(),
GeoTransform::test_default(),
r,
);
let typed_raster_b = TypedRasterTile2D::U16(t);

assert_eq!(
Expand Down Expand Up @@ -391,13 +414,19 @@ mod tests {
}

let r = Grid2D::new([3, 2].into(), vec![1, 2, 3, 4, 5, 6], None).unwrap();
let t =
RasterTile2D::new_without_offset(TimeInterval::default(), GeoTransform::default(), r);
let t = RasterTile2D::new_without_offset(
TimeInterval::default(),
GeoTransform::test_default(),
r,
);
let typed_raster_a = TypedRasterTile2D::U32(t);

let r = Grid2D::new([3, 2].into(), vec![1, 2, 3, 4, 5, 6], None).unwrap();
let t =
RasterTile2D::new_without_offset(TimeInterval::default(), GeoTransform::default(), r);
let t = RasterTile2D::new_without_offset(
TimeInterval::default(),
GeoTransform::test_default(),
r,
);
let typed_raster_b = TypedRasterTile2D::U16(t);

assert_eq!(
Expand Down Expand Up @@ -452,7 +481,7 @@ mod tests {
let typed_raster_a = TypedRasterTile2D::U32(RasterTile2D::new(
TimeInterval::default(),
[0, 0].into(),
[1.0, 1.0, 0.0, 1.0, 0.0, 1.0].into(),
[1.0, 1.0, 0.0, 1.0, 0.0, -1.0].into(),
Grid2D::new([3, 2].into(), vec![1, 2, 3, 4, 5, 6], None)
.unwrap()
.into(),
Expand Down
4 changes: 2 additions & 2 deletions datatypes/src/raster/operations/blit.rs
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,8 @@ impl<T: Pixel> Blit<RasterTile2D<T>> for MaterializedRasterTile2D<T> {
// TODO: ensure pixels are aligned

ensure!(
(self.geo_transform().x_pixel_size == source.geo_transform().x_pixel_size)
&& (self.geo_transform().y_pixel_size == source.geo_transform().y_pixel_size),
(self.geo_transform().x_pixel_size() == source.geo_transform().x_pixel_size())
&& (self.geo_transform().y_pixel_size() == source.geo_transform().y_pixel_size()),
error::Blit {
details: "Incompatible pixel size"
}
Expand Down
Loading

0 comments on commit 88c4ca3

Please sign in to comment.