5006008836

Build Type

push

github

Committed by GitHub

Commit Message

Merge #785 #787

Run Details

936 of 936 new or added lines in 50 files covered. (100.0%)

96010 of 107707 relevant lines covered (89.14%)

72676.46 hits per line

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

92.23

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

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

use crate::adapters::{
    FoldTileAccu, FoldTileAccuMut, RasterSubQueryAdapter, SubQueryTileAggregator,
};
use crate::engine::{
    CanonicOperatorName, ExecutionContext, InitializedRasterOperator, InitializedSources, Operator,
    OperatorName, QueryContext, QueryProcessor, RasterOperator, RasterQueryProcessor,
    RasterResultDescriptor, SingleRasterSource, TypedRasterQueryProcessor, WorkflowOperatorPath,
};
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};

#[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>,
        path: WorkflowOperatorPath,
        context: &dyn ExecutionContext,
    ) -> Result<Box<dyn InitializedRasterOperator>> {
        let name = CanonicOperatorName::from(&self);

        let initialized_sources = self.sources.initialize_sources(path, context).await?;
        let raster_source = initialized_sources.raster;
        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 {
            name,
            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 {
    name: CanonicOperatorName,
    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
    }

    fn canonic_name(&self) -> CanonicOperatorName {
        self.name.clone()
    }
}

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

geo-engine / geoengine / 5006008836

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