6510717572

Committed 13 Oct 2023 03:18PM UTC coverage: 89.501% (-0.03%) from 89.532%

Build # 6510717572

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #885 from geo-engine/permission_txs

check permissions in the same transaction as queries

Run Details

216 of 216 new or added lines in 5 files covered. (100.0%)

109513 of 122359 relevant lines covered (89.5%)

59349.76 hits per line

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

85.78

/operators/src/processing/raster_vector_join/mod.rs

mod aggregated;
mod aggregator;
mod non_aggregated;
mod util;

use crate::engine::{
    CanonicOperatorName, ExecutionContext, InitializedRasterOperator, InitializedVectorOperator,
    Operator, OperatorName, SingleVectorMultipleRasterSources, TypedVectorQueryProcessor,
    VectorColumnInfo, VectorOperator, VectorQueryProcessor, VectorResultDescriptor,
    WorkflowOperatorPath,
};
use crate::error::{self, Error};
use crate::processing::raster_vector_join::non_aggregated::RasterVectorJoinProcessor;
use crate::util::Result;

use crate::processing::raster_vector_join::aggregated::RasterVectorAggregateJoinProcessor;
use async_trait::async_trait;
use geoengine_datatypes::collections::VectorDataType;
use geoengine_datatypes::primitives::FeatureDataType;
use geoengine_datatypes::raster::{Pixel, RasterDataType};
use serde::{Deserialize, Serialize};
use snafu::ensure;

use self::aggregator::{
    Aggregator, FirstValueFloatAggregator, FirstValueIntAggregator, MeanValueAggregator,
    TypedAggregator,
};

/// An operator that attaches raster values to vector data
pub type RasterVectorJoin = Operator<RasterVectorJoinParams, SingleVectorMultipleRasterSources>;

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

const MAX_NUMBER_OF_RASTER_INPUTS: usize = 8;

/// The parameter spec for `RasterVectorJoin`
#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize)]
#[serde(rename_all = "camelCase")]
pub struct RasterVectorJoinParams {
    /// Each name reflects the output column of the join result.
    /// For each raster input, one name must be defined.
    pub names: Vec<String>,

    /// Specifies which method is used for aggregating values for a feature
    pub feature_aggregation: FeatureAggregationMethod,

    /// Whether NO DATA values should be ignored in aggregating the joined feature data
    /// `false` by default
    #[serde(default)]
    pub feature_aggregation_ignore_no_data: bool,

    /// Specifies which method is used for aggregating values over time
    pub temporal_aggregation: TemporalAggregationMethod,

    /// Whether NO DATA values should be ignored in aggregating the joined temporal data
    /// `false` by default
    #[serde(default)]
    pub temporal_aggregation_ignore_no_data: bool,
}

/// How to aggregate the values for the geometries inside a feature e.g.
/// the mean of all the raster values corresponding to the individual
/// points inside a `MultiPoint` feature.
#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize, Copy)]
#[serde(rename_all = "camelCase")]
pub enum FeatureAggregationMethod {
    First,
    Mean,
}

/// How to aggregate the values over time
/// If there are multiple rasters valid during the validity of a feature
/// the featuer is either split into multiple (None-aggregation) or the
/// values are aggreagated
#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize, Copy)]
#[serde(rename_all = "camelCase")]
pub enum TemporalAggregationMethod {
    None,
    First,
    Mean,
}

#[typetag::serde]
#[async_trait]
impl VectorOperator for RasterVectorJoin {
    async fn _initialize(
        mut self: Box<Self>,
        path: WorkflowOperatorPath,
        context: &dyn ExecutionContext,
    ) -> Result<Box<dyn InitializedVectorOperator>> {
        ensure!(
            (1..=MAX_NUMBER_OF_RASTER_INPUTS).contains(&self.sources.rasters.len()),
            error::InvalidNumberOfRasterInputs {
                expected: 1..MAX_NUMBER_OF_RASTER_INPUTS,
                found: self.sources.rasters.len()
            }
        );
        ensure!(
            self.sources.rasters.len() == self.params.names.len(),
            error::InvalidOperatorSpec {
                reason: "`rasters` must be of equal length as `names`"
            }
        );

        let name = CanonicOperatorName::from(&self);

        let vector_source = self
            .sources
            .vector
            .initialize(path.clone_and_append(0), context)
            .await?;

        let vector_rd = vector_source.result_descriptor();

        ensure!(
            vector_rd.data_type != VectorDataType::Data,
            error::InvalidType {
                expected: format!(
                    "{}, {} or {}",
                    VectorDataType::MultiPoint,
                    VectorDataType::MultiLineString,
                    VectorDataType::MultiPolygon
                ),
                found: VectorDataType::Data.to_string()
            },
        );

        let raster_sources = futures::future::try_join_all(
            self.sources
                .rasters
                .into_iter()
                .enumerate()
                .map(|(i, op)| op.initialize(path.clone_and_append(i as u8 + 1), context)),
        )
        .await?;

        let spatial_reference = vector_rd.spatial_reference;

        for other_spatial_reference in raster_sources
            .iter()
            .map(|source| source.result_descriptor().spatial_reference)
        {
            ensure!(
                spatial_reference == other_spatial_reference,
                crate::error::InvalidSpatialReference {
                    expected: spatial_reference,
                    found: other_spatial_reference,
                }
            );
        }

        let params = self.params;

        let result_descriptor = vector_rd.map_columns(|columns| {
            let mut columns = columns.clone();
            for (i, new_column_name) in params.names.iter().enumerate() {
                let feature_data_type = match params.temporal_aggregation {
                    TemporalAggregationMethod::First | TemporalAggregationMethod::None => {
                        match raster_sources[i].result_descriptor().data_type {
                            RasterDataType::U8
                            | RasterDataType::U16
                            | RasterDataType::U32
                            | RasterDataType::U64
                            | RasterDataType::I8
                            | RasterDataType::I16
                            | RasterDataType::I32
                            | RasterDataType::I64 => FeatureDataType::Int,
                            RasterDataType::F32 | RasterDataType::F64 => FeatureDataType::Float,
                        }
                    }
                    TemporalAggregationMethod::Mean => FeatureDataType::Float,
                };
                columns.insert(
                    new_column_name.clone(),
                    VectorColumnInfo {
                        data_type: feature_data_type,
                        measurement: raster_sources[i].result_descriptor().measurement.clone(),
                    },
                );
            }
            columns
        });

        Ok(InitializedRasterVectorJoin {
            name,
            result_descriptor,
            vector_source,
            raster_sources,
            state: params,
        }
        .boxed())
    }

    span_fn!(RasterVectorJoin);
}

pub struct InitializedRasterVectorJoin {
    name: CanonicOperatorName,
    result_descriptor: VectorResultDescriptor,
    vector_source: Box<dyn InitializedVectorOperator>,
    raster_sources: Vec<Box<dyn InitializedRasterOperator>>,
    state: RasterVectorJoinParams,
}

impl InitializedVectorOperator for InitializedRasterVectorJoin {
    fn result_descriptor(&self) -> &VectorResultDescriptor {
        &self.result_descriptor
    }

    fn query_processor(&self) -> Result<TypedVectorQueryProcessor> {
        let typed_raster_processors = self
            .raster_sources
            .iter()
            .map(InitializedRasterOperator::query_processor)
            .collect::<Result<Vec<_>>>()?;

        Ok(match self.vector_source.query_processor()? {
            TypedVectorQueryProcessor::Data(_) => unreachable!(),
            TypedVectorQueryProcessor::MultiPoint(points) => {
                TypedVectorQueryProcessor::MultiPoint(match self.state.temporal_aggregation {
                    TemporalAggregationMethod::None => RasterVectorJoinProcessor::new(
                        points,
                        typed_raster_processors,
                        self.state.names.clone(),
                        self.state.feature_aggregation,
                        self.state.feature_aggregation_ignore_no_data,
                    )
                    .boxed(),
                    TemporalAggregationMethod::First | TemporalAggregationMethod::Mean => {
                        RasterVectorAggregateJoinProcessor::new(
                            points,
                            typed_raster_processors,
                            self.state.names.clone(),
                            self.state.feature_aggregation,
                            self.state.feature_aggregation_ignore_no_data,
                            self.state.temporal_aggregation,
                            self.state.temporal_aggregation_ignore_no_data,
                        )
                        .boxed()
                    }
                })
            }
            TypedVectorQueryProcessor::MultiPolygon(polygons) => {
                TypedVectorQueryProcessor::MultiPolygon(match self.state.temporal_aggregation {
                    TemporalAggregationMethod::None => RasterVectorJoinProcessor::new(
                        polygons,
                        typed_raster_processors,
                        self.state.names.clone(),
                        self.state.feature_aggregation,
                        self.state.feature_aggregation_ignore_no_data,
                    )
                    .boxed(),
                    TemporalAggregationMethod::First | TemporalAggregationMethod::Mean => {
                        RasterVectorAggregateJoinProcessor::new(
                            polygons,
                            typed_raster_processors,
                            self.state.names.clone(),
                            self.state.feature_aggregation,
                            self.state.feature_aggregation_ignore_no_data,
                            self.state.temporal_aggregation,
                            self.state.temporal_aggregation_ignore_no_data,
                        )
                        .boxed()
                    }
                })
            }
            TypedVectorQueryProcessor::MultiLineString(_) => return Err(Error::NotYetImplemented),
        })
    }

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

pub fn create_feature_aggregator<P: Pixel>(
    number_of_features: usize,
    aggregation: FeatureAggregationMethod,
    ignore_no_data: bool,
) -> TypedAggregator {
    match aggregation {
        FeatureAggregationMethod::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, ignore_no_data).into_typed()
            }
            RasterDataType::F32 | RasterDataType::F64 => {
                FirstValueFloatAggregator::new(number_of_features, ignore_no_data).into_typed()
            }
        },
        FeatureAggregationMethod::Mean => {
            MeanValueAggregator::new(number_of_features, ignore_no_data).into_typed()
        }
    }
}

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

    use crate::engine::{
        ChunkByteSize, MockExecutionContext, MockQueryContext, QueryProcessor, RasterOperator,
    };
    use crate::mock::MockFeatureCollectionSource;
    use crate::source::{GdalSource, GdalSourceParameters};
    use crate::util::gdal::add_ndvi_dataset;
    use futures::StreamExt;
    use geoengine_datatypes::collections::{FeatureCollectionInfos, MultiPointCollection};
    use geoengine_datatypes::dataset::NamedData;
    use geoengine_datatypes::primitives::CacheHint;
    use geoengine_datatypes::primitives::{
        BoundingBox2D, DataRef, DateTime, FeatureDataRef, MultiPoint, SpatialResolution,
        TimeInterval, VectorQueryRectangle,
    };
    use geoengine_datatypes::spatial_reference::SpatialReference;
    use geoengine_datatypes::util::{gdal::hide_gdal_errors, test::TestDefault};
    use serde_json::json;

    #[test]
    fn serialization() {
        let raster_vector_join = RasterVectorJoin {
            params: RasterVectorJoinParams {
                names: ["foo", "bar"].iter().copied().map(str::to_string).collect(),
                feature_aggregation: FeatureAggregationMethod::First,
                feature_aggregation_ignore_no_data: false,
                temporal_aggregation: TemporalAggregationMethod::Mean,
                temporal_aggregation_ignore_no_data: false,
            },
            sources: SingleVectorMultipleRasterSources {
                vector: MockFeatureCollectionSource::<MultiPoint>::multiple(vec![]).boxed(),
                rasters: vec![],
            },
        };

        let serialized = json!({
            "type": "RasterVectorJoin",
            "params": {
                "names": ["foo", "bar"],
                "featureAggregation": "first",
                "temporalAggregation": "mean",
            },
            "sources": {
                "vector": {
                    "type": "MockFeatureCollectionSourceMultiPoint",
                    "params": {
                        "collections": [],
                        "spatialReference": "EPSG:4326",
                        "measurements": {},
                    }
                },
                "rasters": [],
            },
        })
        .to_string();

        let deserialized: RasterVectorJoin = serde_json::from_str(&serialized).unwrap();

        assert_eq!(deserialized.params, raster_vector_join.params);
    }

    fn ndvi_source(name: NamedData) -> Box<dyn RasterOperator> {
        let gdal_source = GdalSource {
            params: GdalSourceParameters { data: name },
        };

        gdal_source.boxed()
    }

    #[tokio::test]
    async fn ndvi_time_point() {
        let point_source = MockFeatureCollectionSource::single(
            MultiPointCollection::from_data(
                MultiPoint::many(vec![
                    (-13.95, 20.05),
                    (-14.05, 20.05),
                    (-13.95, 19.95),
                    (-14.05, 19.95),
                ])
                .unwrap(),
                vec![
                    TimeInterval::new(
                        DateTime::new_utc(2014, 1, 1, 0, 0, 0),
                        DateTime::new_utc(2014, 1, 1, 0, 0, 0),
                    )
                    .unwrap();
                    4
                ],
                Default::default(),
                CacheHint::default(),
            )
            .unwrap(),
        )
        .boxed();

        let mut exe_ctc = MockExecutionContext::test_default();
        let ndvi_id = add_ndvi_dataset(&mut exe_ctc);

        let operator = RasterVectorJoin {
            params: RasterVectorJoinParams {
                names: vec!["ndvi".to_string()],
                feature_aggregation: FeatureAggregationMethod::First,
                feature_aggregation_ignore_no_data: false,
                temporal_aggregation: TemporalAggregationMethod::First,
                temporal_aggregation_ignore_no_data: false,
            },
            sources: SingleVectorMultipleRasterSources {
                vector: point_source,
                rasters: vec![ndvi_source(ndvi_id.clone())],
            },
        };

        let operator = operator
            .boxed()
            .initialize(WorkflowOperatorPath::initialize_root(), &exe_ctc)
            .await
            .unwrap();

        let query_processor = operator.query_processor().unwrap().multi_point().unwrap();

        let result = query_processor
            .query(
                VectorQueryRectangle {
                    spatial_bounds: BoundingBox2D::new((-180., -90.).into(), (180., 90.).into())
                        .unwrap(),
                    time_interval: TimeInterval::default(),
                    spatial_resolution: SpatialResolution::new(0.1, 0.1).unwrap(),
                },
                &MockQueryContext::new(ChunkByteSize::MIN),
            )
            .await
            .unwrap()
            .map(Result::unwrap)
            .collect::<Vec<MultiPointCollection>>()
            .await;

        assert_eq!(result.len(), 1);

        let FeatureDataRef::Int(data) = result[0].data("ndvi").unwrap() else {
            unreachable!();
        };

        // these values are taken from loading the tiff in QGIS
        assert_eq!(data.as_ref(), &[54, 55, 51, 55]);
    }

    #[tokio::test]
    #[allow(clippy::float_cmp)]
    async fn ndvi_time_range() {
        let point_source = MockFeatureCollectionSource::single(
            MultiPointCollection::from_data(
                MultiPoint::many(vec![
                    (-13.95, 20.05),
                    (-14.05, 20.05),
                    (-13.95, 19.95),
                    (-14.05, 19.95),
                ])
                .unwrap(),
                vec![
                    TimeInterval::new(
                        DateTime::new_utc(2014, 1, 1, 0, 0, 0),
                        DateTime::new_utc(2014, 3, 1, 0, 0, 0),
                    )
                    .unwrap();
                    4
                ],
                Default::default(),
                CacheHint::default(),
            )
            .unwrap(),
        )
        .boxed();

        let mut exe_ctc = MockExecutionContext::test_default();
        let ndvi_id = add_ndvi_dataset(&mut exe_ctc);

        let operator = RasterVectorJoin {
            params: RasterVectorJoinParams {
                names: vec!["ndvi".to_string()],
                feature_aggregation: FeatureAggregationMethod::First,
                feature_aggregation_ignore_no_data: false,
                temporal_aggregation: TemporalAggregationMethod::Mean,
                temporal_aggregation_ignore_no_data: false,
            },
            sources: SingleVectorMultipleRasterSources {
                vector: point_source,
                rasters: vec![ndvi_source(ndvi_id.clone())],
            },
        };

        let operator = operator
            .boxed()
            .initialize(WorkflowOperatorPath::initialize_root(), &exe_ctc)
            .await
            .unwrap();

        let query_processor = operator.query_processor().unwrap().multi_point().unwrap();

        let result = query_processor
            .query(
                VectorQueryRectangle {
                    spatial_bounds: BoundingBox2D::new((-180., -90.).into(), (180., 90.).into())
                        .unwrap(),
                    time_interval: TimeInterval::default(),
                    spatial_resolution: SpatialResolution::new(0.1, 0.1).unwrap(),
                },
                &MockQueryContext::new(ChunkByteSize::MIN),
            )
            .await
            .unwrap()
            .map(Result::unwrap)
            .collect::<Vec<MultiPointCollection>>()
            .await;

        assert_eq!(result.len(), 1);

        let FeatureDataRef::Float(data) = result[0].data("ndvi").unwrap() else {
            unreachable!();
        };

        // these values are taken from loading the tiff in QGIS
        assert_eq!(
            data.as_ref(),
            &[
                (54. + 52.) / 2.,
                (55. + 55.) / 2.,
                (51. + 50.) / 2.,
                (55. + 53.) / 2.,
            ]
        );
    }

    #[tokio::test]
    #[allow(clippy::float_cmp)]
    async fn ndvi_with_default_time() {
        hide_gdal_errors();

        let point_source = MockFeatureCollectionSource::single(
            MultiPointCollection::from_data(
                MultiPoint::many(vec![
                    (-13.95, 20.05),
                    (-14.05, 20.05),
                    (-13.95, 19.95),
                    (-14.05, 19.95),
                ])
                .unwrap(),
                vec![TimeInterval::default(); 4],
                Default::default(),
                CacheHint::default(),
            )
            .unwrap(),
        )
        .boxed();

        let mut exe_ctc = MockExecutionContext::test_default();
        let ndvi_id = add_ndvi_dataset(&mut exe_ctc);

        let operator = RasterVectorJoin {
            params: RasterVectorJoinParams {
                names: vec!["ndvi".to_string()],
                feature_aggregation: FeatureAggregationMethod::First,
                feature_aggregation_ignore_no_data: false,
                temporal_aggregation: TemporalAggregationMethod::Mean,
                temporal_aggregation_ignore_no_data: false,
            },
            sources: SingleVectorMultipleRasterSources {
                vector: point_source,
                rasters: vec![ndvi_source(ndvi_id.clone())],
            },
        };

        let operator = operator
            .boxed()
            .initialize(WorkflowOperatorPath::initialize_root(), &exe_ctc)
            .await
            .unwrap();

        let query_processor = operator.query_processor().unwrap().multi_point().unwrap();

        let result = query_processor
            .query(
                VectorQueryRectangle {
                    spatial_bounds: BoundingBox2D::new((-180., -90.).into(), (180., 90.).into())
                        .unwrap(),
                    time_interval: TimeInterval::default(),
                    spatial_resolution: SpatialResolution::new(0.1, 0.1).unwrap(),
                },
                &MockQueryContext::new(ChunkByteSize::MIN),
            )
            .await
            .unwrap()
            .map(Result::unwrap)
            .collect::<Vec<MultiPointCollection>>()
            .await;

        assert_eq!(result.len(), 1);

        let FeatureDataRef::Float(data) = result[0].data("ndvi").unwrap() else {
            unreachable!();
        };

        assert_eq!(data.as_ref(), &[0., 0., 0., 0.]);

        assert_eq!(data.nulls(), vec![true, true, true, true]);
    }

    #[tokio::test]
    async fn it_checks_sref() {
        let point_source = MockFeatureCollectionSource::with_collections_and_sref(
            vec![MultiPointCollection::from_data(
                MultiPoint::many(vec![
                    (-13.95, 20.05),
                    (-14.05, 20.05),
                    (-13.95, 19.95),
                    (-14.05, 19.95),
                ])
                .unwrap(),
                vec![TimeInterval::default(); 4],
                Default::default(),
                CacheHint::default(),
            )
            .unwrap()],
            SpatialReference::from_str("EPSG:3857").unwrap(),
        )
        .boxed();

        let mut exe_ctc = MockExecutionContext::test_default();
        let ndvi_id = add_ndvi_dataset(&mut exe_ctc);

        let operator = RasterVectorJoin {
            params: RasterVectorJoinParams {
                names: vec!["ndvi".to_string()],
                feature_aggregation: FeatureAggregationMethod::First,
                feature_aggregation_ignore_no_data: false,
                temporal_aggregation: TemporalAggregationMethod::Mean,
                temporal_aggregation_ignore_no_data: false,
            },
            sources: SingleVectorMultipleRasterSources {
                vector: point_source,
                rasters: vec![ndvi_source(ndvi_id.clone())],
            },
        }
        .boxed()
        .initialize(WorkflowOperatorPath::initialize_root(), &exe_ctc)
        .await;

        assert!(matches!(
            operator,
            Err(Error::InvalidSpatialReference {
                expected: _,
                found: _,
            })
        ));
    }
}

1	mod aggregated;
2	mod aggregator;
3	mod non_aggregated;
4	mod util;
5
6	use crate::engine::{
7	CanonicOperatorName, ExecutionContext, InitializedRasterOperator, InitializedVectorOperator,
8	Operator, OperatorName, SingleVectorMultipleRasterSources, TypedVectorQueryProcessor,
9	VectorColumnInfo, VectorOperator, VectorQueryProcessor, VectorResultDescriptor,
10	WorkflowOperatorPath,
11	};
12	use crate::error::{self, Error};
13	use crate::processing::raster_vector_join::non_aggregated::RasterVectorJoinProcessor;
14	use crate::util::Result;
15
16	use crate::processing::raster_vector_join::aggregated::RasterVectorAggregateJoinProcessor;
17	use async_trait::async_trait;
18	use geoengine_datatypes::collections::VectorDataType;
19	use geoengine_datatypes::primitives::FeatureDataType;
20	use geoengine_datatypes::raster::{Pixel, RasterDataType};
21	use serde::{Deserialize, Serialize};
22	use snafu::ensure;
23
24	use self::aggregator::{
25	Aggregator, FirstValueFloatAggregator, FirstValueIntAggregator, MeanValueAggregator,
26	TypedAggregator,
27	};
28
29	/// An operator that attaches raster values to vector data
30	pub type RasterVectorJoin = Operator<RasterVectorJoinParams, SingleVectorMultipleRasterSources>;
31
32	impl OperatorName for RasterVectorJoin {
33	const TYPE_NAME: &'static str = "RasterVectorJoin";
34	}
35
36	const MAX_NUMBER_OF_RASTER_INPUTS: usize = 8;
37
38	/// The parameter spec for `RasterVectorJoin`
39	#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize)]	25✔
40	#[serde(rename_all = "camelCase")]
41	pub struct RasterVectorJoinParams {
42	/// Each name reflects the output column of the join result.
43	/// For each raster input, one name must be defined.
44	pub names: Vec<String>,
45
46	/// Specifies which method is used for aggregating values for a feature
47	pub feature_aggregation: FeatureAggregationMethod,
48
49	/// Whether NO DATA values should be ignored in aggregating the joined feature data
50	/// `false` by default
51	#[serde(default)]
52	pub feature_aggregation_ignore_no_data: bool,
53
54	/// Specifies which method is used for aggregating values over time
55	pub temporal_aggregation: TemporalAggregationMethod,
56
57	/// Whether NO DATA values should be ignored in aggregating the joined temporal data
58	/// `false` by default
59	#[serde(default)]
60	pub temporal_aggregation_ignore_no_data: bool,
61	}
62
63	/// How to aggregate the values for the geometries inside a feature e.g.
64	/// the mean of all the raster values corresponding to the individual
65	/// points inside a `MultiPoint` feature.
66	#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize, Copy)]	7✔
67	#[serde(rename_all = "camelCase")]
68	pub enum FeatureAggregationMethod {
69	First,
70	Mean,
71	}
72
73	/// How to aggregate the values over time
74	/// If there are multiple rasters valid during the validity of a feature
75	/// the featuer is either split into multiple (None-aggregation) or the
76	/// values are aggreagated
77	#[derive(Debug, Clone, PartialEq, Eq, Serialize, Deserialize, Copy)]	7✔
78	#[serde(rename_all = "camelCase")]
79	pub enum TemporalAggregationMethod {
80	None,
81	First,
82	Mean,
83	}
84
85	#[typetag::serde]	2✔
86	#[async_trait]
87	impl VectorOperator for RasterVectorJoin {
88	async fn _initialize(	5✔
89	mut self: Box<Self>,	5✔
90	path: WorkflowOperatorPath,	5✔
91	context: &dyn ExecutionContext,	5✔
92	) -> Result<Box<dyn InitializedVectorOperator>> {	5✔
93	ensure!(	5✔
94	(1..=MAX_NUMBER_OF_RASTER_INPUTS).contains(&self.sources.rasters.len()),	5✔
95	error::InvalidNumberOfRasterInputs {	×
96	expected: 1..MAX_NUMBER_OF_RASTER_INPUTS,	×
97	found: self.sources.rasters.len()	×
98	}	×
99	);
100	ensure!(	5✔
101	self.sources.rasters.len() == self.params.names.len(),	5✔
102	error::InvalidOperatorSpec {	×
103	reason: "`rasters` must be of equal length as `names`"	×
104	}	×
105	);
106
107	let name = CanonicOperatorName::from(&self);	5✔
108
109	let vector_source = self	5✔
110	.sources	5✔
111	.vector	5✔
112	.initialize(path.clone_and_append(0), context)	5✔
113	.await?;	×
114
115	let vector_rd = vector_source.result_descriptor();	5✔
116		5✔
117	ensure!(	5✔
118	vector_rd.data_type != VectorDataType::Data,	5✔
119	error::InvalidType {	×
120	expected: format!(	×
121	"{}, {} or {}",	×
122	VectorDataType::MultiPoint,	×
123	VectorDataType::MultiLineString,	×
124	VectorDataType::MultiPolygon	×
125	),	×
126	found: VectorDataType::Data.to_string()	×
127	},	×
128	);
129
130	let raster_sources = futures::future::try_join_all(	5✔
131	self.sources	5✔
132	.rasters	5✔
133	.into_iter()	5✔
134	.enumerate()	5✔
135	.map(\|(i, op)\| op.initialize(path.clone_and_append(i as u8 + 1), context)),	5✔
136	)	5✔
137	.await?;	×
138
139	let spatial_reference = vector_rd.spatial_reference;	5✔
140
141	for other_spatial_reference in raster_sources	5✔
142	.iter()	5✔
143	.map(\|source\| source.result_descriptor().spatial_reference)	5✔
144	{
145	ensure!(	5✔
146	spatial_reference == other_spatial_reference,	5✔
147	crate::error::InvalidSpatialReference {	1✔
148	expected: spatial_reference,	1✔
149	found: other_spatial_reference,	1✔
150	}	1✔
151	);
152	}
153
154	let params = self.params;	4✔
155		4✔
156	let result_descriptor = vector_rd.map_columns(\|columns\| {	4✔
157	let mut columns = columns.clone();	4✔
158	for (i, new_column_name) in params.names.iter().enumerate() {	4✔
159	let feature_data_type = match params.temporal_aggregation {	4✔
160	TemporalAggregationMethod::First \| TemporalAggregationMethod::None => {
161	match raster_sources[i].result_descriptor().data_type {	2✔
162	RasterDataType::U8
163	\| RasterDataType::U16
164	\| RasterDataType::U32
165	\| RasterDataType::U64
166	\| RasterDataType::I8
167	\| RasterDataType::I16
168	\| RasterDataType::I32
169	\| RasterDataType::I64 => FeatureDataType::Int,	2✔
170	RasterDataType::F32 \| RasterDataType::F64 => FeatureDataType::Float,	×
171	}
172	}
173	TemporalAggregationMethod::Mean => FeatureDataType::Float,	2✔
174	};
175	columns.insert(	4✔
176	new_column_name.clone(),	4✔
177	VectorColumnInfo {	4✔
178	data_type: feature_data_type,	4✔
179	measurement: raster_sources[i].result_descriptor().measurement.clone(),	4✔
180	},	4✔
181	);	4✔
182	}
183	columns	4✔
184	});	4✔
185		4✔
186	Ok(InitializedRasterVectorJoin {	4✔
187	name,	4✔
188	result_descriptor,	4✔
189	vector_source,	4✔
190	raster_sources,	4✔
191	state: params,	4✔
192	}	4✔
193	.boxed())	4✔
194	}	10✔
195
196	span_fn!(RasterVectorJoin);	×
197	}
198
199	pub struct InitializedRasterVectorJoin {
200	name: CanonicOperatorName,
201	result_descriptor: VectorResultDescriptor,
202	vector_source: Box<dyn InitializedVectorOperator>,
203	raster_sources: Vec<Box<dyn InitializedRasterOperator>>,
204	state: RasterVectorJoinParams,
205	}
206
207	impl InitializedVectorOperator for InitializedRasterVectorJoin {
208	fn result_descriptor(&self) -> &VectorResultDescriptor {	1✔
209	&self.result_descriptor	1✔
210	}	1✔
211
212	fn query_processor(&self) -> Result<TypedVectorQueryProcessor> {	4✔
213	let typed_raster_processors = self	4✔
214	.raster_sources	4✔
215	.iter()	4✔
216	.map(InitializedRasterOperator::query_processor)	4✔
217	.collect::<Result<Vec<_>>>()?;	4✔
218
219	Ok(match self.vector_source.query_processor()? {	4✔
220	TypedVectorQueryProcessor::Data(_) => unreachable!(),	×
221	TypedVectorQueryProcessor::MultiPoint(points) => {	4✔
222	TypedVectorQueryProcessor::MultiPoint(match self.state.temporal_aggregation {	4✔
223	TemporalAggregationMethod::None => RasterVectorJoinProcessor::new(	×
224	points,	×
225	typed_raster_processors,	×
226	self.state.names.clone(),	×
227	self.state.feature_aggregation,	×
228	self.state.feature_aggregation_ignore_no_data,	×
229	)	×
230	.boxed(),	×
231	TemporalAggregationMethod::First \| TemporalAggregationMethod::Mean => {
232	RasterVectorAggregateJoinProcessor::new(	4✔
233	points,	4✔
234	typed_raster_processors,	4✔
235	self.state.names.clone(),	4✔
236	self.state.feature_aggregation,	4✔
237	self.state.feature_aggregation_ignore_no_data,	4✔
238	self.state.temporal_aggregation,	4✔
239	self.state.temporal_aggregation_ignore_no_data,	4✔
240	)	4✔
241	.boxed()	4✔
242	}
243	})
244	}
245	TypedVectorQueryProcessor::MultiPolygon(polygons) => {	×
246	TypedVectorQueryProcessor::MultiPolygon(match self.state.temporal_aggregation {	×
247	TemporalAggregationMethod::None => RasterVectorJoinProcessor::new(	×
248	polygons,	×
249	typed_raster_processors,	×
250	self.state.names.clone(),	×
251	self.state.feature_aggregation,	×
252	self.state.feature_aggregation_ignore_no_data,	×
253	)	×
254	.boxed(),	×
255	TemporalAggregationMethod::First \| TemporalAggregationMethod::Mean => {
256	RasterVectorAggregateJoinProcessor::new(	×
257	polygons,	×
258	typed_raster_processors,	×
259	self.state.names.clone(),	×
260	self.state.feature_aggregation,	×
261	self.state.feature_aggregation_ignore_no_data,	×
262	self.state.temporal_aggregation,	×
263	self.state.temporal_aggregation_ignore_no_data,	×
264	)	×
265	.boxed()	×
266	}
267	})
268	}
269	TypedVectorQueryProcessor::MultiLineString(_) => return Err(Error::NotYetImplemented),	×
270	})
271	}	4✔
272
273	fn canonic_name(&self) -> CanonicOperatorName {	×
274	self.name.clone()	×
275	}	×
276	}
277
278	pub fn create_feature_aggregator<P: Pixel>(	23✔
279	number_of_features: usize,	23✔
280	aggregation: FeatureAggregationMethod,	23✔
281	ignore_no_data: bool,	23✔
282	) -> TypedAggregator {	23✔
283	match aggregation {	23✔
284	FeatureAggregationMethod::First => match P::TYPE {	15✔
285	RasterDataType::U8
286	\| RasterDataType::U16
287	\| RasterDataType::U32
288	\| RasterDataType::U64
289	\| RasterDataType::I8
290	\| RasterDataType::I16
291	\| RasterDataType::I32
292	\| RasterDataType::I64 => {
293	FirstValueIntAggregator::new(number_of_features, ignore_no_data).into_typed()	15✔
294	}
295	RasterDataType::F32 \| RasterDataType::F64 => {
296	FirstValueFloatAggregator::new(number_of_features, ignore_no_data).into_typed()	×
297	}
298	},
299	FeatureAggregationMethod::Mean => {
300	MeanValueAggregator::new(number_of_features, ignore_no_data).into_typed()	8✔
301	}
302	}
303	}	23✔
304
305	#[cfg(test)]
306	mod tests {
307	use super::*;
308	use std::str::FromStr;
309
310	use crate::engine::{
311	ChunkByteSize, MockExecutionContext, MockQueryContext, QueryProcessor, RasterOperator,
312	};
313	use crate::mock::MockFeatureCollectionSource;
314	use crate::source::{GdalSource, GdalSourceParameters};
315	use crate::util::gdal::add_ndvi_dataset;
316	use futures::StreamExt;
317	use geoengine_datatypes::collections::{FeatureCollectionInfos, MultiPointCollection};
318	use geoengine_datatypes::dataset::NamedData;
319	use geoengine_datatypes::primitives::CacheHint;
320	use geoengine_datatypes::primitives::{
321	BoundingBox2D, DataRef, DateTime, FeatureDataRef, MultiPoint, SpatialResolution,
322	TimeInterval, VectorQueryRectangle,
323	};
324	use geoengine_datatypes::spatial_reference::SpatialReference;
325	use geoengine_datatypes::util::{gdal::hide_gdal_errors, test::TestDefault};
326	use serde_json::json;
327
328	#[test]	1✔
329	fn serialization() {	1✔
330	let raster_vector_join = RasterVectorJoin {	1✔
331	params: RasterVectorJoinParams {	1✔
332	names: ["foo", "bar"].iter().copied().map(str::to_string).collect(),	1✔
333	feature_aggregation: FeatureAggregationMethod::First,	1✔
334	feature_aggregation_ignore_no_data: false,	1✔
335	temporal_aggregation: TemporalAggregationMethod::Mean,	1✔
336	temporal_aggregation_ignore_no_data: false,	1✔
337	},	1✔
338	sources: SingleVectorMultipleRasterSources {	1✔
339	vector: MockFeatureCollectionSource::<MultiPoint>::multiple(vec![]).boxed(),	1✔
340	rasters: vec![],	1✔
341	},	1✔
342	};	1✔
343		1✔
344	let serialized = json!({	1✔
345	"type": "RasterVectorJoin",	1✔
346	"params": {	1✔
347	"names": ["foo", "bar"],	1✔
348	"featureAggregation": "first",	1✔
349	"temporalAggregation": "mean",	1✔
350	},	1✔
351	"sources": {	1✔
352	"vector": {	1✔
353	"type": "MockFeatureCollectionSourceMultiPoint",	1✔
354	"params": {	1✔
355	"collections": [],	1✔
356	"spatialReference": "EPSG:4326",	1✔
357	"measurements": {},	1✔
358	}	1✔
359	},	1✔
360	"rasters": [],	1✔
361	},	1✔
362	})	1✔
363	.to_string();	1✔
364		1✔
365	let deserialized: RasterVectorJoin = serde_json::from_str(&serialized).unwrap();	1✔
366		1✔
367	assert_eq!(deserialized.params, raster_vector_join.params);	1✔
368	}	1✔
369
370	fn ndvi_source(name: NamedData) -> Box<dyn RasterOperator> {	4✔
371	let gdal_source = GdalSource {	4✔
372	params: GdalSourceParameters { data: name },	4✔
373	};	4✔
374		4✔
375	gdal_source.boxed()	4✔
376	}	4✔
377
378	#[tokio::test]	1✔
379	async fn ndvi_time_point() {	1✔
380	let point_source = MockFeatureCollectionSource::single(	1✔
381	MultiPointCollection::from_data(	1✔
382	MultiPoint::many(vec![	1✔
383	(-13.95, 20.05),	1✔
384	(-14.05, 20.05),	1✔
385	(-13.95, 19.95),	1✔
386	(-14.05, 19.95),	1✔
387	])	1✔
388	.unwrap(),	1✔
389	vec![	1✔
390	TimeInterval::new(	1✔
391	DateTime::new_utc(2014, 1, 1, 0, 0, 0),	1✔
392	DateTime::new_utc(2014, 1, 1, 0, 0, 0),	1✔
393	)	1✔
394	.unwrap();	1✔
395	4	1✔
396	],	1✔
397	Default::default(),	1✔
398	CacheHint::default(),	1✔
399	)	1✔
400	.unwrap(),	1✔
401	)	1✔
402	.boxed();	1✔
403		1✔
404	let mut exe_ctc = MockExecutionContext::test_default();	1✔
405	let ndvi_id = add_ndvi_dataset(&mut exe_ctc);	1✔
406		1✔
407	let operator = RasterVectorJoin {	1✔
408	params: RasterVectorJoinParams {	1✔
409	names: vec!["ndvi".to_string()],	1✔
410	feature_aggregation: FeatureAggregationMethod::First,	1✔
411	feature_aggregation_ignore_no_data: false,	1✔
412	temporal_aggregation: TemporalAggregationMethod::First,	1✔
413	temporal_aggregation_ignore_no_data: false,	1✔
414	},	1✔
415	sources: SingleVectorMultipleRasterSources {	1✔
416	vector: point_source,	1✔
417	rasters: vec![ndvi_source(ndvi_id.clone())],	1✔
418	},	1✔
419	};	1✔
420
421	let operator = operator	1✔
422	.boxed()	1✔
423	.initialize(WorkflowOperatorPath::initialize_root(), &exe_ctc)	1✔
424	.await	×
425	.unwrap();	1✔
426		1✔
427	let query_processor = operator.query_processor().unwrap().multi_point().unwrap();	1✔
428
429	let result = query_processor	1✔
430	.query(	1✔
431	VectorQueryRectangle {	1✔
432	spatial_bounds: BoundingBox2D::new((-180., -90.).into(), (180., 90.).into())	1✔
433	.unwrap(),	1✔
434	time_interval: TimeInterval::default(),	1✔
435	spatial_resolution: SpatialResolution::new(0.1, 0.1).unwrap(),	1✔
436	},	1✔
437	&MockQueryContext::new(ChunkByteSize::MIN),	1✔
438	)	1✔
439	.await	×
440	.unwrap()	1✔
441	.map(Result::unwrap)	1✔
442	.collect::<Vec<MultiPointCollection>>()	1✔
443	.await;	33✔
444
445	assert_eq!(result.len(), 1);	1✔
446
447	let FeatureDataRef::Int(data) = result[0].data("ndvi").unwrap() else {	1✔
448	unreachable!();	×
449	};
450
451	// these values are taken from loading the tiff in QGIS
452	assert_eq!(data.as_ref(), &[54, 55, 51, 55]);	1✔
453	}
454
455	#[tokio::test]	1✔
456	#[allow(clippy::float_cmp)]
457	async fn ndvi_time_range() {	1✔
458	let point_source = MockFeatureCollectionSource::single(	1✔
459	MultiPointCollection::from_data(	1✔
460	MultiPoint::many(vec![	1✔
461	(-13.95, 20.05),	1✔
462	(-14.05, 20.05),	1✔
463	(-13.95, 19.95),	1✔
464	(-14.05, 19.95),	1✔
465	])	1✔
466	.unwrap(),	1✔
467	vec![	1✔
468	TimeInterval::new(	1✔
469	DateTime::new_utc(2014, 1, 1, 0, 0, 0),	1✔
470	DateTime::new_utc(2014, 3, 1, 0, 0, 0),	1✔
471	)	1✔
472	.unwrap();	1✔
473	4	1✔
474	],	1✔
475	Default::default(),	1✔
476	CacheHint::default(),	1✔
477	)	1✔
478	.unwrap(),	1✔
479	)	1✔
480	.boxed();	1✔
481		1✔
482	let mut exe_ctc = MockExecutionContext::test_default();	1✔
483	let ndvi_id = add_ndvi_dataset(&mut exe_ctc);	1✔
484		1✔
485	let operator = RasterVectorJoin {	1✔
486	params: RasterVectorJoinParams {	1✔
487	names: vec!["ndvi".to_string()],	1✔
488	feature_aggregation: FeatureAggregationMethod::First,	1✔
489	feature_aggregation_ignore_no_data: false,	1✔
490	temporal_aggregation: TemporalAggregationMethod::Mean,	1✔
491	temporal_aggregation_ignore_no_data: false,	1✔
492	},	1✔
493	sources: SingleVectorMultipleRasterSources {	1✔
494	vector: point_source,	1✔
495	rasters: vec![ndvi_source(ndvi_id.clone())],	1✔
496	},	1✔
497	};	1✔
498
499	let operator = operator	1✔
500	.boxed()	1✔
501	.initialize(WorkflowOperatorPath::initialize_root(), &exe_ctc)	1✔
502	.await	×
503	.unwrap();	1✔
504		1✔
505	let query_processor = operator.query_processor().unwrap().multi_point().unwrap();	1✔
506
507	let result = query_processor	1✔
508	.query(	1✔
509	VectorQueryRectangle {	1✔
510	spatial_bounds: BoundingBox2D::new((-180., -90.).into(), (180., 90.).into())	1✔
511	.unwrap(),	1✔
512	time_interval: TimeInterval::default(),	1✔
513	spatial_resolution: SpatialResolution::new(0.1, 0.1).unwrap(),	1✔
514	},	1✔
515	&MockQueryContext::new(ChunkByteSize::MIN),	1✔
516	)	1✔
517	.await	×
518	.unwrap()	1✔
519	.map(Result::unwrap)	1✔
520	.collect::<Vec<MultiPointCollection>>()	1✔
521	.await;	65✔
522
523	assert_eq!(result.len(), 1);	1✔
524
525	let FeatureDataRef::Float(data) = result[0].data("ndvi").unwrap() else {	1✔
526	unreachable!();	×
527	};
528
529	// these values are taken from loading the tiff in QGIS
530	assert_eq!(	1✔
531	data.as_ref(),	1✔
532	&[	1✔
533	(54. + 52.) / 2.,	1✔
534	(55. + 55.) / 2.,	1✔
535	(51. + 50.) / 2.,	1✔
536	(55. + 53.) / 2.,	1✔
537	]	1✔
538	);	1✔
539	}
540
541	#[tokio::test]	1✔
542	#[allow(clippy::float_cmp)]
543	async fn ndvi_with_default_time() {	1✔
544	hide_gdal_errors();	1✔
545		1✔
546	let point_source = MockFeatureCollectionSource::single(	1✔
547	MultiPointCollection::from_data(	1✔
548	MultiPoint::many(vec![	1✔
549	(-13.95, 20.05),	1✔
550	(-14.05, 20.05),	1✔
551	(-13.95, 19.95),	1✔
552	(-14.05, 19.95),	1✔
553	])	1✔
554	.unwrap(),	1✔
555	vec![TimeInterval::default(); 4],	1✔
556	Default::default(),	1✔
557	CacheHint::default(),	1✔
558	)	1✔
559	.unwrap(),	1✔
560	)	1✔
561	.boxed();	1✔
562		1✔
563	let mut exe_ctc = MockExecutionContext::test_default();	1✔
564	let ndvi_id = add_ndvi_dataset(&mut exe_ctc);	1✔
565		1✔
566	let operator = RasterVectorJoin {	1✔
567	params: RasterVectorJoinParams {	1✔
568	names: vec!["ndvi".to_string()],	1✔
569	feature_aggregation: FeatureAggregationMethod::First,	1✔
570	feature_aggregation_ignore_no_data: false,	1✔
571	temporal_aggregation: TemporalAggregationMethod::Mean,	1✔
572	temporal_aggregation_ignore_no_data: false,	1✔
573	},	1✔
574	sources: SingleVectorMultipleRasterSources {	1✔
575	vector: point_source,	1✔
576	rasters: vec![ndvi_source(ndvi_id.clone())],	1✔
577	},	1✔
578	};	1✔
579
580	let operator = operator	1✔
581	.boxed()	1✔
582	.initialize(WorkflowOperatorPath::initialize_root(), &exe_ctc)	1✔
583	.await	×
584	.unwrap();	1✔
585		1✔
586	let query_processor = operator.query_processor().unwrap().multi_point().unwrap();	1✔
587
588	let result = query_processor	1✔
589	.query(	1✔
590	VectorQueryRectangle {	1✔
591	spatial_bounds: BoundingBox2D::new((-180., -90.).into(), (180., 90.).into())	1✔
592	.unwrap(),	1✔
593	time_interval: TimeInterval::default(),	1✔
594	spatial_resolution: SpatialResolution::new(0.1, 0.1).unwrap(),	1✔
595	},	1✔
596	&MockQueryContext::new(ChunkByteSize::MIN),	1✔
597	)	1✔
598	.await	×
599	.unwrap()	1✔
600	.map(Result::unwrap)	1✔
601	.collect::<Vec<MultiPointCollection>>()	1✔
602	.await;	3✔
603
604	assert_eq!(result.len(), 1);	1✔
605
606	let FeatureDataRef::Float(data) = result[0].data("ndvi").unwrap() else {	1✔
607	unreachable!();	×
608	};
609
610	assert_eq!(data.as_ref(), &[0., 0., 0., 0.]);	1✔
611
612	assert_eq!(data.nulls(), vec![true, true, true, true]);	1✔
613	}
614
615	#[tokio::test]	1✔
616	async fn it_checks_sref() {	1✔
617	let point_source = MockFeatureCollectionSource::with_collections_and_sref(	1✔
618	vec![MultiPointCollection::from_data(	1✔
619	MultiPoint::many(vec![	1✔
620	(-13.95, 20.05),	1✔
621	(-14.05, 20.05),	1✔
622	(-13.95, 19.95),	1✔
623	(-14.05, 19.95),	1✔
624	])	1✔
625	.unwrap(),	1✔
626	vec![TimeInterval::default(); 4],	1✔
627	Default::default(),	1✔
628	CacheHint::default(),	1✔
629	)	1✔
630	.unwrap()],	1✔
631	SpatialReference::from_str("EPSG:3857").unwrap(),	1✔
632	)	1✔
633	.boxed();	1✔
634		1✔
635	let mut exe_ctc = MockExecutionContext::test_default();	1✔
636	let ndvi_id = add_ndvi_dataset(&mut exe_ctc);	1✔
637
638	let operator = RasterVectorJoin {	1✔
639	params: RasterVectorJoinParams {	1✔
640	names: vec!["ndvi".to_string()],	1✔
641	feature_aggregation: FeatureAggregationMethod::First,	1✔
642	feature_aggregation_ignore_no_data: false,	1✔
643	temporal_aggregation: TemporalAggregationMethod::Mean,	1✔
644	temporal_aggregation_ignore_no_data: false,	1✔
645	},	1✔
646	sources: SingleVectorMultipleRasterSources {	1✔
647	vector: point_source,	1✔
648	rasters: vec![ndvi_source(ndvi_id.clone())],	1✔
649	},	1✔
650	}	1✔
651	.boxed()	1✔
652	.initialize(WorkflowOperatorPath::initialize_root(), &exe_ctc)	1✔
653	.await;	×
654
655	assert!(matches!(	1✔
656	operator,	1✔
657	Err(Error::InvalidSpatialReference {
658	expected: _,
659	found: _,
660	})
661	));
662	}
663	}

geo-engine / geoengine / 6510717572

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