3929938005

Build Type

push

github

Committed by GitHub

Commit Message

Merge #713

Run Details

84930 of 96741 relevant lines covered (87.79%)

79640.1 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

92.08

/operators/src/processing/interpolation/mod.rs

use std::marker::PhantomData;
use std::sync::Arc;

use crate::adapters::{
    FoldTileAccu, FoldTileAccuMut, RasterSubQueryAdapter, SubQueryTileAggregator,
};
use crate::engine::{
    CreateSpan, ExecutionContext, InitializedRasterOperator, Operator, OperatorName, QueryContext,
    QueryProcessor, RasterOperator, RasterQueryProcessor, RasterResultDescriptor,
    SingleRasterSource, TypedRasterQueryProcessor,
};
use crate::util::Result;
use async_trait::async_trait;
use futures::future::BoxFuture;
use futures::stream::BoxStream;
use futures::{Future, FutureExt, TryFuture, TryFutureExt};
use geoengine_datatypes::primitives::{
    AxisAlignedRectangle, Coordinate2D, RasterQueryRectangle, SpatialPartition2D,
    SpatialPartitioned, SpatialResolution, TimeInstance, TimeInterval,
};
use geoengine_datatypes::raster::{
    Bilinear, Blit, EmptyGrid2D, GeoTransform, GridOrEmpty, GridSize, InterpolationAlgorithm,
    NearestNeighbor, Pixel, RasterTile2D, TileInformation, TilingSpecification,
};
use rayon::ThreadPool;
use serde::{Deserialize, Serialize};
use snafu::{ensure, Snafu};
use tracing::{span, Level};

#[derive(Debug, Serialize, Deserialize, Clone)]
#[serde(rename_all = "camelCase")]
pub struct InterpolationParams {
    pub interpolation: InterpolationMethod,
    pub input_resolution: InputResolution,
}

#[derive(Debug, Serialize, Deserialize, Clone, Copy)]
#[serde(rename_all = "camelCase", tag = "type")]
pub enum InputResolution {
    Value(SpatialResolution),
    Source,
}

#[derive(Debug, Serialize, Deserialize, Clone, Copy)]
#[serde(rename_all = "camelCase")]
pub enum InterpolationMethod {
    NearestNeighbor,
    BiLinear,
}

#[derive(Debug, Snafu)]
#[snafu(visibility(pub(crate)), context(suffix(false)), module(error))]
pub enum InterpolationError {
    #[snafu(display(
        "The input resolution was defined as `source` but the source resolution is unknown.",
    ))]
    UnknownInputResolution,
}

pub type Interpolation = Operator<InterpolationParams, SingleRasterSource>;

impl OperatorName for Interpolation {
    const TYPE_NAME: &'static str = "Interpolation";
}

#[typetag::serde]
#[async_trait]
impl RasterOperator for Interpolation {
    async fn _initialize(
        self: Box<Self>,
        context: &dyn ExecutionContext,
    ) -> Result<Box<dyn InitializedRasterOperator>> {
        let raster_source = self.sources.raster.initialize(context).await?;
        let in_descriptor = raster_source.result_descriptor();

        ensure!(
            matches!(self.params.input_resolution, InputResolution::Value(_))
                || in_descriptor.resolution.is_some(),
            error::UnknownInputResolution
        );

        let input_resolution = if let InputResolution::Value(res) = self.params.input_resolution {
            res
        } else {
            in_descriptor.resolution.expect("checked in ensure")
        };

        let out_descriptor = RasterResultDescriptor {
            spatial_reference: in_descriptor.spatial_reference,
            data_type: in_descriptor.data_type,
            measurement: in_descriptor.measurement.clone(),
            bbox: in_descriptor.bbox,
            time: in_descriptor.time,
            resolution: None, // after interpolation the resolution is uncapped
        };

        let initialized_operator = InitializedInterpolation {
            result_descriptor: out_descriptor,
            raster_source,
            interpolation_method: self.params.interpolation,
            input_resolution,
            tiling_specification: context.tiling_specification(),
        };

        Ok(initialized_operator.boxed())
    }

    span_fn!(Interpolation);
}

pub struct InitializedInterpolation {
    result_descriptor: RasterResultDescriptor,
    raster_source: Box<dyn InitializedRasterOperator>,
    interpolation_method: InterpolationMethod,
    input_resolution: SpatialResolution,
    tiling_specification: TilingSpecification,
}

impl InitializedRasterOperator for InitializedInterpolation {
    fn query_processor(&self) -> Result<TypedRasterQueryProcessor> {
        let source_processor = self.raster_source.query_processor()?;

        let res = call_on_generic_raster_processor!(
            source_processor, p => match self.interpolation_method  {
                InterpolationMethod::NearestNeighbor => InterploationProcessor::<_,_, NearestNeighbor>::new(
                        p,
                        self.input_resolution,
                        self.tiling_specification,
                    ).boxed()
                    .into(),
                InterpolationMethod::BiLinear =>InterploationProcessor::<_,_, Bilinear>::new(
                        p,
                        self.input_resolution,
                        self.tiling_specification,
                    ).boxed()
                    .into(),
            }
        );

        Ok(res)
    }

    fn result_descriptor(&self) -> &RasterResultDescriptor {
        &self.result_descriptor
    }
}

pub struct InterploationProcessor<Q, P, I>
where
    Q: RasterQueryProcessor<RasterType = P>,
    P: Pixel,
    I: InterpolationAlgorithm<P>,
{
    source: Q,
    input_resolution: SpatialResolution,
    tiling_specification: TilingSpecification,
    interpolation: PhantomData<I>,
}

impl<Q, P, I> InterploationProcessor<Q, P, I>
where
    Q: RasterQueryProcessor<RasterType = P>,
    P: Pixel,
    I: InterpolationAlgorithm<P>,
{
    pub fn new(
        source: Q,
        input_resolution: SpatialResolution,
        tiling_specification: TilingSpecification,
    ) -> Self {
        Self {
            source,
            input_resolution,
            tiling_specification,
            interpolation: PhantomData,
        }
    }
}

#[async_trait]
impl<Q, P, I> QueryProcessor for InterploationProcessor<Q, P, I>
where
    Q: QueryProcessor<Output = RasterTile2D<P>, SpatialBounds = SpatialPartition2D>,
    P: Pixel,
    I: InterpolationAlgorithm<P>,
{
    type Output = RasterTile2D<P>;
    type SpatialBounds = SpatialPartition2D;

    async fn _query<'a>(
        &'a self,
        query: RasterQueryRectangle,
        ctx: &'a dyn QueryContext,
    ) -> Result<BoxStream<'a, Result<Self::Output>>> {
        // do not interpolate if the source resolution is already fine enough
        if query.spatial_resolution.x >= self.input_resolution.x
            && query.spatial_resolution.y >= self.input_resolution.y
        {
            // TODO: should we use the query or the input resolution here?
            return self.source.query(query, ctx).await;
        }

        let sub_query = InterpolationSubQuery::<_, P, I> {
            input_resolution: self.input_resolution,
            fold_fn: fold_future,
            tiling_specification: self.tiling_specification,
            phantom: PhantomData,
            _phantom_pixel_type: PhantomData,
        };

        Ok(RasterSubQueryAdapter::<'a, P, _, _>::new(
            &self.source,
            query,
            self.tiling_specification,
            ctx,
            sub_query,
        )
        .filter_and_fill())
    }
}

#[derive(Debug, Clone)]
pub struct InterpolationSubQuery<F, T, I> {
    input_resolution: SpatialResolution,
    fold_fn: F,
    tiling_specification: TilingSpecification,
    phantom: PhantomData<I>,
    _phantom_pixel_type: PhantomData<T>,
}

impl<'a, T, FoldM, FoldF, I> SubQueryTileAggregator<'a, T> for InterpolationSubQuery<FoldM, T, I>
where
    T: Pixel,
    FoldM: Send + Sync + 'a + Clone + Fn(InterpolationAccu<T, I>, RasterTile2D<T>) -> FoldF,
    FoldF: Send + TryFuture<Ok = InterpolationAccu<T, I>, Error = crate::error::Error>,
    I: InterpolationAlgorithm<T>,
{
    type FoldFuture = FoldF;

    type FoldMethod = FoldM;

    type TileAccu = InterpolationAccu<T, I>;
    type TileAccuFuture = BoxFuture<'a, Result<Self::TileAccu>>;

    fn new_fold_accu(
        &self,
        tile_info: TileInformation,
        query_rect: RasterQueryRectangle,
        pool: &Arc<ThreadPool>,
    ) -> Self::TileAccuFuture {
        create_accu(
            tile_info,
            query_rect,
            pool.clone(),
            self.tiling_specification,
        )
        .boxed()
    }

    fn tile_query_rectangle(
        &self,
        tile_info: TileInformation,
        _query_rect: RasterQueryRectangle,
        start_time: TimeInstance,
    ) -> Result<Option<RasterQueryRectangle>> {
        // enlarge the spatial bounds in order to have the neighbor pixels for the interpolation
        let spatial_bounds = tile_info.spatial_partition();
        let enlarge: Coordinate2D = (self.input_resolution.x, -self.input_resolution.y).into();
        let spatial_bounds = SpatialPartition2D::new(
            spatial_bounds.upper_left(),
            spatial_bounds.lower_right() + enlarge,
        )?;

        Ok(Some(RasterQueryRectangle {
            spatial_bounds,
            time_interval: TimeInterval::new_instant(start_time)?,
            spatial_resolution: self.input_resolution,
        }))
    }

    fn fold_method(&self) -> Self::FoldMethod {
        self.fold_fn.clone()
    }
}

#[derive(Clone, Debug)]
pub struct InterpolationAccu<T: Pixel, I: InterpolationAlgorithm<T>> {
    pub output_info: TileInformation,
    pub input_tile: RasterTile2D<T>,
    pub pool: Arc<ThreadPool>,
    phantom: PhantomData<I>,
}

impl<T: Pixel, I: InterpolationAlgorithm<T>> InterpolationAccu<T, I> {
    pub fn new(
        input_tile: RasterTile2D<T>,
        output_info: TileInformation,
        pool: Arc<ThreadPool>,
    ) -> Self {
        InterpolationAccu {
            input_tile,
            output_info,
            pool,
            phantom: Default::default(),
        }
    }
}

#[async_trait]
impl<T: Pixel, I: InterpolationAlgorithm<T>> FoldTileAccu for InterpolationAccu<T, I> {
    type RasterType = T;

    async fn into_tile(self) -> Result<RasterTile2D<Self::RasterType>> {
        // now that we collected all the input tile pixels we perform the actual interpolation

        let output_tile = crate::util::spawn_blocking_with_thread_pool(self.pool, move || {
            I::interpolate(&self.input_tile, &self.output_info)
        })
        .await??;

        Ok(output_tile)
    }

    fn thread_pool(&self) -> &Arc<ThreadPool> {
        &self.pool
    }
}

impl<T: Pixel, I: InterpolationAlgorithm<T>> FoldTileAccuMut for InterpolationAccu<T, I> {
    fn tile_mut(&mut self) -> &mut RasterTile2D<T> {
        &mut self.input_tile
    }
}

pub fn create_accu<T: Pixel, I: InterpolationAlgorithm<T>>(
    tile_info: TileInformation,
    query_rect: RasterQueryRectangle,
    pool: Arc<ThreadPool>,
    tiling_specification: TilingSpecification,
) -> impl Future<Output = Result<InterpolationAccu<T, I>>> {
    // create an accumulator as a single tile that fits all the input tiles
    crate::util::spawn_blocking(move || {
        let tiling = tiling_specification.strategy(
            query_rect.spatial_resolution.x,
            -query_rect.spatial_resolution.y,
        );

        let origin_coordinate = tiling
            .tile_information_iterator(query_rect.spatial_bounds)
            .next()
            .expect("a query contains at least one tile")
            .spatial_partition()
            .upper_left();

        let geo_transform = GeoTransform::new(
            origin_coordinate,
            query_rect.spatial_resolution.x,
            -query_rect.spatial_resolution.y,
        );

        let bbox = tiling.tile_grid_box(query_rect.spatial_bounds);

        let shape = [
            bbox.axis_size_y() * tiling.tile_size_in_pixels.axis_size_y(),
            bbox.axis_size_x() * tiling.tile_size_in_pixels.axis_size_x(),
        ];

        // create a non-aligned (w.r.t. the tiling specification) grid by setting the origin to the top-left of the tile and the tile-index to [0, 0]
        let grid = EmptyGrid2D::new(shape.into());

        let input_tile = RasterTile2D::new(
            query_rect.time_interval,
            [0, 0].into(),
            geo_transform,
            GridOrEmpty::from(grid),
        );

        InterpolationAccu::new(input_tile, tile_info, pool)
    })
    .map_err(From::from)
}

pub fn fold_future<T, I>(
    accu: InterpolationAccu<T, I>,
    tile: RasterTile2D<T>,
) -> impl Future<Output = Result<InterpolationAccu<T, I>>>
where
    T: Pixel,
    I: InterpolationAlgorithm<T>,
{
    crate::util::spawn_blocking(|| fold_impl(accu, tile)).then(|x| async move {
        match x {
            Ok(r) => r,
            Err(e) => Err(e.into()),
        }
    })
}

pub fn fold_impl<T, I>(
    mut accu: InterpolationAccu<T, I>,
    tile: RasterTile2D<T>,
) -> Result<InterpolationAccu<T, I>>
where
    T: Pixel,
    I: InterpolationAlgorithm<T>,
{
    // get the time now because it is not known when the accu was created
    accu.input_tile.time = tile.time;

    // TODO: add a skip if both tiles are empty?

    // copy all input tiles into the accu to have all data for interpolation
    let mut accu_input_tile = accu.input_tile.into_materialized_tile();
    accu_input_tile.blit(tile)?;

    Ok(InterpolationAccu::new(
        accu_input_tile.into(),
        accu.output_info,
        accu.pool,
    ))
}

#[cfg(test)]
mod tests {
    use super::*;
    use futures::StreamExt;
    use geoengine_datatypes::{
        primitives::{
            Measurement, RasterQueryRectangle, SpatialPartition2D, SpatialResolution, TimeInterval,
        },
        raster::{
            Grid2D, GridOrEmpty, RasterDataType, RasterTile2D, TileInformation, TilingSpecification,
        },
        spatial_reference::SpatialReference,
        util::test::TestDefault,
    };

    use crate::{
        engine::{MockExecutionContext, MockQueryContext, RasterOperator, RasterResultDescriptor},
        mock::{MockRasterSource, MockRasterSourceParams},
    };

    #[tokio::test]
    async fn nearest_neighbor_operator() -> Result<()> {
        let exe_ctx = MockExecutionContext::new_with_tiling_spec(TilingSpecification::new(
            (0., 0.).into(),
            [2, 2].into(),
        ));

        let raster = make_raster();

        let operator = Interpolation {
            params: InterpolationParams {
                interpolation: InterpolationMethod::NearestNeighbor,
                input_resolution: InputResolution::Value(SpatialResolution::one()),
            },
            sources: SingleRasterSource { raster },
        }
        .boxed()
        .initialize(&exe_ctx)
        .await?;

        let processor = operator.query_processor()?.get_i8().unwrap();

        let query_rect = RasterQueryRectangle {
            spatial_bounds: SpatialPartition2D::new_unchecked((0., 2.).into(), (4., 0.).into()),
            time_interval: TimeInterval::new_unchecked(0, 20),
            spatial_resolution: SpatialResolution::zero_point_five(),
        };
        let query_ctx = MockQueryContext::test_default();

        let result_stream = processor.query(query_rect, &query_ctx).await?;

        let result: Vec<Result<RasterTile2D<i8>>> = result_stream.collect().await;
        let result = result.into_iter().collect::<Result<Vec<_>>>()?;

        let mut times: Vec<TimeInterval> = vec![TimeInterval::new_unchecked(0, 10); 8];
        times.append(&mut vec![TimeInterval::new_unchecked(10, 20); 8]);

        let data = vec![
            vec![1, 2, 5, 6],
            vec![2, 3, 6, 7],
            vec![3, 4, 7, 8],
            vec![4, 0, 8, 0],
            vec![5, 6, 0, 0],
            vec![6, 7, 0, 0],
            vec![7, 8, 0, 0],
            vec![8, 0, 0, 0],
            vec![8, 7, 4, 3],
            vec![7, 6, 3, 2],
            vec![6, 5, 2, 1],
            vec![5, 0, 1, 0],
            vec![4, 3, 0, 0],
            vec![3, 2, 0, 0],
            vec![2, 1, 0, 0],
            vec![1, 0, 0, 0],
        ];

        let valid = vec![
            vec![true; 4],
            vec![true; 4],
            vec![true; 4],
            vec![true, false, true, false],
            vec![true, true, false, false],
            vec![true, true, false, false],
            vec![true, true, false, false],
            vec![true, false, false, false],
            vec![true; 4],
            vec![true; 4],
            vec![true; 4],
            vec![true, false, true, false],
            vec![true, true, false, false],
            vec![true, true, false, false],
            vec![true, true, false, false],
            vec![true, false, false, false],
        ];

        for (i, tile) in result.into_iter().enumerate() {
            let tile = tile.into_materialized_tile();
            assert_eq!(tile.time, times[i]);
            assert_eq!(tile.grid_array.inner_grid.data, data[i]);
            assert_eq!(tile.grid_array.validity_mask.data, valid[i]);
        }

        Ok(())
    }

    fn make_raster() -> Box<dyn RasterOperator> {
        // test raster:
        // [0, 10)
        // || 1 | 2 || 3 | 4 ||
        // || 5 | 6 || 7 | 8 ||
        //
        // [10, 20)
        // || 8 | 7 || 6 | 5 ||
        // || 4 | 3 || 2 | 1 ||
        let raster_tiles = vec![
            RasterTile2D::<i8>::new_with_tile_info(
                TimeInterval::new_unchecked(0, 10),
                TileInformation {
                    global_tile_position: [-1, 0].into(),
                    tile_size_in_pixels: [2, 2].into(),
                    global_geo_transform: TestDefault::test_default(),
                },
                GridOrEmpty::from(Grid2D::new([2, 2].into(), vec![1, 2, 5, 6]).unwrap()),
            ),
            RasterTile2D::new_with_tile_info(
                TimeInterval::new_unchecked(0, 10),
                TileInformation {
                    global_tile_position: [-1, 1].into(),
                    tile_size_in_pixels: [2, 2].into(),
                    global_geo_transform: TestDefault::test_default(),
                },
                GridOrEmpty::from(Grid2D::new([2, 2].into(), vec![3, 4, 7, 8]).unwrap()),
            ),
            RasterTile2D::new_with_tile_info(
                TimeInterval::new_unchecked(10, 20),
                TileInformation {
                    global_tile_position: [-1, 0].into(),
                    tile_size_in_pixels: [2, 2].into(),
                    global_geo_transform: TestDefault::test_default(),
                },
                GridOrEmpty::from(Grid2D::new([2, 2].into(), vec![8, 7, 4, 3]).unwrap()),
            ),
            RasterTile2D::new_with_tile_info(
                TimeInterval::new_unchecked(10, 20),
                TileInformation {
                    global_tile_position: [-1, 1].into(),
                    tile_size_in_pixels: [2, 2].into(),
                    global_geo_transform: TestDefault::test_default(),
                },
                GridOrEmpty::from(Grid2D::new([2, 2].into(), vec![6, 5, 2, 1]).unwrap()),
            ),
        ];

        MockRasterSource {
            params: MockRasterSourceParams {
                data: raster_tiles,
                result_descriptor: RasterResultDescriptor {
                    data_type: RasterDataType::I8,
                    spatial_reference: SpatialReference::epsg_4326().into(),
                    measurement: Measurement::Unitless,
                    time: None,
                    bbox: None,
                    resolution: None,
                },
            },
        }
        .boxed()
    }
}

1	use std::marker::PhantomData;
2	use std::sync::Arc;
3
4	use crate::adapters::{
5	FoldTileAccu, FoldTileAccuMut, RasterSubQueryAdapter, SubQueryTileAggregator,
6	};
7	use crate::engine::{
8	CreateSpan, ExecutionContext, InitializedRasterOperator, Operator, OperatorName, QueryContext,
9	QueryProcessor, RasterOperator, RasterQueryProcessor, RasterResultDescriptor,
10	SingleRasterSource, TypedRasterQueryProcessor,
11	};
12	use crate::util::Result;
13	use async_trait::async_trait;
14	use futures::future::BoxFuture;
15	use futures::stream::BoxStream;
16	use futures::{Future, FutureExt, TryFuture, TryFutureExt};
17	use geoengine_datatypes::primitives::{
18	AxisAlignedRectangle, Coordinate2D, RasterQueryRectangle, SpatialPartition2D,
19	SpatialPartitioned, SpatialResolution, TimeInstance, TimeInterval,
20	};
21	use geoengine_datatypes::raster::{
22	Bilinear, Blit, EmptyGrid2D, GeoTransform, GridOrEmpty, GridSize, InterpolationAlgorithm,
23	NearestNeighbor, Pixel, RasterTile2D, TileInformation, TilingSpecification,
24	};
25	use rayon::ThreadPool;
26	use serde::{Deserialize, Serialize};
27	use snafu::{ensure, Snafu};
28	use tracing::{span, Level};
29
30	#[derive(Debug, Serialize, Deserialize, Clone)]	×
31	#[serde(rename_all = "camelCase")]
32	pub struct InterpolationParams {
33	pub interpolation: InterpolationMethod,
34	pub input_resolution: InputResolution,
35	}
36
37	#[derive(Debug, Serialize, Deserialize, Clone, Copy)]	×
38	#[serde(rename_all = "camelCase", tag = "type")]
39	pub enum InputResolution {
40	Value(SpatialResolution),
41	Source,
42	}
43
44	#[derive(Debug, Serialize, Deserialize, Clone, Copy)]	×
45	#[serde(rename_all = "camelCase")]
46	pub enum InterpolationMethod {
47	NearestNeighbor,
48	BiLinear,
49	}
50
51	#[derive(Debug, Snafu)]	×
52	#[snafu(visibility(pub(crate)), context(suffix(false)), module(error))]	×
53	pub enum InterpolationError {
54	#[snafu(display(
55	"The input resolution was defined as `source` but the source resolution is unknown.",
56	))]
57	UnknownInputResolution,
58	}
59
60	pub type Interpolation = Operator<InterpolationParams, SingleRasterSource>;
61
62	impl OperatorName for Interpolation {
63	const TYPE_NAME: &'static str = "Interpolation";
64	}
65
66	#[typetag::serde]	×
67	#[async_trait]
68	impl RasterOperator for Interpolation {
69	async fn _initialize(	1✔
70	self: Box<Self>,	1✔
71	context: &dyn ExecutionContext,	1✔
72	) -> Result<Box<dyn InitializedRasterOperator>> {	1✔
73	let raster_source = self.sources.raster.initialize(context).await?;	1✔
74	let in_descriptor = raster_source.result_descriptor();	1✔
75
76	ensure!(	1✔
77	matches!(self.params.input_resolution, InputResolution::Value(_))	1✔
78	\|\| in_descriptor.resolution.is_some(),	×
79	error::UnknownInputResolution	×
80	);
81
82	let input_resolution = if let InputResolution::Value(res) = self.params.input_resolution {	1✔
83	res	1✔
84	} else {
85	in_descriptor.resolution.expect("checked in ensure")	×
86	};
87
88	let out_descriptor = RasterResultDescriptor {	1✔
89	spatial_reference: in_descriptor.spatial_reference,	1✔
90	data_type: in_descriptor.data_type,	1✔
91	measurement: in_descriptor.measurement.clone(),	1✔
92	bbox: in_descriptor.bbox,	1✔
93	time: in_descriptor.time,	1✔
94	resolution: None, // after interpolation the resolution is uncapped	1✔
95	};	1✔
96		1✔
97	let initialized_operator = InitializedInterpolation {	1✔
98	result_descriptor: out_descriptor,	1✔
99	raster_source,	1✔
100	interpolation_method: self.params.interpolation,	1✔
101	input_resolution,	1✔
102	tiling_specification: context.tiling_specification(),	1✔
103	};	1✔
104		1✔
105	Ok(initialized_operator.boxed())	1✔
106	}	2✔
107
108	span_fn!(Interpolation);	×
109	}
110
111	pub struct InitializedInterpolation {
112	result_descriptor: RasterResultDescriptor,
113	raster_source: Box<dyn InitializedRasterOperator>,
114	interpolation_method: InterpolationMethod,
115	input_resolution: SpatialResolution,
116	tiling_specification: TilingSpecification,
117	}
118
119	impl InitializedRasterOperator for InitializedInterpolation {
120	fn query_processor(&self) -> Result<TypedRasterQueryProcessor> {	1✔
121	let source_processor = self.raster_source.query_processor()?;	1✔
122
123	let res = call_on_generic_raster_processor!(	1✔
124	source_processor, p => match self.interpolation_method {	1✔
125	InterpolationMethod::NearestNeighbor => InterploationProcessor::<_,_, NearestNeighbor>::new(
126	p,	1✔
127	self.input_resolution,	1✔
128	self.tiling_specification,	1✔
129	).boxed()	1✔
130	.into(),	1✔
131	InterpolationMethod::BiLinear =>InterploationProcessor::<_,_, Bilinear>::new(
132	p,	×
133	self.input_resolution,	×
134	self.tiling_specification,	×
135	).boxed()	×
136	.into(),	×
137	}
138	);
139
140	Ok(res)	1✔
141	}	1✔
142
143	fn result_descriptor(&self) -> &RasterResultDescriptor {	×
144	&self.result_descriptor	×
145	}	×
146	}
147
148	pub struct InterploationProcessor<Q, P, I>
149	where
150	Q: RasterQueryProcessor<RasterType = P>,
151	P: Pixel,
152	I: InterpolationAlgorithm<P>,
153	{
154	source: Q,
155	input_resolution: SpatialResolution,
156	tiling_specification: TilingSpecification,
157	interpolation: PhantomData<I>,
158	}
159
160	impl<Q, P, I> InterploationProcessor<Q, P, I>
161	where
162	Q: RasterQueryProcessor<RasterType = P>,
163	P: Pixel,
164	I: InterpolationAlgorithm<P>,
165	{
166	pub fn new(	1✔
167	source: Q,	1✔
168	input_resolution: SpatialResolution,	1✔
169	tiling_specification: TilingSpecification,	1✔
170	) -> Self {	1✔
171	Self {	1✔
172	source,	1✔
173	input_resolution,	1✔
174	tiling_specification,	1✔
175	interpolation: PhantomData,	1✔
176	}	1✔
177	}	1✔
178	}
179
180	#[async_trait]
181	impl<Q, P, I> QueryProcessor for InterploationProcessor<Q, P, I>
182	where
183	Q: QueryProcessor<Output = RasterTile2D<P>, SpatialBounds = SpatialPartition2D>,
184	P: Pixel,
185	I: InterpolationAlgorithm<P>,
186	{
187	type Output = RasterTile2D<P>;
188	type SpatialBounds = SpatialPartition2D;
189
190	async fn _query<'a>(	1✔
191	&'a self,	1✔
192	query: RasterQueryRectangle,	1✔
193	ctx: &'a dyn QueryContext,	1✔
194	) -> Result<BoxStream<'a, Result<Self::Output>>> {	1✔
195	// do not interpolate if the source resolution is already fine enough
196	if query.spatial_resolution.x >= self.input_resolution.x	1✔
197	&& query.spatial_resolution.y >= self.input_resolution.y	×
198	{
199	// TODO: should we use the query or the input resolution here?
200	return self.source.query(query, ctx).await;	×
201	}	1✔
202		1✔
203	let sub_query = InterpolationSubQuery::<_, P, I> {	1✔
204	input_resolution: self.input_resolution,	1✔
205	fold_fn: fold_future,	1✔
206	tiling_specification: self.tiling_specification,	1✔
207	phantom: PhantomData,	1✔
208	_phantom_pixel_type: PhantomData,	1✔
209	};	1✔
210		1✔
211	Ok(RasterSubQueryAdapter::<'a, P, _, _>::new(	1✔
212	&self.source,	1✔
213	query,	1✔
214	self.tiling_specification,	1✔
215	ctx,	1✔
216	sub_query,	1✔
217	)	1✔
218	.filter_and_fill())	1✔
219	}	2✔
220	}
221
222	#[derive(Debug, Clone)]	×
223	pub struct InterpolationSubQuery<F, T, I> {
224	input_resolution: SpatialResolution,
225	fold_fn: F,
226	tiling_specification: TilingSpecification,
227	phantom: PhantomData<I>,
228	_phantom_pixel_type: PhantomData<T>,
229	}
230
231	impl<'a, T, FoldM, FoldF, I> SubQueryTileAggregator<'a, T> for InterpolationSubQuery<FoldM, T, I>
232	where
233	T: Pixel,
234	FoldM: Send + Sync + 'a + Clone + Fn(InterpolationAccu<T, I>, RasterTile2D<T>) -> FoldF,
235	FoldF: Send + TryFuture<Ok = InterpolationAccu<T, I>, Error = crate::error::Error>,
236	I: InterpolationAlgorithm<T>,
237	{
238	type FoldFuture = FoldF;
239
240	type FoldMethod = FoldM;
241
242	type TileAccu = InterpolationAccu<T, I>;
243	type TileAccuFuture = BoxFuture<'a, Result<Self::TileAccu>>;
244
245	fn new_fold_accu(	16✔
246	&self,	16✔
247	tile_info: TileInformation,	16✔
248	query_rect: RasterQueryRectangle,	16✔
249	pool: &Arc<ThreadPool>,	16✔
250	) -> Self::TileAccuFuture {	16✔
251	create_accu(	16✔
252	tile_info,	16✔
253	query_rect,	16✔
254	pool.clone(),	16✔
255	self.tiling_specification,	16✔
256	)	16✔
257	.boxed()	16✔
258	}	16✔
259
260	fn tile_query_rectangle(	16✔
261	&self,	16✔
262	tile_info: TileInformation,	16✔
263	_query_rect: RasterQueryRectangle,	16✔
264	start_time: TimeInstance,	16✔
265	) -> Result<Option<RasterQueryRectangle>> {	16✔
266	// enlarge the spatial bounds in order to have the neighbor pixels for the interpolation	16✔
267	let spatial_bounds = tile_info.spatial_partition();	16✔
268	let enlarge: Coordinate2D = (self.input_resolution.x, -self.input_resolution.y).into();	16✔
269	let spatial_bounds = SpatialPartition2D::new(	16✔
270	spatial_bounds.upper_left(),	16✔
271	spatial_bounds.lower_right() + enlarge,	16✔
272	)?;	16✔
273
274	Ok(Some(RasterQueryRectangle {
275	spatial_bounds,	16✔
276	time_interval: TimeInterval::new_instant(start_time)?,	16✔
277	spatial_resolution: self.input_resolution,	16✔
278	}))
279	}	16✔
280
281	fn fold_method(&self) -> Self::FoldMethod {	16✔
282	self.fold_fn.clone()	16✔
283	}	16✔
284	}
285
286	#[derive(Clone, Debug)]	×
287	pub struct InterpolationAccu<T: Pixel, I: InterpolationAlgorithm<T>> {
288	pub output_info: TileInformation,
289	pub input_tile: RasterTile2D<T>,
290	pub pool: Arc<ThreadPool>,
291	phantom: PhantomData<I>,
292	}
293
294	impl<T: Pixel, I: InterpolationAlgorithm<T>> InterpolationAccu<T, I> {
295	pub fn new(	52✔
296	input_tile: RasterTile2D<T>,	52✔
297	output_info: TileInformation,	52✔
298	pool: Arc<ThreadPool>,	52✔
299	) -> Self {	52✔
300	InterpolationAccu {	52✔
301	input_tile,	52✔
302	output_info,	52✔
303	pool,	52✔
304	phantom: Default::default(),	52✔
305	}	52✔
306	}	52✔
307	}
308
309	#[async_trait]
310	impl<T: Pixel, I: InterpolationAlgorithm<T>> FoldTileAccu for InterpolationAccu<T, I> {
311	type RasterType = T;
312
313	async fn into_tile(self) -> Result<RasterTile2D<Self::RasterType>> {	16✔
314	// now that we collected all the input tile pixels we perform the actual interpolation
315
316	let output_tile = crate::util::spawn_blocking_with_thread_pool(self.pool, move \|\| {	16✔
317	I::interpolate(&self.input_tile, &self.output_info)	16✔
318	})	16✔
319	.await??;	16✔
320
321	Ok(output_tile)	16✔
322	}	32✔
323
324	fn thread_pool(&self) -> &Arc<ThreadPool> {	×
325	&self.pool	×
326	}	×
327	}
328
329	impl<T: Pixel, I: InterpolationAlgorithm<T>> FoldTileAccuMut for InterpolationAccu<T, I> {
330	fn tile_mut(&mut self) -> &mut RasterTile2D<T> {	×
331	&mut self.input_tile	×
332	}	×
333	}
334
335	pub fn create_accu<T: Pixel, I: InterpolationAlgorithm<T>>(	16✔
336	tile_info: TileInformation,	16✔
337	query_rect: RasterQueryRectangle,	16✔
338	pool: Arc<ThreadPool>,	16✔
339	tiling_specification: TilingSpecification,	16✔
340	) -> impl Future<Output = Result<InterpolationAccu<T, I>>> {	16✔
341	// create an accumulator as a single tile that fits all the input tiles	16✔
342	crate::util::spawn_blocking(move \|\| {	16✔
343	let tiling = tiling_specification.strategy(	16✔
344	query_rect.spatial_resolution.x,	16✔
345	-query_rect.spatial_resolution.y,	16✔
346	);	16✔
347		16✔
348	let origin_coordinate = tiling	16✔
349	.tile_information_iterator(query_rect.spatial_bounds)	16✔
350	.next()	16✔
351	.expect("a query contains at least one tile")	16✔
352	.spatial_partition()	16✔
353	.upper_left();	16✔
354		16✔
355	let geo_transform = GeoTransform::new(	16✔
356	origin_coordinate,	16✔
357	query_rect.spatial_resolution.x,	16✔
358	-query_rect.spatial_resolution.y,	16✔
359	);	16✔
360		16✔
361	let bbox = tiling.tile_grid_box(query_rect.spatial_bounds);	16✔
362		16✔
363	let shape = [	16✔
364	bbox.axis_size_y() * tiling.tile_size_in_pixels.axis_size_y(),	16✔
365	bbox.axis_size_x() * tiling.tile_size_in_pixels.axis_size_x(),	16✔
366	];	16✔
367		16✔
368	// create a non-aligned (w.r.t. the tiling specification) grid by setting the origin to the top-left of the tile and the tile-index to [0, 0]	16✔
369	let grid = EmptyGrid2D::new(shape.into());	16✔
370		16✔
371	let input_tile = RasterTile2D::new(	16✔
372	query_rect.time_interval,	16✔
373	[0, 0].into(),	16✔
374	geo_transform,	16✔
375	GridOrEmpty::from(grid),	16✔
376	);	16✔
377		16✔
378	InterpolationAccu::new(input_tile, tile_info, pool)	16✔
379	})	16✔
380	.map_err(From::from)	16✔
381	}	16✔
382
383	pub fn fold_future<T, I>(	36✔
384	accu: InterpolationAccu<T, I>,	36✔
385	tile: RasterTile2D<T>,	36✔
386	) -> impl Future<Output = Result<InterpolationAccu<T, I>>>	36✔
387	where	36✔
388	T: Pixel,	36✔
389	I: InterpolationAlgorithm<T>,	36✔
390	{	36✔
391	crate::util::spawn_blocking(\|\| fold_impl(accu, tile)).then(\|x\| async move {	36✔
392	match x {	36✔
393	Ok(r) => r,	36✔
394	Err(e) => Err(e.into()),	×
395	}
396	})	36✔
397	}	36✔
398
399	pub fn fold_impl<T, I>(	36✔
400	mut accu: InterpolationAccu<T, I>,	36✔
401	tile: RasterTile2D<T>,	36✔
402	) -> Result<InterpolationAccu<T, I>>	36✔
403	where	36✔
404	T: Pixel,	36✔
405	I: InterpolationAlgorithm<T>,	36✔
406	{	36✔
407	// get the time now because it is not known when the accu was created	36✔
408	accu.input_tile.time = tile.time;	36✔
409		36✔
410	// TODO: add a skip if both tiles are empty?	36✔
411		36✔
412	// copy all input tiles into the accu to have all data for interpolation	36✔
413	let mut accu_input_tile = accu.input_tile.into_materialized_tile();	36✔
414	accu_input_tile.blit(tile)?;	36✔
415
416	Ok(InterpolationAccu::new(	36✔
417	accu_input_tile.into(),	36✔
418	accu.output_info,	36✔
419	accu.pool,	36✔
420	))	36✔
421	}	36✔
422
423	#[cfg(test)]
424	mod tests {
425	use super::*;
426	use futures::StreamExt;
427	use geoengine_datatypes::{
428	primitives::{
429	Measurement, RasterQueryRectangle, SpatialPartition2D, SpatialResolution, TimeInterval,
430	},
431	raster::{
432	Grid2D, GridOrEmpty, RasterDataType, RasterTile2D, TileInformation, TilingSpecification,
433	},
434	spatial_reference::SpatialReference,
435	util::test::TestDefault,
436	};
437
438	use crate::{
439	engine::{MockExecutionContext, MockQueryContext, RasterOperator, RasterResultDescriptor},
440	mock::{MockRasterSource, MockRasterSourceParams},
441	};
442
443	#[tokio::test]	1✔
444	async fn nearest_neighbor_operator() -> Result<()> {	1✔
445	let exe_ctx = MockExecutionContext::new_with_tiling_spec(TilingSpecification::new(	1✔
446	(0., 0.).into(),	1✔
447	[2, 2].into(),	1✔
448	));	1✔
449		1✔
450	let raster = make_raster();	1✔
451
452	let operator = Interpolation {	1✔
453	params: InterpolationParams {	1✔
454	interpolation: InterpolationMethod::NearestNeighbor,	1✔
455	input_resolution: InputResolution::Value(SpatialResolution::one()),	1✔
456	},	1✔
457	sources: SingleRasterSource { raster },	1✔
458	}	1✔
459	.boxed()	1✔
460	.initialize(&exe_ctx)	1✔
461	.await?;	×
462
463	let processor = operator.query_processor()?.get_i8().unwrap();	1✔
464		1✔
465	let query_rect = RasterQueryRectangle {	1✔
466	spatial_bounds: SpatialPartition2D::new_unchecked((0., 2.).into(), (4., 0.).into()),	1✔
467	time_interval: TimeInterval::new_unchecked(0, 20),	1✔
468	spatial_resolution: SpatialResolution::zero_point_five(),	1✔
469	};	1✔
470	let query_ctx = MockQueryContext::test_default();	1✔
471
472	let result_stream = processor.query(query_rect, &query_ctx).await?;	1✔
473
474	let result: Vec<Result<RasterTile2D<i8>>> = result_stream.collect().await;	67✔
475	let result = result.into_iter().collect::<Result<Vec<_>>>()?;	1✔
476
477	let mut times: Vec<TimeInterval> = vec![TimeInterval::new_unchecked(0, 10); 8];	1✔
478	times.append(&mut vec![TimeInterval::new_unchecked(10, 20); 8]);	1✔
479		1✔
480	let data = vec![	1✔
481	vec![1, 2, 5, 6],	1✔
482	vec![2, 3, 6, 7],	1✔
483	vec![3, 4, 7, 8],	1✔
484	vec![4, 0, 8, 0],	1✔
485	vec![5, 6, 0, 0],	1✔
486	vec![6, 7, 0, 0],	1✔
487	vec![7, 8, 0, 0],	1✔
488	vec![8, 0, 0, 0],	1✔
489	vec![8, 7, 4, 3],	1✔
490	vec![7, 6, 3, 2],	1✔
491	vec![6, 5, 2, 1],	1✔
492	vec![5, 0, 1, 0],	1✔
493	vec![4, 3, 0, 0],	1✔
494	vec![3, 2, 0, 0],	1✔
495	vec![2, 1, 0, 0],	1✔
496	vec![1, 0, 0, 0],	1✔
497	];	1✔
498		1✔
499	let valid = vec![	1✔
500	vec![true; 4],	1✔
501	vec![true; 4],	1✔
502	vec![true; 4],	1✔
503	vec![true, false, true, false],	1✔
504	vec![true, true, false, false],	1✔
505	vec![true, true, false, false],	1✔
506	vec![true, true, false, false],	1✔
507	vec![true, false, false, false],	1✔
508	vec![true; 4],	1✔
509	vec![true; 4],	1✔
510	vec![true; 4],	1✔
511	vec![true, false, true, false],	1✔
512	vec![true, true, false, false],	1✔
513	vec![true, true, false, false],	1✔
514	vec![true, true, false, false],	1✔
515	vec![true, false, false, false],	1✔
516	];	1✔
517
518	for (i, tile) in result.into_iter().enumerate() {	16✔
519	let tile = tile.into_materialized_tile();	16✔
520	assert_eq!(tile.time, times[i]);	16✔
521	assert_eq!(tile.grid_array.inner_grid.data, data[i]);	16✔
522	assert_eq!(tile.grid_array.validity_mask.data, valid[i]);	16✔
523	}
524
525	Ok(())	1✔
526	}
527
528	fn make_raster() -> Box<dyn RasterOperator> {	1✔
529	// test raster:	1✔
530	// [0, 10)	1✔
531	// \|\| 1 \| 2 \|\| 3 \| 4 \|\|	1✔
532	// \|\| 5 \| 6 \|\| 7 \| 8 \|\|	1✔
533	//	1✔
534	// [10, 20)	1✔
535	// \|\| 8 \| 7 \|\| 6 \| 5 \|\|	1✔
536	// \|\| 4 \| 3 \|\| 2 \| 1 \|\|	1✔
537	let raster_tiles = vec![	1✔
538	RasterTile2D::<i8>::new_with_tile_info(	1✔
539	TimeInterval::new_unchecked(0, 10),	1✔
540	TileInformation {	1✔
541	global_tile_position: [-1, 0].into(),	1✔
542	tile_size_in_pixels: [2, 2].into(),	1✔
543	global_geo_transform: TestDefault::test_default(),	1✔
544	},	1✔
545	GridOrEmpty::from(Grid2D::new([2, 2].into(), vec![1, 2, 5, 6]).unwrap()),	1✔
546	),	1✔
547	RasterTile2D::new_with_tile_info(	1✔
548	TimeInterval::new_unchecked(0, 10),	1✔
549	TileInformation {	1✔
550	global_tile_position: [-1, 1].into(),	1✔
551	tile_size_in_pixels: [2, 2].into(),	1✔
552	global_geo_transform: TestDefault::test_default(),	1✔
553	},	1✔
554	GridOrEmpty::from(Grid2D::new([2, 2].into(), vec![3, 4, 7, 8]).unwrap()),	1✔
555	),	1✔
556	RasterTile2D::new_with_tile_info(	1✔
557	TimeInterval::new_unchecked(10, 20),	1✔
558	TileInformation {	1✔
559	global_tile_position: [-1, 0].into(),	1✔
560	tile_size_in_pixels: [2, 2].into(),	1✔
561	global_geo_transform: TestDefault::test_default(),	1✔
562	},	1✔
563	GridOrEmpty::from(Grid2D::new([2, 2].into(), vec![8, 7, 4, 3]).unwrap()),	1✔
564	),	1✔
565	RasterTile2D::new_with_tile_info(	1✔
566	TimeInterval::new_unchecked(10, 20),	1✔
567	TileInformation {	1✔
568	global_tile_position: [-1, 1].into(),	1✔
569	tile_size_in_pixels: [2, 2].into(),	1✔
570	global_geo_transform: TestDefault::test_default(),	1✔
571	},	1✔
572	GridOrEmpty::from(Grid2D::new([2, 2].into(), vec![6, 5, 2, 1]).unwrap()),	1✔
573	),	1✔
574	];	1✔
575		1✔
576	MockRasterSource {	1✔
577	params: MockRasterSourceParams {	1✔
578	data: raster_tiles,	1✔
579	result_descriptor: RasterResultDescriptor {	1✔
580	data_type: RasterDataType::I8,	1✔
581	spatial_reference: SpatialReference::epsg_4326().into(),	1✔
582	measurement: Measurement::Unitless,	1✔
583	time: None,	1✔
584	bbox: None,	1✔
585	resolution: None,	1✔
586	},	1✔
587	},	1✔
588	}	1✔
589	.boxed()	1✔
590	}	1✔
591	}

geo-engine / geoengine / 3929938005

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous