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

97.63

/operators/src/processing/raster_vector_join/aggregated.rs

use futures::stream::BoxStream;
use futures::{StreamExt, TryStreamExt};

use geoengine_datatypes::collections::{
    FeatureCollection, FeatureCollectionInfos, FeatureCollectionModifications,
};
use geoengine_datatypes::raster::{GridIndexAccess, Pixel, RasterDataType};
use geoengine_datatypes::util::arrow::ArrowTyped;

use crate::engine::{
    QueryContext, QueryProcessor, RasterQueryProcessor, TypedRasterQueryProcessor,
    VectorQueryProcessor,
};
use crate::processing::raster_vector_join::aggregator::{
    Aggregator, FirstValueFloatAggregator, FirstValueIntAggregator, MeanValueAggregator,
    TypedAggregator,
};
use crate::processing::raster_vector_join::TemporalAggregationMethod;
use crate::util::Result;
use async_trait::async_trait;
use geoengine_datatypes::primitives::{BoundingBox2D, Geometry, VectorQueryRectangle};

use super::util::{CoveredPixels, FeatureTimeSpanIter, PixelCoverCreator};
use super::{create_feature_aggregator, FeatureAggregationMethod};

pub struct RasterVectorAggregateJoinProcessor<G> {
    collection: Box<dyn VectorQueryProcessor<VectorType = FeatureCollection<G>>>,
    raster_processors: Vec<TypedRasterQueryProcessor>,
    column_names: Vec<String>,
    feature_aggregation: FeatureAggregationMethod,
    temporal_aggregation: TemporalAggregationMethod,
}

impl<G> RasterVectorAggregateJoinProcessor<G>
where
    G: Geometry + ArrowTyped,
    FeatureCollection<G>: PixelCoverCreator<G>,
{
    pub fn new(
        collection: Box<dyn VectorQueryProcessor<VectorType = FeatureCollection<G>>>,
        raster_processors: Vec<TypedRasterQueryProcessor>,
        column_names: Vec<String>,
        feature_aggregation: FeatureAggregationMethod,
        temporal_aggregation: TemporalAggregationMethod,
    ) -> Self {
        Self {
            collection,
            raster_processors,
            column_names,
            feature_aggregation,
            temporal_aggregation,
        }
    }

    async fn extract_raster_values<P: Pixel>(
        collection: &FeatureCollection<G>,
        raster_processor: &dyn RasterQueryProcessor<RasterType = P>,
        new_column_name: &str,
        feature_aggreation: FeatureAggregationMethod,
        temporal_aggregation: TemporalAggregationMethod,
        query: VectorQueryRectangle,
        ctx: &dyn QueryContext,
    ) -> Result<FeatureCollection<G>> {
        let mut temporal_aggregator =
            Self::create_aggregator::<P>(collection.len(), temporal_aggregation);

        let collection = collection.sort_by_time_asc()?;

        let covered_pixels = collection.create_covered_pixels();

        let collection = covered_pixels.collection_ref();

        for time_span in FeatureTimeSpanIter::new(collection.time_intervals()) {
            let query = VectorQueryRectangle {
                spatial_bounds: query.spatial_bounds,
                time_interval: time_span.time_interval,
                spatial_resolution: query.spatial_resolution,
            };

            let mut rasters = raster_processor.raster_query(query.into(), ctx).await?;

            // TODO: optimize geo access (only specific tiles, etc.)

            let mut feature_aggregator =
                create_feature_aggregator::<P>(collection.len(), feature_aggreation);

            let mut time_end = None;

            while let Some(raster) = rasters.next().await {
                let raster = raster?;

                if let Some(end) = time_end {
                    if end != raster.time.end() {
                        // new time slice => consume old aggregator and create new one
                        temporal_aggregator.add_feature_data(
                            feature_aggregator.into_data(),
                            time_span.time_interval.duration_ms(), // TODO: use individual feature duration?
                        )?;

                        feature_aggregator =
                            create_feature_aggregator::<P>(collection.len(), feature_aggreation);

                        if temporal_aggregator.is_satisfied() {
                            break;
                        }
                    }
                }
                time_end = Some(raster.time.end());

                for feature_index in time_span.feature_index_start..=time_span.feature_index_end {
                    // TODO: don't do random access but use a single iterator
                    let mut satisfied = false;
                    for grid_idx in covered_pixels.covered_pixels(feature_index, &raster) {
                        // try to get the pixel if the coordinate is within the current tile
                        if let Ok(pixel) = raster.get_at_grid_index(grid_idx) {
                            // finally, attach value to feature
                            if let Some(data) = pixel {
                                feature_aggregator.add_value(feature_index, data, 1);
                            } else {
                                // TODO: weigh by area?
                                feature_aggregator.add_null(feature_index);
                            }

                            if feature_aggregator.is_satisfied() {
                                satisfied = true;
                                break;
                            }
                        }
                    }

                    if satisfied {
                        break;
                    }
                }
            }

            temporal_aggregator.add_feature_data(
                feature_aggregator.into_data(),
                time_span.time_interval.duration_ms(), // TODO: use individual feature duration?
            )?;

            if temporal_aggregator.is_satisfied() {
                break;
            }
        }

        collection
            .add_column(new_column_name, temporal_aggregator.into_data())
            .map_err(Into::into)
    }

    fn create_aggregator<P: Pixel>(
        number_of_features: usize,
        aggregation: TemporalAggregationMethod,
    ) -> TypedAggregator {
        match aggregation {
            TemporalAggregationMethod::First => match P::TYPE {
                RasterDataType::U8
                | RasterDataType::U16
                | RasterDataType::U32
                | RasterDataType::U64
                | RasterDataType::I8
                | RasterDataType::I16
                | RasterDataType::I32
                | RasterDataType::I64 => {
                    FirstValueIntAggregator::new(number_of_features).into_typed()
                }
                RasterDataType::F32 | RasterDataType::F64 => {
                    FirstValueFloatAggregator::new(number_of_features).into_typed()
                }
            },
            TemporalAggregationMethod::Mean => {
                MeanValueAggregator::new(number_of_features).into_typed()
            }
            TemporalAggregationMethod::None => {
                unreachable!("this type of aggregator does not lead to this kind of processor")
            }
        }
    }
}

#[async_trait]
impl<G> QueryProcessor for RasterVectorAggregateJoinProcessor<G>
where
    G: Geometry + ArrowTyped + 'static,
    FeatureCollection<G>: PixelCoverCreator<G>,
{
    type Output = FeatureCollection<G>;
    type SpatialBounds = BoundingBox2D;

    async fn _query<'a>(
        &'a self,
        query: VectorQueryRectangle,
        ctx: &'a dyn QueryContext,
    ) -> Result<BoxStream<'a, Result<Self::Output>>> {
        let stream = self.collection
            .query(query, ctx).await?
            .and_then(move |mut collection| async move {

                for (raster, new_column_name) in self.raster_processors.iter().zip(&self.column_names) {
                    collection = call_on_generic_raster_processor!(raster, raster => {
                        Self::extract_raster_values(&collection, raster, new_column_name, self.feature_aggregation, self.temporal_aggregation, query, ctx).await?
                    });
                }

                Ok(collection)
            })
            .boxed();

        Ok(stream)
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    use crate::engine::{ChunkByteSize, MockExecutionContext, RasterResultDescriptor};
    use crate::engine::{MockQueryContext, RasterOperator};
    use crate::mock::{MockRasterSource, MockRasterSourceParams};
    use geoengine_datatypes::collections::{MultiPointCollection, MultiPolygonCollection};
    use geoengine_datatypes::primitives::MultiPolygon;
    use geoengine_datatypes::raster::{Grid2D, RasterTile2D, TileInformation};
    use geoengine_datatypes::spatial_reference::SpatialReference;
    use geoengine_datatypes::util::test::TestDefault;
    use geoengine_datatypes::{
        primitives::{
            BoundingBox2D, FeatureDataRef, Measurement, MultiPoint, SpatialResolution, TimeInterval,
        },
        raster::TilingSpecification,
    };

    #[tokio::test]
    async fn extract_raster_values_single_raster() {
        let raster_tile = RasterTile2D::<u8>::new_with_tile_info(
            TimeInterval::default(),
            TileInformation {
                global_geo_transform: TestDefault::test_default(),
                global_tile_position: [0, 0].into(),
                tile_size_in_pixels: [3, 2].into(),
            },
            Grid2D::new([3, 2].into(), vec![1, 2, 3, 4, 5, 6])
                .unwrap()
                .into(),
        );

        let raster_source = MockRasterSource {
            params: MockRasterSourceParams {
                data: vec![raster_tile],
                result_descriptor: RasterResultDescriptor {
                    data_type: RasterDataType::U8,
                    spatial_reference: SpatialReference::epsg_4326().into(),
                    measurement: Measurement::Unitless,
                    time: None,
                    bbox: None,
                    resolution: None,
                },
            },
        }
        .boxed();

        let execution_context = MockExecutionContext::new_with_tiling_spec(
            TilingSpecification::new((0., 0.).into(), [3, 2].into()),
        );

        let raster_source = raster_source.initialize(&execution_context).await.unwrap();

        let points = MultiPointCollection::from_data(
            MultiPoint::many(vec![
                (0.0, 0.0),
                (1.0, 0.0),
                (0.0, -1.0),
                (1.0, -1.0),
                (0.0, -2.0),
                (1.0, -2.0),
            ])
            .unwrap(),
            vec![TimeInterval::default(); 6],
            Default::default(),
        )
        .unwrap();

        let result = RasterVectorAggregateJoinProcessor::extract_raster_values(
            &points,
            &raster_source.query_processor().unwrap().get_u8().unwrap(),
            "foo",
            FeatureAggregationMethod::First,
            TemporalAggregationMethod::First,
            VectorQueryRectangle {
                spatial_bounds: BoundingBox2D::new((0.0, -3.0).into(), (2.0, 0.).into()).unwrap(),
                time_interval: Default::default(),
                spatial_resolution: SpatialResolution::new(1., 1.).unwrap(),
            },
            &MockQueryContext::new(ChunkByteSize::MIN),
        )
        .await
        .unwrap();

        if let FeatureDataRef::Int(extracted_data) = result.data("foo").unwrap() {
            assert_eq!(extracted_data.as_ref(), &[1, 2, 3, 4, 5, 6]);
        } else {
            unreachable!();
        }
    }

    #[tokio::test]
    #[allow(clippy::float_cmp)]
    async fn extract_raster_values_two_raster_timesteps() {
        let raster_tile_a = RasterTile2D::<u8>::new_with_tile_info(
            TimeInterval::new(0, 10).unwrap(),
            TileInformation {
                global_geo_transform: TestDefault::test_default(),
                global_tile_position: [0, 0].into(),
                tile_size_in_pixels: [3, 2].into(),
            },
            Grid2D::new([3, 2].into(), vec![6, 5, 4, 3, 2, 1])
                .unwrap()
                .into(),
        );
        let raster_tile_b = RasterTile2D::new_with_tile_info(
            TimeInterval::new(10, 20).unwrap(),
            TileInformation {
                global_geo_transform: TestDefault::test_default(),
                global_tile_position: [0, 0].into(),
                tile_size_in_pixels: [3, 2].into(),
            },
            Grid2D::new([3, 2].into(), vec![1, 2, 3, 4, 5, 6])
                .unwrap()
                .into(),
        );

        let raster_source = MockRasterSource {
            params: MockRasterSourceParams {
                data: vec![raster_tile_a, raster_tile_b],
                result_descriptor: RasterResultDescriptor {
                    data_type: RasterDataType::U8,
                    spatial_reference: SpatialReference::epsg_4326().into(),
                    measurement: Measurement::Unitless,
                    time: None,
                    bbox: None,
                    resolution: None,
                },
            },
        }
        .boxed();

        let execution_context = MockExecutionContext::new_with_tiling_spec(
            TilingSpecification::new((0., 0.).into(), [3, 2].into()),
        );

        let raster_source = raster_source.initialize(&execution_context).await.unwrap();

        let points = MultiPointCollection::from_data(
            MultiPoint::many(vec![
                (0.0, 0.0),
                (1.0, 0.0),
                (0.0, -1.0),
                (1.0, -1.0),
                (0.0, -2.0),
                (1.0, -2.0),
            ])
            .unwrap(),
            vec![TimeInterval::default(); 6],
            Default::default(),
        )
        .unwrap();

        let result = RasterVectorAggregateJoinProcessor::extract_raster_values(
            &points,
            &raster_source.query_processor().unwrap().get_u8().unwrap(),
            "foo",
            FeatureAggregationMethod::First,
            TemporalAggregationMethod::Mean,
            VectorQueryRectangle {
                spatial_bounds: BoundingBox2D::new((0.0, -3.0).into(), (2.0, 0.0).into()).unwrap(),
                time_interval: Default::default(),
                spatial_resolution: SpatialResolution::new(1., 1.).unwrap(),
            },
            &MockQueryContext::new(ChunkByteSize::MIN),
        )
        .await
        .unwrap();

        if let FeatureDataRef::Float(extracted_data) = result.data("foo").unwrap() {
            assert_eq!(extracted_data.as_ref(), &[3.5, 3.5, 3.5, 3.5, 3.5, 3.5]);
        } else {
            unreachable!();
        }
    }

    #[tokio::test]
    #[allow(clippy::float_cmp)]
    async fn extract_raster_values_two_spatial_tiles_per_time_step() {
        let raster_tile_a_0 = RasterTile2D::<u8>::new_with_tile_info(
            TimeInterval::new(0, 10).unwrap(),
            TileInformation {
                global_geo_transform: TestDefault::test_default(),
                global_tile_position: [0, 0].into(),
                tile_size_in_pixels: [3, 2].into(),
            },
            Grid2D::new([3, 2].into(), vec![6, 5, 4, 3, 2, 1])
                .unwrap()
                .into(),
        );
        let raster_tile_a_1 = RasterTile2D::new_with_tile_info(
            TimeInterval::new(0, 10).unwrap(),
            TileInformation {
                global_geo_transform: TestDefault::test_default(),
                global_tile_position: [0, 1].into(),
                tile_size_in_pixels: [3, 2].into(),
            },
            Grid2D::new([3, 2].into(), vec![60, 50, 40, 30, 20, 10])
                .unwrap()
                .into(),
        );
        let raster_tile_b_0 = RasterTile2D::new_with_tile_info(
            TimeInterval::new(10, 20).unwrap(),
            TileInformation {
                global_geo_transform: TestDefault::test_default(),
                global_tile_position: [0, 0].into(),
                tile_size_in_pixels: [3, 2].into(),
            },
            Grid2D::new([3, 2].into(), vec![1, 2, 3, 4, 5, 6])
                .unwrap()
                .into(),
        );
        let raster_tile_b_1 = RasterTile2D::new_with_tile_info(
            TimeInterval::new(10, 20).unwrap(),
            TileInformation {
                global_geo_transform: TestDefault::test_default(),
                global_tile_position: [0, 1].into(),
                tile_size_in_pixels: [3, 2].into(),
            },
            Grid2D::new([3, 2].into(), vec![10, 20, 30, 40, 50, 60])
                .unwrap()
                .into(),
        );

        let raster_source = MockRasterSource {
            params: MockRasterSourceParams {
                data: vec![
                    raster_tile_a_0,
                    raster_tile_a_1,
                    raster_tile_b_0,
                    raster_tile_b_1,
                ],
                result_descriptor: RasterResultDescriptor {
                    data_type: RasterDataType::U8,
                    spatial_reference: SpatialReference::epsg_4326().into(),
                    measurement: Measurement::Unitless,
                    time: None,
                    bbox: None,
                    resolution: None,
                },
            },
        }
        .boxed();

        let execution_context = MockExecutionContext::new_with_tiling_spec(
            TilingSpecification::new((0., 0.).into(), [3, 2].into()),
        );

        let raster_source = raster_source.initialize(&execution_context).await.unwrap();

        let points = MultiPointCollection::from_data(
            MultiPoint::many(vec![
                vec![(0.0, 0.0), (2.0, 0.0)],
                vec![(1.0, 0.0), (3.0, 0.0)],
            ])
            .unwrap(),
            vec![TimeInterval::default(); 2],
            Default::default(),
        )
        .unwrap();

        let result = RasterVectorAggregateJoinProcessor::extract_raster_values(
            &points,
            &raster_source.query_processor().unwrap().get_u8().unwrap(),
            "foo",
            FeatureAggregationMethod::Mean,
            TemporalAggregationMethod::Mean,
            VectorQueryRectangle {
                spatial_bounds: BoundingBox2D::new((0.0, -3.0).into(), (4.0, 0.0).into()).unwrap(),
                time_interval: Default::default(),
                spatial_resolution: SpatialResolution::new(1., 1.).unwrap(),
            },
            &MockQueryContext::new(ChunkByteSize::MIN),
        )
        .await
        .unwrap();

        if let FeatureDataRef::Float(extracted_data) = result.data("foo").unwrap() {
            assert_eq!(
                extracted_data.as_ref(),
                &[(6. + 60. + 1. + 10.) / 4., (5. + 50. + 2. + 20.) / 4.]
            );
        } else {
            unreachable!();
        }
    }

    #[tokio::test]
    #[allow(clippy::too_many_lines)]
    #[allow(clippy::float_cmp)]
    async fn polygons() {
        let raster_tile_a_0 = RasterTile2D::<u8>::new_with_tile_info(
            TimeInterval::new(0, 10).unwrap(),
            TileInformation {
                global_geo_transform: TestDefault::test_default(),
                global_tile_position: [0, 0].into(),
                tile_size_in_pixels: [3, 2].into(),
            },
            Grid2D::new([3, 2].into(), vec![6, 5, 4, 3, 2, 1])
                .unwrap()
                .into(),
        );
        let raster_tile_a_1 = RasterTile2D::new_with_tile_info(
            TimeInterval::new(0, 10).unwrap(),
            TileInformation {
                global_geo_transform: TestDefault::test_default(),
                global_tile_position: [0, 1].into(),
                tile_size_in_pixels: [3, 2].into(),
            },
            Grid2D::new([3, 2].into(), vec![60, 50, 40, 30, 20, 10])
                .unwrap()
                .into(),
        );
        let raster_tile_b_0 = RasterTile2D::new_with_tile_info(
            TimeInterval::new(10, 20).unwrap(),
            TileInformation {
                global_geo_transform: TestDefault::test_default(),
                global_tile_position: [0, 0].into(),
                tile_size_in_pixels: [3, 2].into(),
            },
            Grid2D::new([3, 2].into(), vec![1, 2, 3, 4, 5, 6])
                .unwrap()
                .into(),
        );
        let raster_tile_b_1 = RasterTile2D::new_with_tile_info(
            TimeInterval::new(10, 20).unwrap(),
            TileInformation {
                global_geo_transform: TestDefault::test_default(),
                global_tile_position: [0, 1].into(),
                tile_size_in_pixels: [3, 2].into(),
            },
            Grid2D::new([3, 2].into(), vec![10, 20, 30, 40, 50, 60])
                .unwrap()
                .into(),
        );

        let raster_source = MockRasterSource {
            params: MockRasterSourceParams {
                data: vec![
                    raster_tile_a_0,
                    raster_tile_a_1,
                    raster_tile_b_0,
                    raster_tile_b_1,
                ],
                result_descriptor: RasterResultDescriptor {
                    data_type: RasterDataType::U8,
                    spatial_reference: SpatialReference::epsg_4326().into(),
                    measurement: Measurement::Unitless,
                    time: None,
                    bbox: None,
                    resolution: None,
                },
            },
        }
        .boxed();

        let execution_context = MockExecutionContext::new_with_tiling_spec(
            TilingSpecification::new((0., 0.).into(), [3, 2].into()),
        );

        let raster_source = raster_source.initialize(&execution_context).await.unwrap();

        let polygons = MultiPolygonCollection::from_data(
            vec![MultiPolygon::new(vec![vec![vec![
                (0.5, -0.5).into(),
                (4., -1.).into(),
                (0.5, -2.5).into(),
                (0.5, -0.5).into(),
            ]]])
            .unwrap()],
            vec![TimeInterval::default(); 1],
            Default::default(),
        )
        .unwrap();

        let result = RasterVectorAggregateJoinProcessor::extract_raster_values(
            &polygons,
            &raster_source.query_processor().unwrap().get_u8().unwrap(),
            "foo",
            FeatureAggregationMethod::Mean,
            TemporalAggregationMethod::Mean,
            VectorQueryRectangle {
                spatial_bounds: BoundingBox2D::new((0.0, -3.0).into(), (4.0, 0.0).into()).unwrap(),
                time_interval: Default::default(),
                spatial_resolution: SpatialResolution::new(1., 1.).unwrap(),
            },
            &MockQueryContext::new(ChunkByteSize::MIN),
        )
        .await
        .unwrap();

        if let FeatureDataRef::Float(extracted_data) = result.data("foo").unwrap() {
            assert_eq!(
                extracted_data.as_ref(),
                &[(3. + 1. + 40. + 30. + 4. + 6. + 30. + 40.) / 8.]
            );
        } else {
            unreachable!();
        }
    }
}

1	use futures::stream::BoxStream;
2	use futures::{StreamExt, TryStreamExt};
3
4	use geoengine_datatypes::collections::{
5	FeatureCollection, FeatureCollectionInfos, FeatureCollectionModifications,
6	};
7	use geoengine_datatypes::raster::{GridIndexAccess, Pixel, RasterDataType};
8	use geoengine_datatypes::util::arrow::ArrowTyped;
9
10	use crate::engine::{
11	QueryContext, QueryProcessor, RasterQueryProcessor, TypedRasterQueryProcessor,
12	VectorQueryProcessor,
13	};
14	use crate::processing::raster_vector_join::aggregator::{
15	Aggregator, FirstValueFloatAggregator, FirstValueIntAggregator, MeanValueAggregator,
16	TypedAggregator,
17	};
18	use crate::processing::raster_vector_join::TemporalAggregationMethod;
19	use crate::util::Result;
20	use async_trait::async_trait;
21	use geoengine_datatypes::primitives::{BoundingBox2D, Geometry, VectorQueryRectangle};
22
23	use super::util::{CoveredPixels, FeatureTimeSpanIter, PixelCoverCreator};
24	use super::{create_feature_aggregator, FeatureAggregationMethod};
25
26	pub struct RasterVectorAggregateJoinProcessor<G> {
27	collection: Box<dyn VectorQueryProcessor<VectorType = FeatureCollection<G>>>,
28	raster_processors: Vec<TypedRasterQueryProcessor>,
29	column_names: Vec<String>,
30	feature_aggregation: FeatureAggregationMethod,
31	temporal_aggregation: TemporalAggregationMethod,
32	}
33
34	impl<G> RasterVectorAggregateJoinProcessor<G>
35	where
36	G: Geometry + ArrowTyped,
37	FeatureCollection<G>: PixelCoverCreator<G>,
38	{
39	pub fn new(	4✔
40	collection: Box<dyn VectorQueryProcessor<VectorType = FeatureCollection<G>>>,	4✔
41	raster_processors: Vec<TypedRasterQueryProcessor>,	4✔
42	column_names: Vec<String>,	4✔
43	feature_aggregation: FeatureAggregationMethod,	4✔
44	temporal_aggregation: TemporalAggregationMethod,	4✔
45	) -> Self {	4✔
46	Self {	4✔
47	collection,	4✔
48	raster_processors,	4✔
49	column_names,	4✔
50	feature_aggregation,	4✔
51	temporal_aggregation,	4✔
52	}	4✔
53	}	4✔
54
55	async fn extract_raster_values<P: Pixel>(	8✔
56	collection: &FeatureCollection<G>,	8✔
57	raster_processor: &dyn RasterQueryProcessor<RasterType = P>,	8✔
58	new_column_name: &str,	8✔
59	feature_aggreation: FeatureAggregationMethod,	8✔
60	temporal_aggregation: TemporalAggregationMethod,	8✔
61	query: VectorQueryRectangle,	8✔
62	ctx: &dyn QueryContext,	8✔
63	) -> Result<FeatureCollection<G>> {	8✔
64	let mut temporal_aggregator =	8✔
65	Self::create_aggregator::<P>(collection.len(), temporal_aggregation);	8✔
66
67	let collection = collection.sort_by_time_asc()?;	8✔
68
69	let covered_pixels = collection.create_covered_pixels();	8✔
70		8✔
71	let collection = covered_pixels.collection_ref();	8✔
72
73	for time_span in FeatureTimeSpanIter::new(collection.time_intervals()) {	8✔
74	let query = VectorQueryRectangle {	8✔
75	spatial_bounds: query.spatial_bounds,	8✔
76	time_interval: time_span.time_interval,	8✔
77	spatial_resolution: query.spatial_resolution,	8✔
78	};	8✔
79
80	let mut rasters = raster_processor.raster_query(query.into(), ctx).await?;	8✔
81
82	// TODO: optimize geo access (only specific tiles, etc.)
83
84	let mut feature_aggregator =	8✔
85	create_feature_aggregator::<P>(collection.len(), feature_aggreation);	8✔
86		8✔
87	let mut time_end = None;	8✔
88
89	while let Some(raster) = rasters.next().await {	171✔
90	let raster = raster?;	165✔
91
92	if let Some(end) = time_end {	165✔
93	if end != raster.time.end() {	157✔
94	// new time slice => consume old aggregator and create new one
95	temporal_aggregator.add_feature_data(	6✔
96	feature_aggregator.into_data(),	6✔
97	time_span.time_interval.duration_ms(), // TODO: use individual feature duration?	6✔
98	)?;	6✔
99
100	feature_aggregator =
101	create_feature_aggregator::<P>(collection.len(), feature_aggreation);	6✔
102		6✔
103	if temporal_aggregator.is_satisfied() {	6✔
104	break;	2✔
105	}	4✔
106	}	151✔
107	}	8✔
108	time_end = Some(raster.time.end());	163✔
109
110	for feature_index in time_span.feature_index_start..=time_span.feature_index_end {	566✔
111	// TODO: don't do random access but use a single iterator
112	let mut satisfied = false;	566✔
113	for grid_idx in covered_pixels.covered_pixels(feature_index, &raster) {	566✔
114	// try to get the pixel if the coordinate is within the current tile
115	if let Ok(pixel) = raster.get_at_grid_index(grid_idx) {	51✔
116	// finally, attach value to feature
117	if let Some(data) = pixel {	51✔
118	feature_aggregator.add_value(feature_index, data, 1);	47✔
119	} else {	47✔
120	// TODO: weigh by area?	4✔
121	feature_aggregator.add_null(feature_index);	4✔
122	}	4✔
123
124	if feature_aggregator.is_satisfied() {	51✔
125	satisfied = true;	8✔
126	break;	8✔
127	}	43✔
128	}	×
129	}
130
131	if satisfied {	566✔
132	break;	8✔
133	}	558✔
134	}
135	}
136
137	temporal_aggregator.add_feature_data(	8✔
138	feature_aggregator.into_data(),	8✔
139	time_span.time_interval.duration_ms(), // TODO: use individual feature duration?	8✔
140	)?;	8✔
141
142	if temporal_aggregator.is_satisfied() {	8✔
143	break;	4✔
144	}	4✔
145	}
146
147	collection	8✔
148	.add_column(new_column_name, temporal_aggregator.into_data())	8✔
149	.map_err(Into::into)	8✔
150	}	8✔
151
152	fn create_aggregator<P: Pixel>(	8✔
153	number_of_features: usize,	8✔
154	aggregation: TemporalAggregationMethod,	8✔
155	) -> TypedAggregator {	8✔
156	match aggregation {	8✔
157	TemporalAggregationMethod::First => match P::TYPE {	3✔
158	RasterDataType::U8
159	\| RasterDataType::U16
160	\| RasterDataType::U32
161	\| RasterDataType::U64
162	\| RasterDataType::I8
163	\| RasterDataType::I16
164	\| RasterDataType::I32
165	\| RasterDataType::I64 => {
166	FirstValueIntAggregator::new(number_of_features).into_typed()	3✔
167	}
168	RasterDataType::F32 \| RasterDataType::F64 => {
169	FirstValueFloatAggregator::new(number_of_features).into_typed()	×
170	}
171	},
172	TemporalAggregationMethod::Mean => {
173	MeanValueAggregator::new(number_of_features).into_typed()	5✔
174	}
175	TemporalAggregationMethod::None => {
176	unreachable!("this type of aggregator does not lead to this kind of processor")	×
177	}
178	}
179	}	8✔
180	}
181
182	#[async_trait]
183	impl<G> QueryProcessor for RasterVectorAggregateJoinProcessor<G>
184	where
185	G: Geometry + ArrowTyped + 'static,
186	FeatureCollection<G>: PixelCoverCreator<G>,
187	{
188	type Output = FeatureCollection<G>;
189	type SpatialBounds = BoundingBox2D;
190
191	async fn _query<'a>(	4✔
192	&'a self,	4✔
193	query: VectorQueryRectangle,	4✔
194	ctx: &'a dyn QueryContext,	4✔
195	) -> Result<BoxStream<'a, Result<Self::Output>>> {	4✔
196	let stream = self.collection	4✔
197	.query(query, ctx).await?	4✔
198	.and_then(move \|mut collection\| async move {	4✔
199
200	for (raster, new_column_name) in self.raster_processors.iter().zip(&self.column_names) {	4✔
201	collection = call_on_generic_raster_processor!(raster, raster => {	4✔
202	Self::extract_raster_values(&collection, raster, new_column_name, self.feature_aggregation, self.temporal_aggregation, query, ctx).await?	117✔
203	});
204	}
205
206	Ok(collection)	4✔
207	})	4✔
208	.boxed();	4✔
209		4✔
210	Ok(stream)	4✔
211	}	8✔
212	}
213
214	#[cfg(test)]
215	mod tests {
216	use super::*;
217
218	use crate::engine::{ChunkByteSize, MockExecutionContext, RasterResultDescriptor};
219	use crate::engine::{MockQueryContext, RasterOperator};
220	use crate::mock::{MockRasterSource, MockRasterSourceParams};
221	use geoengine_datatypes::collections::{MultiPointCollection, MultiPolygonCollection};
222	use geoengine_datatypes::primitives::MultiPolygon;
223	use geoengine_datatypes::raster::{Grid2D, RasterTile2D, TileInformation};
224	use geoengine_datatypes::spatial_reference::SpatialReference;
225	use geoengine_datatypes::util::test::TestDefault;
226	use geoengine_datatypes::{
227	primitives::{
228	BoundingBox2D, FeatureDataRef, Measurement, MultiPoint, SpatialResolution, TimeInterval,
229	},
230	raster::TilingSpecification,
231	};
232
233	#[tokio::test]	1✔
234	async fn extract_raster_values_single_raster() {	1✔
235	let raster_tile = RasterTile2D::<u8>::new_with_tile_info(	1✔
236	TimeInterval::default(),	1✔
237	TileInformation {	1✔
238	global_geo_transform: TestDefault::test_default(),	1✔
239	global_tile_position: [0, 0].into(),	1✔
240	tile_size_in_pixels: [3, 2].into(),	1✔
241	},	1✔
242	Grid2D::new([3, 2].into(), vec![1, 2, 3, 4, 5, 6])	1✔
243	.unwrap()	1✔
244	.into(),	1✔
245	);	1✔
246		1✔
247	let raster_source = MockRasterSource {	1✔
248	params: MockRasterSourceParams {	1✔
249	data: vec![raster_tile],	1✔
250	result_descriptor: RasterResultDescriptor {	1✔
251	data_type: RasterDataType::U8,	1✔
252	spatial_reference: SpatialReference::epsg_4326().into(),	1✔
253	measurement: Measurement::Unitless,	1✔
254	time: None,	1✔
255	bbox: None,	1✔
256	resolution: None,	1✔
257	},	1✔
258	},	1✔
259	}	1✔
260	.boxed();	1✔
261		1✔
262	let execution_context = MockExecutionContext::new_with_tiling_spec(	1✔
263	TilingSpecification::new((0., 0.).into(), [3, 2].into()),	1✔
264	);	1✔
265
266	let raster_source = raster_source.initialize(&execution_context).await.unwrap();	1✔
267		1✔
268	let points = MultiPointCollection::from_data(	1✔
269	MultiPoint::many(vec![	1✔
270	(0.0, 0.0),	1✔
271	(1.0, 0.0),	1✔
272	(0.0, -1.0),	1✔
273	(1.0, -1.0),	1✔
274	(0.0, -2.0),	1✔
275	(1.0, -2.0),	1✔
276	])	1✔
277	.unwrap(),	1✔
278	vec![TimeInterval::default(); 6],	1✔
279	Default::default(),	1✔
280	)	1✔
281	.unwrap();	1✔
282
283	let result = RasterVectorAggregateJoinProcessor::extract_raster_values(	1✔
284	&points,	1✔
285	&raster_source.query_processor().unwrap().get_u8().unwrap(),	1✔
286	"foo",	1✔
287	FeatureAggregationMethod::First,	1✔
288	TemporalAggregationMethod::First,	1✔
289	VectorQueryRectangle {	1✔
290	spatial_bounds: BoundingBox2D::new((0.0, -3.0).into(), (2.0, 0.).into()).unwrap(),	1✔
291	time_interval: Default::default(),	1✔
292	spatial_resolution: SpatialResolution::new(1., 1.).unwrap(),	1✔
293	},	1✔
294	&MockQueryContext::new(ChunkByteSize::MIN),	1✔
295	)	1✔
296	.await	×
297	.unwrap();	1✔
298
299	if let FeatureDataRef::Int(extracted_data) = result.data("foo").unwrap() {	1✔
300	assert_eq!(extracted_data.as_ref(), &[1, 2, 3, 4, 5, 6]);	1✔
301	} else {
302	unreachable!();	×
303	}
304	}
305
306	#[tokio::test]	1✔
307	#[allow(clippy::float_cmp)]
308	async fn extract_raster_values_two_raster_timesteps() {	1✔
309	let raster_tile_a = RasterTile2D::<u8>::new_with_tile_info(	1✔
310	TimeInterval::new(0, 10).unwrap(),	1✔
311	TileInformation {	1✔
312	global_geo_transform: TestDefault::test_default(),	1✔
313	global_tile_position: [0, 0].into(),	1✔
314	tile_size_in_pixels: [3, 2].into(),	1✔
315	},	1✔
316	Grid2D::new([3, 2].into(), vec![6, 5, 4, 3, 2, 1])	1✔
317	.unwrap()	1✔
318	.into(),	1✔
319	);	1✔
320	let raster_tile_b = RasterTile2D::new_with_tile_info(	1✔
321	TimeInterval::new(10, 20).unwrap(),	1✔
322	TileInformation {	1✔
323	global_geo_transform: TestDefault::test_default(),	1✔
324	global_tile_position: [0, 0].into(),	1✔
325	tile_size_in_pixels: [3, 2].into(),	1✔
326	},	1✔
327	Grid2D::new([3, 2].into(), vec![1, 2, 3, 4, 5, 6])	1✔
328	.unwrap()	1✔
329	.into(),	1✔
330	);	1✔
331		1✔
332	let raster_source = MockRasterSource {	1✔
333	params: MockRasterSourceParams {	1✔
334	data: vec![raster_tile_a, raster_tile_b],	1✔
335	result_descriptor: RasterResultDescriptor {	1✔
336	data_type: RasterDataType::U8,	1✔
337	spatial_reference: SpatialReference::epsg_4326().into(),	1✔
338	measurement: Measurement::Unitless,	1✔
339	time: None,	1✔
340	bbox: None,	1✔
341	resolution: None,	1✔
342	},	1✔
343	},	1✔
344	}	1✔
345	.boxed();	1✔
346		1✔
347	let execution_context = MockExecutionContext::new_with_tiling_spec(	1✔
348	TilingSpecification::new((0., 0.).into(), [3, 2].into()),	1✔
349	);	1✔
350
351	let raster_source = raster_source.initialize(&execution_context).await.unwrap();	1✔
352		1✔
353	let points = MultiPointCollection::from_data(	1✔
354	MultiPoint::many(vec![	1✔
355	(0.0, 0.0),	1✔
356	(1.0, 0.0),	1✔
357	(0.0, -1.0),	1✔
358	(1.0, -1.0),	1✔
359	(0.0, -2.0),	1✔
360	(1.0, -2.0),	1✔
361	])	1✔
362	.unwrap(),	1✔
363	vec![TimeInterval::default(); 6],	1✔
364	Default::default(),	1✔
365	)	1✔
366	.unwrap();	1✔
367
368	let result = RasterVectorAggregateJoinProcessor::extract_raster_values(	1✔
369	&points,	1✔
370	&raster_source.query_processor().unwrap().get_u8().unwrap(),	1✔
371	"foo",	1✔
372	FeatureAggregationMethod::First,	1✔
373	TemporalAggregationMethod::Mean,	1✔
374	VectorQueryRectangle {	1✔
375	spatial_bounds: BoundingBox2D::new((0.0, -3.0).into(), (2.0, 0.0).into()).unwrap(),	1✔
376	time_interval: Default::default(),	1✔
377	spatial_resolution: SpatialResolution::new(1., 1.).unwrap(),	1✔
378	},	1✔
379	&MockQueryContext::new(ChunkByteSize::MIN),	1✔
380	)	1✔
381	.await	×
382	.unwrap();	1✔
383
384	if let FeatureDataRef::Float(extracted_data) = result.data("foo").unwrap() {	1✔
385	assert_eq!(extracted_data.as_ref(), &[3.5, 3.5, 3.5, 3.5, 3.5, 3.5]);	1✔
386	} else {
387	unreachable!();	×
388	}
389	}
390
391	#[tokio::test]	1✔
392	#[allow(clippy::float_cmp)]
393	async fn extract_raster_values_two_spatial_tiles_per_time_step() {	1✔
394	let raster_tile_a_0 = RasterTile2D::<u8>::new_with_tile_info(	1✔
395	TimeInterval::new(0, 10).unwrap(),	1✔
396	TileInformation {	1✔
397	global_geo_transform: TestDefault::test_default(),	1✔
398	global_tile_position: [0, 0].into(),	1✔
399	tile_size_in_pixels: [3, 2].into(),	1✔
400	},	1✔
401	Grid2D::new([3, 2].into(), vec![6, 5, 4, 3, 2, 1])	1✔
402	.unwrap()	1✔
403	.into(),	1✔
404	);	1✔
405	let raster_tile_a_1 = RasterTile2D::new_with_tile_info(	1✔
406	TimeInterval::new(0, 10).unwrap(),	1✔
407	TileInformation {	1✔
408	global_geo_transform: TestDefault::test_default(),	1✔
409	global_tile_position: [0, 1].into(),	1✔
410	tile_size_in_pixels: [3, 2].into(),	1✔
411	},	1✔
412	Grid2D::new([3, 2].into(), vec![60, 50, 40, 30, 20, 10])	1✔
413	.unwrap()	1✔
414	.into(),	1✔
415	);	1✔
416	let raster_tile_b_0 = RasterTile2D::new_with_tile_info(	1✔
417	TimeInterval::new(10, 20).unwrap(),	1✔
418	TileInformation {	1✔
419	global_geo_transform: TestDefault::test_default(),	1✔
420	global_tile_position: [0, 0].into(),	1✔
421	tile_size_in_pixels: [3, 2].into(),	1✔
422	},	1✔
423	Grid2D::new([3, 2].into(), vec![1, 2, 3, 4, 5, 6])	1✔
424	.unwrap()	1✔
425	.into(),	1✔
426	);	1✔
427	let raster_tile_b_1 = RasterTile2D::new_with_tile_info(	1✔
428	TimeInterval::new(10, 20).unwrap(),	1✔
429	TileInformation {	1✔
430	global_geo_transform: TestDefault::test_default(),	1✔
431	global_tile_position: [0, 1].into(),	1✔
432	tile_size_in_pixels: [3, 2].into(),	1✔
433	},	1✔
434	Grid2D::new([3, 2].into(), vec![10, 20, 30, 40, 50, 60])	1✔
435	.unwrap()	1✔
436	.into(),	1✔
437	);	1✔
438		1✔
439	let raster_source = MockRasterSource {	1✔
440	params: MockRasterSourceParams {	1✔
441	data: vec![	1✔
442	raster_tile_a_0,	1✔
443	raster_tile_a_1,	1✔
444	raster_tile_b_0,	1✔
445	raster_tile_b_1,	1✔
446	],	1✔
447	result_descriptor: RasterResultDescriptor {	1✔
448	data_type: RasterDataType::U8,	1✔
449	spatial_reference: SpatialReference::epsg_4326().into(),	1✔
450	measurement: Measurement::Unitless,	1✔
451	time: None,	1✔
452	bbox: None,	1✔
453	resolution: None,	1✔
454	},	1✔
455	},	1✔
456	}	1✔
457	.boxed();	1✔
458		1✔
459	let execution_context = MockExecutionContext::new_with_tiling_spec(	1✔
460	TilingSpecification::new((0., 0.).into(), [3, 2].into()),	1✔
461	);	1✔
462
463	let raster_source = raster_source.initialize(&execution_context).await.unwrap();	1✔
464		1✔
465	let points = MultiPointCollection::from_data(	1✔
466	MultiPoint::many(vec![	1✔
467	vec![(0.0, 0.0), (2.0, 0.0)],	1✔
468	vec![(1.0, 0.0), (3.0, 0.0)],	1✔
469	])	1✔
470	.unwrap(),	1✔
471	vec![TimeInterval::default(); 2],	1✔
472	Default::default(),	1✔
473	)	1✔
474	.unwrap();	1✔
475
476	let result = RasterVectorAggregateJoinProcessor::extract_raster_values(	1✔
477	&points,	1✔
478	&raster_source.query_processor().unwrap().get_u8().unwrap(),	1✔
479	"foo",	1✔
480	FeatureAggregationMethod::Mean,	1✔
481	TemporalAggregationMethod::Mean,	1✔
482	VectorQueryRectangle {	1✔
483	spatial_bounds: BoundingBox2D::new((0.0, -3.0).into(), (4.0, 0.0).into()).unwrap(),	1✔
484	time_interval: Default::default(),	1✔
485	spatial_resolution: SpatialResolution::new(1., 1.).unwrap(),	1✔
486	},	1✔
487	&MockQueryContext::new(ChunkByteSize::MIN),	1✔
488	)	1✔
489	.await	×
490	.unwrap();	1✔
491
492	if let FeatureDataRef::Float(extracted_data) = result.data("foo").unwrap() {	1✔
493	assert_eq!(	1✔
494	extracted_data.as_ref(),	1✔
495	&[(6. + 60. + 1. + 10.) / 4., (5. + 50. + 2. + 20.) / 4.]	1✔
496	);	1✔
497	} else {
498	unreachable!();	×
499	}
500	}
501
502	#[tokio::test]	1✔
503	#[allow(clippy::too_many_lines)]
504	#[allow(clippy::float_cmp)]
505	async fn polygons() {	1✔
506	let raster_tile_a_0 = RasterTile2D::<u8>::new_with_tile_info(	1✔
507	TimeInterval::new(0, 10).unwrap(),	1✔
508	TileInformation {	1✔
509	global_geo_transform: TestDefault::test_default(),	1✔
510	global_tile_position: [0, 0].into(),	1✔
511	tile_size_in_pixels: [3, 2].into(),	1✔
512	},	1✔
513	Grid2D::new([3, 2].into(), vec![6, 5, 4, 3, 2, 1])	1✔
514	.unwrap()	1✔
515	.into(),	1✔
516	);	1✔
517	let raster_tile_a_1 = RasterTile2D::new_with_tile_info(	1✔
518	TimeInterval::new(0, 10).unwrap(),	1✔
519	TileInformation {	1✔
520	global_geo_transform: TestDefault::test_default(),	1✔
521	global_tile_position: [0, 1].into(),	1✔
522	tile_size_in_pixels: [3, 2].into(),	1✔
523	},	1✔
524	Grid2D::new([3, 2].into(), vec![60, 50, 40, 30, 20, 10])	1✔
525	.unwrap()	1✔
526	.into(),	1✔
527	);	1✔
528	let raster_tile_b_0 = RasterTile2D::new_with_tile_info(	1✔
529	TimeInterval::new(10, 20).unwrap(),	1✔
530	TileInformation {	1✔
531	global_geo_transform: TestDefault::test_default(),	1✔
532	global_tile_position: [0, 0].into(),	1✔
533	tile_size_in_pixels: [3, 2].into(),	1✔
534	},	1✔
535	Grid2D::new([3, 2].into(), vec![1, 2, 3, 4, 5, 6])	1✔
536	.unwrap()	1✔
537	.into(),	1✔
538	);	1✔
539	let raster_tile_b_1 = RasterTile2D::new_with_tile_info(	1✔
540	TimeInterval::new(10, 20).unwrap(),	1✔
541	TileInformation {	1✔
542	global_geo_transform: TestDefault::test_default(),	1✔
543	global_tile_position: [0, 1].into(),	1✔
544	tile_size_in_pixels: [3, 2].into(),	1✔
545	},	1✔
546	Grid2D::new([3, 2].into(), vec![10, 20, 30, 40, 50, 60])	1✔
547	.unwrap()	1✔
548	.into(),	1✔
549	);	1✔
550		1✔
551	let raster_source = MockRasterSource {	1✔
552	params: MockRasterSourceParams {	1✔
553	data: vec![	1✔
554	raster_tile_a_0,	1✔
555	raster_tile_a_1,	1✔
556	raster_tile_b_0,	1✔
557	raster_tile_b_1,	1✔
558	],	1✔
559	result_descriptor: RasterResultDescriptor {	1✔
560	data_type: RasterDataType::U8,	1✔
561	spatial_reference: SpatialReference::epsg_4326().into(),	1✔
562	measurement: Measurement::Unitless,	1✔
563	time: None,	1✔
564	bbox: None,	1✔
565	resolution: None,	1✔
566	},	1✔
567	},	1✔
568	}	1✔
569	.boxed();	1✔
570		1✔
571	let execution_context = MockExecutionContext::new_with_tiling_spec(	1✔
572	TilingSpecification::new((0., 0.).into(), [3, 2].into()),	1✔
573	);	1✔
574
575	let raster_source = raster_source.initialize(&execution_context).await.unwrap();	1✔
576		1✔
577	let polygons = MultiPolygonCollection::from_data(	1✔
578	vec![MultiPolygon::new(vec![vec![vec![	1✔
579	(0.5, -0.5).into(),	1✔
580	(4., -1.).into(),	1✔
581	(0.5, -2.5).into(),	1✔
582	(0.5, -0.5).into(),	1✔
583	]]])	1✔
584	.unwrap()],	1✔
585	vec![TimeInterval::default(); 1],	1✔
586	Default::default(),	1✔
587	)	1✔
588	.unwrap();	1✔
589
590	let result = RasterVectorAggregateJoinProcessor::extract_raster_values(	1✔
591	&polygons,	1✔
592	&raster_source.query_processor().unwrap().get_u8().unwrap(),	1✔
593	"foo",	1✔
594	FeatureAggregationMethod::Mean,	1✔
595	TemporalAggregationMethod::Mean,	1✔
596	VectorQueryRectangle {	1✔
597	spatial_bounds: BoundingBox2D::new((0.0, -3.0).into(), (4.0, 0.0).into()).unwrap(),	1✔
598	time_interval: Default::default(),	1✔
599	spatial_resolution: SpatialResolution::new(1., 1.).unwrap(),	1✔
600	},	1✔
601	&MockQueryContext::new(ChunkByteSize::MIN),	1✔
602	)	1✔
603	.await	×
604	.unwrap();	1✔
605
606	if let FeatureDataRef::Float(extracted_data) = result.data("foo").unwrap() {	1✔
607	assert_eq!(	1✔
608	extracted_data.as_ref(),	1✔
609	&[(3. + 1. + 40. + 30. + 4. + 6. + 30. + 40.) / 8.]	1✔
610	);	1✔
611	} else {
612	unreachable!();	×
613	}
614	}
615	}

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