7006568925

Committed 27 Nov 2023 02:07PM UTC coverage: 89.651% (+0.2%) from 89.498%

Build # 7006568925

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #888 from geo-engine/raster_stacks

raster stacking

Run Details

4032 of 4274 new or added lines in 107 files covered. (94.34%)

12 existing lines in 8 files now uncovered.

113020 of 126066 relevant lines covered (89.65%)

59901.79 hits per line

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

96.4

/operators/src/processing/meteosat/reflectance.rs

use std::sync::Arc;

use crate::engine::{
    CanonicOperatorName, ExecutionContext, InitializedRasterOperator, InitializedSources, Operator,
    OperatorName, QueryContext, QueryProcessor, RasterBandDescriptor, RasterBandDescriptors,
    RasterOperator, RasterQueryProcessor, RasterResultDescriptor, SingleRasterSource,
    TypedRasterQueryProcessor, WorkflowOperatorPath,
};
use crate::util::Result;
use async_trait::async_trait;
use num_traits::AsPrimitive;
use rayon::ThreadPool;
use snafu::ensure;
use TypedRasterQueryProcessor::F32 as QueryProcessorOut;

use crate::error::Error;
use futures::stream::BoxStream;
use futures::{StreamExt, TryStreamExt};
use geoengine_datatypes::primitives::{
    BandSelection, ClassificationMeasurement, ContinuousMeasurement, DateTime, Measurement,
    RasterQueryRectangle, SpatialPartition2D,
};
use geoengine_datatypes::raster::{
    GridIdx2D, MapIndexedElementsParallel, RasterDataType, RasterPropertiesKey, RasterTile2D,
};
use serde::{Deserialize, Serialize};

// Output type is always f32
type PixelOut = f32;
use crate::processing::meteosat::satellite::{Channel, Satellite};
use crate::processing::meteosat::{new_channel_key, new_satellite_key};
use crate::util::sunpos::SunPos;
use RasterDataType::F32 as RasterOut;

/// Parameters for the `Reflectance` operator.
/// * `solar_correction` switch to enable solar correction.
/// * `force_hrv` switch to force the use of the hrv channel.
/// * `force_satellite` forces the use of the satellite with the given name.
#[derive(Debug, Serialize, Deserialize, PartialEq, Eq, Clone, Default, Copy)]
#[serde(rename_all = "camelCase")]
pub struct ReflectanceParams {
    pub solar_correction: bool,
    #[serde(rename = "forceHRV")]
    pub force_hrv: bool,
    pub force_satellite: Option<u8>,
}

/// The reflectance operator consumes an MSG image preprocessed
/// via the radiance operator and computes the reflectance value
/// from a given radiance raster.
pub type Reflectance = Operator<ReflectanceParams, SingleRasterSource>;

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

pub struct InitializedReflectance {
    name: CanonicOperatorName,
    result_descriptor: RasterResultDescriptor,
    source: Box<dyn InitializedRasterOperator>,
    params: ReflectanceParams,
}

#[typetag::serde]
#[async_trait]
impl RasterOperator for Reflectance {
    async fn _initialize(
        self: Box<Self>,
        path: WorkflowOperatorPath,
        context: &dyn ExecutionContext,
    ) -> Result<Box<dyn InitializedRasterOperator>> {
        let name = CanonicOperatorName::from(&self);

        let initialized_source = self.sources.initialize_sources(path, context).await?;
        let input = initialized_source.raster;

        let in_desc = input.result_descriptor();

        // TODO: implement multi-band functionality and remove this check
        ensure!(
            in_desc.bands.len() == 1,
            crate::error::OperatorDoesNotSupportMultiBandsSourcesYet {
                operator: Reflectance::TYPE_NAME
            }
        );

        match &in_desc.bands[0].measurement {
            Measurement::Continuous(ContinuousMeasurement {
                measurement: m,
                unit: _,
            }) if m != "radiance" => {
                return Err(Error::InvalidMeasurement {
                    expected: "radiance".into(),
                    found: m.clone(),
                })
            }
            Measurement::Classification(ClassificationMeasurement {
                measurement: m,
                classes: _,
            }) => {
                return Err(Error::InvalidMeasurement {
                    expected: "radiance".into(),
                    found: m.clone(),
                })
            }
            Measurement::Unitless => {
                return Err(Error::InvalidMeasurement {
                    expected: "radiance".into(),
                    found: "unitless".into(),
                })
            }
            // OK Case
            Measurement::Continuous(ContinuousMeasurement {
                measurement: _,
                unit: _,
            }) => {}
        }

        let out_desc = RasterResultDescriptor {
            spatial_reference: in_desc.spatial_reference,
            data_type: RasterOut,
            time: in_desc.time,
            bbox: in_desc.bbox,
            resolution: in_desc.resolution,
            bands: RasterBandDescriptors::new(vec![RasterBandDescriptor::new(
                in_desc.bands[0].name.clone(),
                Measurement::Continuous(ContinuousMeasurement {
                    measurement: "reflectance".into(),
                    unit: Some("fraction".into()),
                }),
            )])
            .unwrap(),
        };

        let initialized_operator = InitializedReflectance {
            name,
            result_descriptor: out_desc,
            source: input,
            params: self.params,
        };

        Ok(initialized_operator.boxed())
    }

    span_fn!(Reflectance);
}

impl InitializedRasterOperator for InitializedReflectance {
    fn result_descriptor(&self) -> &RasterResultDescriptor {
        &self.result_descriptor
    }

    fn query_processor(&self) -> Result<TypedRasterQueryProcessor, Error> {
        let q = self.source.query_processor()?;

        // We only support f32 input rasters
        if let TypedRasterQueryProcessor::F32(p) = q {
            Ok(QueryProcessorOut(Box::new(ReflectanceProcessor::new(
                p,
                self.params,
            ))))
        } else {
            Err(Error::InvalidRasterDataType)
        }
    }

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

struct ReflectanceProcessor<Q>
where
    Q: RasterQueryProcessor<RasterType = PixelOut>,
{
    source: Q,
    params: ReflectanceParams,
    channel_key: RasterPropertiesKey,
    satellite_key: RasterPropertiesKey,
}

impl<Q> ReflectanceProcessor<Q>
where
    Q: RasterQueryProcessor<RasterType = PixelOut>,
{
    pub fn new(source: Q, params: ReflectanceParams) -> Self {
        Self {
            source,
            params,
            channel_key: new_channel_key(),
            satellite_key: new_satellite_key(),
        }
    }

    fn channel<'a>(
        &self,
        tile: &RasterTile2D<PixelOut>,
        satellite: &'a Satellite,
    ) -> Result<&'a Channel> {
        if self.params.force_hrv {
            Ok(satellite.hrv())
        } else {
            let channel_id = tile
                .properties
                .number_property::<usize>(&self.channel_key)?
                - 1;
            satellite.channel(channel_id)
        }
    }

    fn satellite(&self, tile: &RasterTile2D<PixelOut>) -> Result<&'static Satellite> {
        let id = match self.params.force_satellite {
            Some(id) => id,
            _ => tile.properties.number_property(&self.satellite_key)?,
        };
        Satellite::satellite_by_msg_id(id)
    }

    async fn process_tile_async(
        &self,
        tile: RasterTile2D<PixelOut>,
        pool: Arc<ThreadPool>,
    ) -> Result<RasterTile2D<PixelOut>> {
        if tile.is_empty() {
            return Ok(tile);
        }

        let satellite = self.satellite(&tile)?;
        let channel = self.channel(&tile, satellite)?;
        let solar_correction = self.params.solar_correction;
        let timestamp = tile
            .time
            .start()
            .as_date_time()
            .ok_or(Error::InvalidUTCTimestamp)?;
        let sun_pos_option = if solar_correction {
            Some(SunPos::new(&timestamp))
        } else {
            None
        };
        let etsr = channel.etsr / std::f64::consts::PI;
        let esd = calculate_esd(&timestamp);
        let tile_geo_transform = tile.tile_geo_transform();

        let map_fn = move |grid_idx: GridIdx2D, pixel_option: Option<PixelOut>| {
            pixel_option.map(|p| {
                if let Some(sun_pos) = sun_pos_option {
                    let geos_coord =
                        tile_geo_transform.grid_idx_to_pixel_center_coordinate_2d(grid_idx);

                    let (lat, lon) = channel.view_angle_lat_lon(geos_coord, 0.0);
                    let (_, zenith) = sun_pos.solar_azimuth_zenith(lat, lon);

                    (f64::from(p) * esd * esd / (etsr * zenith.min(80.0).to_radians().cos()))
                        as PixelOut
                } else {
                    (f64::from(p) * esd * esd / etsr).as_()
                }
            })
        };

        let refl_tile = crate::util::spawn_blocking_with_thread_pool(pool, move || {
            tile.map_indexed_elements_parallel(map_fn)
        })
        .await?;

        Ok(refl_tile)
    }
}

fn calculate_esd(timestamp: &DateTime) -> f64 {
    let perihelion = f64::from(timestamp.day_of_year()) - 3.0;
    let e = 0.0167;
    let theta = std::f64::consts::TAU * (perihelion / 365.0);
    1.0 - e * theta.cos()
}

#[async_trait]
impl<Q> QueryProcessor for ReflectanceProcessor<Q>
where
    Q: QueryProcessor<
        Output = RasterTile2D<PixelOut>,
        SpatialBounds = SpatialPartition2D,
        Selection = BandSelection,
    >,
{
    type Output = RasterTile2D<PixelOut>;
    type SpatialBounds = SpatialPartition2D;
    type Selection = BandSelection;

    async fn _query<'a>(
        &'a self,
        query: RasterQueryRectangle,
        ctx: &'a dyn QueryContext,
    ) -> Result<BoxStream<'a, Result<Self::Output>>> {
        let src = self.source.query(query, ctx).await?;
        let rs = src.and_then(move |tile| self.process_tile_async(tile, ctx.thread_pool().clone()));
        Ok(rs.boxed())
    }
}

#[cfg(test)]
mod tests {
    use crate::engine::{MockExecutionContext, RasterOperator, SingleRasterSource};
    use crate::processing::meteosat::reflectance::{Reflectance, ReflectanceParams};
    use crate::processing::meteosat::test_util;
    use crate::util::Result;
    use geoengine_datatypes::primitives::{
        ClassificationMeasurement, ContinuousMeasurement, Measurement,
    };
    use geoengine_datatypes::raster::{
        EmptyGrid2D, Grid2D, GridOrEmpty, MaskedGrid2D, RasterTile2D, TilingSpecification,
    };
    use std::collections::HashMap;

    async fn process_mock(
        params: ReflectanceParams,
        channel: Option<u8>,
        satellite: Option<u8>,
        empty: bool,
        measurement: Option<Measurement>,
    ) -> Result<RasterTile2D<f32>> {
        let tile_size_in_pixels = [3, 2].into();
        let tiling_specification = TilingSpecification {
            origin_coordinate: [0.0, 0.0].into(),
            tile_size_in_pixels,
        };

        let ctx = MockExecutionContext::new_with_tiling_spec(tiling_specification);

        test_util::process(
            || {
                let props = test_util::create_properties(channel, satellite, None, None);
                let cc = if empty {
                    Some(EmptyGrid2D::new([3, 2].into()).into())
                } else {
                    None
                };

                let m = measurement.or_else(|| {
                    Some(Measurement::Continuous(ContinuousMeasurement {
                        measurement: "radiance".into(),
                        unit: Some("W·m^(-2)·sr^(-1)·cm^(-1)".into()),
                    }))
                });

                let src = test_util::create_mock_source::<f32>(props, cc, m);

                RasterOperator::boxed(Reflectance {
                    params,
                    sources: SingleRasterSource {
                        raster: src.boxed(),
                    },
                })
            },
            test_util::create_mock_query(),
            &ctx,
        )
        .await
    }

    // #[tokio::test]
    // async fn test_msg_raster() {
    //     let mut ctx = MockExecutionContext::test_default();
    //     let src = test_util::_create_gdal_src(&mut ctx);
    //
    //     let rad = Radiance {
    //         sources: SingleRasterSource {
    //             raster: src.boxed(),
    //         },
    //         params: RadianceParams {},
    //     };
    //
    //     let result = test_util::process(
    //         move || {
    //             RasterOperator::boxed(Reflectance {
    //                 params: ReflectanceParams::default(),
    //                 sources: SingleRasterSource {
    //                     raster: RasterOperator::boxed(rad),
    //                 },
    //             })
    //         },
    //         test_util::_create_gdal_query(),
    //         &ctx,
    //     )
    //     .await;
    //     assert!(result.as_ref().is_ok());
    // }

    #[tokio::test]
    async fn test_empty_ok() {
        let result = process_mock(ReflectanceParams::default(), Some(1), Some(1), true, None).await;

        assert!(result.is_ok());

        assert!(geoengine_datatypes::util::test::grid_or_empty_grid_eq(
            &result.unwrap().grid_array,
            &GridOrEmpty::from(EmptyGrid2D::new([3, 2].into()))
        ));
    }

    #[tokio::test]
    async fn test_ok_no_solar_correction() {
        let result =
            process_mock(ReflectanceParams::default(), Some(1), Some(1), false, None).await;

        assert!(result.is_ok());
        assert!(geoengine_datatypes::util::test::grid_or_empty_grid_eq(
            &result.as_ref().unwrap().grid_array,
            &MaskedGrid2D::new(
                Grid2D::new(
                    [3, 2].into(),
                    vec![
                        0.046_567_827_f32,
                        0.093_135_655_f32,
                        0.139_703_48_f32,
                        0.186_271_31_f32,
                        0.232_839_14_f32,
                        0.
                    ],
                )
                .unwrap(),
                Grid2D::new([3, 2].into(), vec![true, true, true, true, true, false,],).unwrap(),
            )
            .unwrap()
            .into()
        ));
    }

    #[tokio::test]
    async fn test_ok_force_satellite() {
        let params = ReflectanceParams {
            force_satellite: Some(4),
            ..Default::default()
        };
        let result = process_mock(params, Some(1), Some(1), false, None).await;

        assert!(result.is_ok());
        assert!(geoengine_datatypes::util::test::grid_or_empty_grid_eq(
            &result.as_ref().unwrap().grid_array,
            &MaskedGrid2D::new(
                Grid2D::new(
                    [3, 2].into(),
                    vec![
                        0.046_542_14_f32,
                        0.093_084_28_f32,
                        0.139_626_43_f32,
                        0.186_168_57_f32,
                        0.232_710_7_f32,
                        0. // TODO: check nodata mask
                    ],
                )
                .unwrap(),
                Grid2D::new([3, 2].into(), vec![true, true, true, true, true, false,],).unwrap()
            )
            .unwrap()
            .into()
        ));
    }

    #[tokio::test]
    async fn test_ok_force_hrv() {
        let params = ReflectanceParams {
            force_hrv: true,
            ..Default::default()
        };
        let result = process_mock(params, Some(1), Some(1), false, None).await;

        assert!(result.is_ok());
        assert!(geoengine_datatypes::util::test::grid_or_empty_grid_eq(
            &result.as_ref().unwrap().grid_array,
            &MaskedGrid2D::new(
                Grid2D::new(
                    [3, 2].into(),
                    vec![
                        0.038_567_86_f32,
                        0.077_135_72_f32,
                        0.115_703_575_f32,
                        0.154_271_44_f32,
                        0.192_839_3_f32,
                        0. // TODO: check nodata mask
                    ],
                )
                .unwrap(),
                Grid2D::new([3, 2].into(), vec![true, true, true, true, true, false,],).unwrap(),
            )
            .unwrap()
            .into()
        ));
    }

    #[tokio::test]
    async fn test_ok_solar_correction() {
        let params = ReflectanceParams {
            solar_correction: true,
            ..Default::default()
        };
        let result = process_mock(params, Some(1), Some(1), false, None).await;

        assert!(result.is_ok());
        assert!(geoengine_datatypes::util::test::grid_or_empty_grid_eq(
            &result.as_ref().unwrap().grid_array,
            &MaskedGrid2D::new(
                Grid2D::new(
                    [3, 2].into(),
                    vec![
                        0.268_173_43_f32,
                        0.536_346_85_f32,
                        0.804_520_3_f32,
                        1.072_693_7_f32,
                        1.340_867_2_f32,
                        0.
                    ],
                )
                .unwrap(),
                Grid2D::new([3, 2].into(), vec![true, true, true, true, true, false],).unwrap(),
            )
            .unwrap()
            .into()
        ));
    }

    #[tokio::test]
    async fn test_invalid_force_satellite() {
        let params = ReflectanceParams {
            force_satellite: Some(42),
            ..Default::default()
        };
        let result = process_mock(params, Some(1), Some(1), false, None).await;
        assert!(result.is_err());
    }

    #[tokio::test]
    async fn test_missing_satellite() {
        let params = ReflectanceParams::default();
        let result = process_mock(params, Some(1), None, false, None).await;
        assert!(result.is_err());
    }

    #[tokio::test]
    async fn test_invalid_satellite() {
        let params = ReflectanceParams::default();
        let result = process_mock(params, Some(42), None, false, None).await;
        assert!(result.is_err());
    }

    #[tokio::test]
    async fn test_missing_channel() {
        let params = ReflectanceParams::default();
        let result = process_mock(params, None, Some(1), false, None).await;
        assert!(result.is_err());
    }

    #[tokio::test]
    async fn test_invalid_channel() {
        let params = ReflectanceParams::default();
        let result = process_mock(params, Some(42), Some(1), false, None).await;
        assert!(result.is_err());
    }

    #[tokio::test]
    async fn test_invalid_measurement_unitless() {
        let params = ReflectanceParams::default();
        let result =
            process_mock(params, Some(1), Some(1), false, Some(Measurement::Unitless)).await;
        assert!(result.is_err());
    }

    #[tokio::test]
    async fn test_invalid_measurement_continuous() {
        let params = ReflectanceParams::default();
        let result = process_mock(
            params,
            Some(1),
            Some(1),
            false,
            Some(Measurement::Continuous(ContinuousMeasurement {
                measurement: "invalid".into(),
                unit: None,
            })),
        )
        .await;
        assert!(result.is_err());
    }

    #[tokio::test]
    async fn test_invalid_measurement_classification() {
        let params = ReflectanceParams::default();
        let result = process_mock(
            params,
            Some(1),
            Some(1),
            false,
            Some(Measurement::Classification(ClassificationMeasurement {
                measurement: "invalid".into(),
                classes: HashMap::new(),
            })),
        )
        .await;
        assert!(result.is_err());
    }
}

1	use std::sync::Arc;
2
3	use crate::engine::{
4	CanonicOperatorName, ExecutionContext, InitializedRasterOperator, InitializedSources, Operator,
5	OperatorName, QueryContext, QueryProcessor, RasterBandDescriptor, RasterBandDescriptors,
6	RasterOperator, RasterQueryProcessor, RasterResultDescriptor, SingleRasterSource,
7	TypedRasterQueryProcessor, WorkflowOperatorPath,
8	};
9	use crate::util::Result;
10	use async_trait::async_trait;
11	use num_traits::AsPrimitive;
12	use rayon::ThreadPool;
13	use snafu::ensure;
14	use TypedRasterQueryProcessor::F32 as QueryProcessorOut;
15
16	use crate::error::Error;
17	use futures::stream::BoxStream;
18	use futures::{StreamExt, TryStreamExt};
19	use geoengine_datatypes::primitives::{
20	BandSelection, ClassificationMeasurement, ContinuousMeasurement, DateTime, Measurement,
21	RasterQueryRectangle, SpatialPartition2D,
22	};
23	use geoengine_datatypes::raster::{
24	GridIdx2D, MapIndexedElementsParallel, RasterDataType, RasterPropertiesKey, RasterTile2D,
25	};
26	use serde::{Deserialize, Serialize};
27
28	// Output type is always f32
29	type PixelOut = f32;
30	use crate::processing::meteosat::satellite::{Channel, Satellite};
31	use crate::processing::meteosat::{new_channel_key, new_satellite_key};
32	use crate::util::sunpos::SunPos;
33	use RasterDataType::F32 as RasterOut;
34
35	/// Parameters for the `Reflectance` operator.
36	/// * `solar_correction` switch to enable solar correction.
37	/// * `force_hrv` switch to force the use of the hrv channel.
38	/// * `force_satellite` forces the use of the satellite with the given name.
39	#[derive(Debug, Serialize, Deserialize, PartialEq, Eq, Clone, Default, Copy)]	13✔
40	#[serde(rename_all = "camelCase")]
41	pub struct ReflectanceParams {
42	pub solar_correction: bool,
43	#[serde(rename = "forceHRV")]
44	pub force_hrv: bool,
45	pub force_satellite: Option<u8>,
46	}
47
48	/// The reflectance operator consumes an MSG image preprocessed
49	/// via the radiance operator and computes the reflectance value
50	/// from a given radiance raster.
51	pub type Reflectance = Operator<ReflectanceParams, SingleRasterSource>;
52
53	impl OperatorName for Reflectance {
54	const TYPE_NAME: &'static str = "Reflectance";
55	}
56
57	pub struct InitializedReflectance {
58	name: CanonicOperatorName,
59	result_descriptor: RasterResultDescriptor,
60	source: Box<dyn InitializedRasterOperator>,
61	params: ReflectanceParams,
62	}
63
64	#[typetag::serde]	×
65	#[async_trait]
66	impl RasterOperator for Reflectance {
67	async fn _initialize(	13✔
68	self: Box<Self>,	13✔
69	path: WorkflowOperatorPath,	13✔
70	context: &dyn ExecutionContext,	13✔
71	) -> Result<Box<dyn InitializedRasterOperator>> {	13✔
72	let name = CanonicOperatorName::from(&self);	13✔
73
74	let initialized_source = self.sources.initialize_sources(path, context).await?;	13✔
75	let input = initialized_source.raster;	13✔
76		13✔
77	let in_desc = input.result_descriptor();	13✔
78		13✔
79	// TODO: implement multi-band functionality and remove this check	13✔
80	ensure!(	13✔
81	in_desc.bands.len() == 1,	13✔
NEW 82	crate::error::OperatorDoesNotSupportMultiBandsSourcesYet {	×
NEW 83	operator: Reflectance::TYPE_NAME	×
NEW 84	}	×
85	);
86
87	match &in_desc.bands[0].measurement {	13✔
88	Measurement::Continuous(ContinuousMeasurement {
89	measurement: m,	11✔
90	unit: _,	11✔
91	}) if m != "radiance" => {	11✔
92	return Err(Error::InvalidMeasurement {	1✔
93	expected: "radiance".into(),	1✔
94	found: m.clone(),	1✔
95	})	1✔
96	}
97	Measurement::Classification(ClassificationMeasurement {
98	measurement: m,	1✔
99	classes: _,	1✔
100	}) => {	1✔
101	return Err(Error::InvalidMeasurement {	1✔
102	expected: "radiance".into(),	1✔
103	found: m.clone(),	1✔
104	})	1✔
105	}
106	Measurement::Unitless => {
107	return Err(Error::InvalidMeasurement {	1✔
108	expected: "radiance".into(),	1✔
109	found: "unitless".into(),	1✔
110	})	1✔
111	}
112	// OK Case
113	Measurement::Continuous(ContinuousMeasurement {
114	measurement: _,
115	unit: _,
116	}) => {}	10✔
117	}	10✔
118		10✔
119	let out_desc = RasterResultDescriptor {	10✔
120	spatial_reference: in_desc.spatial_reference,	10✔
121	data_type: RasterOut,	10✔
122	time: in_desc.time,	10✔
123	bbox: in_desc.bbox,	10✔
124	resolution: in_desc.resolution,	10✔
125	bands: RasterBandDescriptors::new(vec![RasterBandDescriptor::new(	10✔
126	in_desc.bands[0].name.clone(),	10✔
127	Measurement::Continuous(ContinuousMeasurement {	10✔
128	measurement: "reflectance".into(),	10✔
129	unit: Some("fraction".into()),	10✔
130	}),	10✔
131	)])	10✔
132	.unwrap(),	10✔
133	};	10✔
134		10✔
135	let initialized_operator = InitializedReflectance {	10✔
136	name,	10✔
137	result_descriptor: out_desc,	10✔
138	source: input,	10✔
139	params: self.params,	10✔
140	};	10✔
141		10✔
142	Ok(initialized_operator.boxed())	10✔
143	}	26✔
144
145	span_fn!(Reflectance);	×
146	}
147
148	impl InitializedRasterOperator for InitializedReflectance {
149	fn result_descriptor(&self) -> &RasterResultDescriptor {	×
150	&self.result_descriptor	×
151	}	×
152
153	fn query_processor(&self) -> Result<TypedRasterQueryProcessor, Error> {	10✔
154	let q = self.source.query_processor()?;	10✔
155
156	// We only support f32 input rasters
157	if let TypedRasterQueryProcessor::F32(p) = q {	10✔
158	Ok(QueryProcessorOut(Box::new(ReflectanceProcessor::new(	10✔
159	p,	10✔
160	self.params,	10✔
161	))))	10✔
162	} else {
163	Err(Error::InvalidRasterDataType)	×
164	}
165	}	10✔
166
167	fn canonic_name(&self) -> CanonicOperatorName {	×
168	self.name.clone()	×
169	}	×
170	}
171
172	struct ReflectanceProcessor<Q>
173	where
174	Q: RasterQueryProcessor<RasterType = PixelOut>,
175	{
176	source: Q,
177	params: ReflectanceParams,
178	channel_key: RasterPropertiesKey,
179	satellite_key: RasterPropertiesKey,
180	}
181
182	impl<Q> ReflectanceProcessor<Q>
183	where
184	Q: RasterQueryProcessor<RasterType = PixelOut>,
185	{
186	pub fn new(source: Q, params: ReflectanceParams) -> Self {	10✔
187	Self {	10✔
188	source,	10✔
189	params,	10✔
190	channel_key: new_channel_key(),	10✔
191	satellite_key: new_satellite_key(),	10✔
192	}	10✔
193	}	10✔
194
195	fn channel<'a>(	6✔
196	&self,	6✔
197	tile: &RasterTile2D<PixelOut>,	6✔
198	satellite: &'a Satellite,	6✔
199	) -> Result<&'a Channel> {	6✔
200	if self.params.force_hrv {	6✔
201	Ok(satellite.hrv())	1✔
202	} else {
203	let channel_id = tile	5✔
204	.properties	5✔
205	.number_property::<usize>(&self.channel_key)?	5✔
206	- 1;
207	satellite.channel(channel_id)	4✔
208	}
209	}	6✔
210
211	fn satellite(&self, tile: &RasterTile2D<PixelOut>) -> Result<&'static Satellite> {	9✔
212	let id = match self.params.force_satellite {	9✔
213	Some(id) => id,	2✔
214	_ => tile.properties.number_property(&self.satellite_key)?,	7✔
215	};
216	Satellite::satellite_by_msg_id(id)	7✔
217	}	9✔
218
219	async fn process_tile_async(	10✔
220	&self,	10✔
221	tile: RasterTile2D<PixelOut>,	10✔
222	pool: Arc<ThreadPool>,	10✔
223	) -> Result<RasterTile2D<PixelOut>> {	10✔
224	if tile.is_empty() {	10✔
225	return Ok(tile);	1✔
226	}	9✔
227
228	let satellite = self.satellite(&tile)?;	9✔
229	let channel = self.channel(&tile, satellite)?;	6✔
230	let solar_correction = self.params.solar_correction;	4✔
231	let timestamp = tile	4✔
232	.time	4✔
233	.start()	4✔
234	.as_date_time()	4✔
235	.ok_or(Error::InvalidUTCTimestamp)?;	4✔
236	let sun_pos_option = if solar_correction {	4✔
237	Some(SunPos::new(&timestamp))	1✔
238	} else {
239	None	3✔
240	};
241	let etsr = channel.etsr / std::f64::consts::PI;	4✔
242	let esd = calculate_esd(&timestamp);	4✔
243	let tile_geo_transform = tile.tile_geo_transform();	4✔
244		4✔
245	let map_fn = move \|grid_idx: GridIdx2D, pixel_option: Option<PixelOut>\| {	24✔
246	pixel_option.map(\|p\| {	24✔
247	if let Some(sun_pos) = sun_pos_option {	20✔
248	let geos_coord =	5✔
249	tile_geo_transform.grid_idx_to_pixel_center_coordinate_2d(grid_idx);	5✔
250		5✔
251	let (lat, lon) = channel.view_angle_lat_lon(geos_coord, 0.0);	5✔
252	let (_, zenith) = sun_pos.solar_azimuth_zenith(lat, lon);	5✔
253		5✔
254	(f64::from(p) * esd * esd / (etsr * zenith.min(80.0).to_radians().cos()))	5✔
255	as PixelOut	5✔
256	} else {
257	(f64::from(p) * esd * esd / etsr).as_()	15✔
258	}
259	})	24✔
260	};	24✔
261
262	let refl_tile = crate::util::spawn_blocking_with_thread_pool(pool, move \|\| {	4✔
263	tile.map_indexed_elements_parallel(map_fn)	4✔
264	})	4✔
265	.await?;	4✔
266
267	Ok(refl_tile)	4✔
268	}	10✔
269	}
270
271	fn calculate_esd(timestamp: &DateTime) -> f64 {	4✔
272	let perihelion = f64::from(timestamp.day_of_year()) - 3.0;	4✔
273	let e = 0.0167;	4✔
274	let theta = std::f64::consts::TAU * (perihelion / 365.0);	4✔
275	1.0 - e * theta.cos()	4✔
276	}	4✔
277
278	#[async_trait]
279	impl<Q> QueryProcessor for ReflectanceProcessor<Q>
280	where
281	Q: QueryProcessor<
282	Output = RasterTile2D<PixelOut>,
283	SpatialBounds = SpatialPartition2D,
284	Selection = BandSelection,
285	>,
286	{
287	type Output = RasterTile2D<PixelOut>;
288	type SpatialBounds = SpatialPartition2D;
289	type Selection = BandSelection;
290
291	async fn _query<'a>(	10✔
292	&'a self,	10✔
293	query: RasterQueryRectangle,	10✔
294	ctx: &'a dyn QueryContext,	10✔
295	) -> Result<BoxStream<'a, Result<Self::Output>>> {	10✔
296	let src = self.source.query(query, ctx).await?;	10✔
297	let rs = src.and_then(move \|tile\| self.process_tile_async(tile, ctx.thread_pool().clone()));	10✔
298	Ok(rs.boxed())	10✔
299	}	20✔
300	}
301
302	#[cfg(test)]
303	mod tests {
304	use crate::engine::{MockExecutionContext, RasterOperator, SingleRasterSource};
305	use crate::processing::meteosat::reflectance::{Reflectance, ReflectanceParams};
306	use crate::processing::meteosat::test_util;
307	use crate::util::Result;
308	use geoengine_datatypes::primitives::{
309	ClassificationMeasurement, ContinuousMeasurement, Measurement,
310	};
311	use geoengine_datatypes::raster::{
312	EmptyGrid2D, Grid2D, GridOrEmpty, MaskedGrid2D, RasterTile2D, TilingSpecification,
313	};
314	use std::collections::HashMap;
315
316	async fn process_mock(	13✔
317	params: ReflectanceParams,	13✔
318	channel: Option<u8>,	13✔
319	satellite: Option<u8>,	13✔
320	empty: bool,	13✔
321	measurement: Option<Measurement>,	13✔
322	) -> Result<RasterTile2D<f32>> {	13✔
323	let tile_size_in_pixels = [3, 2].into();	13✔
324	let tiling_specification = TilingSpecification {	13✔
325	origin_coordinate: [0.0, 0.0].into(),	13✔
326	tile_size_in_pixels,	13✔
327	};	13✔
328		13✔
329	let ctx = MockExecutionContext::new_with_tiling_spec(tiling_specification);	13✔
330		13✔
331	test_util::process(	13✔
332	\|\| {	13✔
333	let props = test_util::create_properties(channel, satellite, None, None);	13✔
334	let cc = if empty {	13✔
335	Some(EmptyGrid2D::new([3, 2].into()).into())	1✔
336	} else {
337	None	12✔
338	};
339
340	let m = measurement.or_else(\|\| {	13✔
341	Some(Measurement::Continuous(ContinuousMeasurement {	10✔
342	measurement: "radiance".into(),	10✔
343	unit: Some("W·m^(-2)·sr^(-1)·cm^(-1)".into()),	10✔
344	}))	10✔
345	});	13✔
346		13✔
347	let src = test_util::create_mock_source::<f32>(props, cc, m);	13✔
348		13✔
349	RasterOperator::boxed(Reflectance {	13✔
350	params,	13✔
351	sources: SingleRasterSource {	13✔
352	raster: src.boxed(),	13✔
353	},	13✔
354	})	13✔
355	},	13✔
356	test_util::create_mock_query(),	13✔
357	&ctx,	13✔
358	)	13✔
359	.await	4✔
360	}	13✔
361
362	// #[tokio::test]
363	// async fn test_msg_raster() {
364	// let mut ctx = MockExecutionContext::test_default();
365	// let src = test_util::_create_gdal_src(&mut ctx);
366	//
367	// let rad = Radiance {
368	// sources: SingleRasterSource {
369	// raster: src.boxed(),
370	// },
371	// params: RadianceParams {},
372	// };
373	//
374	// let result = test_util::process(
375	// move \|\| {
376	// RasterOperator::boxed(Reflectance {
377	// params: ReflectanceParams::default(),
378	// sources: SingleRasterSource {
379	// raster: RasterOperator::boxed(rad),
380	// },
381	// })
382	// },
383	// test_util::_create_gdal_query(),
384	// &ctx,
385	// )
386	// .await;
387	// assert!(result.as_ref().is_ok());
388	// }
389
390	#[tokio::test]	1✔
391	async fn test_empty_ok() {	1✔
392	let result = process_mock(ReflectanceParams::default(), Some(1), Some(1), true, None).await;	1✔
393
394	assert!(result.is_ok());	1✔
395
396	assert!(geoengine_datatypes::util::test::grid_or_empty_grid_eq(	1✔
397	&result.unwrap().grid_array,	1✔
398	&GridOrEmpty::from(EmptyGrid2D::new([3, 2].into()))	1✔
399	));	1✔
400	}
401
402	#[tokio::test]	1✔
403	async fn test_ok_no_solar_correction() {	1✔
404	let result =	1✔
405	process_mock(ReflectanceParams::default(), Some(1), Some(1), false, None).await;	1✔
406
407	assert!(result.is_ok());	1✔
408	assert!(geoengine_datatypes::util::test::grid_or_empty_grid_eq(	1✔
409	&result.as_ref().unwrap().grid_array,	1✔
410	&MaskedGrid2D::new(	1✔
411	Grid2D::new(	1✔
412	[3, 2].into(),	1✔
413	vec![	1✔
414	0.046_567_827_f32,	1✔
415	0.093_135_655_f32,	1✔
416	0.139_703_48_f32,	1✔
417	0.186_271_31_f32,	1✔
418	0.232_839_14_f32,	1✔
419	0.	1✔
420	],	1✔
421	)	1✔
422	.unwrap(),	1✔
423	Grid2D::new([3, 2].into(), vec![true, true, true, true, true, false,],).unwrap(),	1✔
424	)	1✔
425	.unwrap()	1✔
426	.into()	1✔
427	));	1✔
428	}
429
430	#[tokio::test]	1✔
431	async fn test_ok_force_satellite() {	1✔
432	let params = ReflectanceParams {	1✔
433	force_satellite: Some(4),	1✔
434	..Default::default()	1✔
435	};	1✔
436	let result = process_mock(params, Some(1), Some(1), false, None).await;	1✔
437
438	assert!(result.is_ok());	1✔
439	assert!(geoengine_datatypes::util::test::grid_or_empty_grid_eq(	1✔
440	&result.as_ref().unwrap().grid_array,	1✔
441	&MaskedGrid2D::new(	1✔
442	Grid2D::new(	1✔
443	[3, 2].into(),	1✔
444	vec![	1✔
445	0.046_542_14_f32,	1✔
446	0.093_084_28_f32,	1✔
447	0.139_626_43_f32,	1✔
448	0.186_168_57_f32,	1✔
449	0.232_710_7_f32,	1✔
450	0. // TODO: check nodata mask	1✔
451	],	1✔
452	)	1✔
453	.unwrap(),	1✔
454	Grid2D::new([3, 2].into(), vec![true, true, true, true, true, false,],).unwrap()	1✔
455	)	1✔
456	.unwrap()	1✔
457	.into()	1✔
458	));	1✔
459	}
460
461	#[tokio::test]	1✔
462	async fn test_ok_force_hrv() {	1✔
463	let params = ReflectanceParams {	1✔
464	force_hrv: true,	1✔
465	..Default::default()	1✔
466	};	1✔
467	let result = process_mock(params, Some(1), Some(1), false, None).await;	1✔
468
469	assert!(result.is_ok());	1✔
470	assert!(geoengine_datatypes::util::test::grid_or_empty_grid_eq(	1✔
471	&result.as_ref().unwrap().grid_array,	1✔
472	&MaskedGrid2D::new(	1✔
473	Grid2D::new(	1✔
474	[3, 2].into(),	1✔
475	vec![	1✔
476	0.038_567_86_f32,	1✔
477	0.077_135_72_f32,	1✔
478	0.115_703_575_f32,	1✔
479	0.154_271_44_f32,	1✔
480	0.192_839_3_f32,	1✔
481	0. // TODO: check nodata mask	1✔
482	],	1✔
483	)	1✔
484	.unwrap(),	1✔
485	Grid2D::new([3, 2].into(), vec![true, true, true, true, true, false,],).unwrap(),	1✔
486	)	1✔
487	.unwrap()	1✔
488	.into()	1✔
489	));	1✔
490	}
491
492	#[tokio::test]	1✔
493	async fn test_ok_solar_correction() {	1✔
494	let params = ReflectanceParams {	1✔
495	solar_correction: true,	1✔
496	..Default::default()	1✔
497	};	1✔
498	let result = process_mock(params, Some(1), Some(1), false, None).await;	1✔
499
500	assert!(result.is_ok());	1✔
501	assert!(geoengine_datatypes::util::test::grid_or_empty_grid_eq(	1✔
502	&result.as_ref().unwrap().grid_array,	1✔
503	&MaskedGrid2D::new(	1✔
504	Grid2D::new(	1✔
505	[3, 2].into(),	1✔
506	vec![	1✔
507	0.268_173_43_f32,	1✔
508	0.536_346_85_f32,	1✔
509	0.804_520_3_f32,	1✔
510	1.072_693_7_f32,	1✔
511	1.340_867_2_f32,	1✔
512	0.	1✔
513	],	1✔
514	)	1✔
515	.unwrap(),	1✔
516	Grid2D::new([3, 2].into(), vec![true, true, true, true, true, false],).unwrap(),	1✔
517	)	1✔
518	.unwrap()	1✔
519	.into()	1✔
520	));	1✔
521	}
522
523	#[tokio::test]	1✔
524	async fn test_invalid_force_satellite() {	1✔
525	let params = ReflectanceParams {	1✔
526	force_satellite: Some(42),	1✔
527	..Default::default()	1✔
528	};	1✔
529	let result = process_mock(params, Some(1), Some(1), false, None).await;	1✔
530	assert!(result.is_err());	1✔
531	}
532
533	#[tokio::test]	1✔
534	async fn test_missing_satellite() {	1✔
535	let params = ReflectanceParams::default();	1✔
536	let result = process_mock(params, Some(1), None, false, None).await;	1✔
537	assert!(result.is_err());	1✔
538	}
539
540	#[tokio::test]	1✔
541	async fn test_invalid_satellite() {	1✔
542	let params = ReflectanceParams::default();	1✔
543	let result = process_mock(params, Some(42), None, false, None).await;	1✔
544	assert!(result.is_err());	1✔
545	}
546
547	#[tokio::test]	1✔
548	async fn test_missing_channel() {	1✔
549	let params = ReflectanceParams::default();	1✔
550	let result = process_mock(params, None, Some(1), false, None).await;	1✔
551	assert!(result.is_err());	1✔
552	}
553
554	#[tokio::test]	1✔
555	async fn test_invalid_channel() {	1✔
556	let params = ReflectanceParams::default();	1✔
557	let result = process_mock(params, Some(42), Some(1), false, None).await;	1✔
558	assert!(result.is_err());	1✔
559	}
560
561	#[tokio::test]	1✔
562	async fn test_invalid_measurement_unitless() {	1✔
563	let params = ReflectanceParams::default();	1✔
564	let result =	1✔
565	process_mock(params, Some(1), Some(1), false, Some(Measurement::Unitless)).await;	1✔
566	assert!(result.is_err());	1✔
567	}
568
569	#[tokio::test]	1✔
570	async fn test_invalid_measurement_continuous() {	1✔
571	let params = ReflectanceParams::default();	1✔
572	let result = process_mock(	1✔
573	params,	1✔
574	Some(1),	1✔
575	Some(1),	1✔
576	false,	1✔
577	Some(Measurement::Continuous(ContinuousMeasurement {	1✔
578	measurement: "invalid".into(),	1✔
579	unit: None,	1✔
580	})),	1✔
581	)	1✔
582	.await;	×
583	assert!(result.is_err());	1✔
584	}
585
586	#[tokio::test]	1✔
587	async fn test_invalid_measurement_classification() {	1✔
588	let params = ReflectanceParams::default();	1✔
589	let result = process_mock(	1✔
590	params,	1✔
591	Some(1),	1✔
592	Some(1),	1✔
593	false,	1✔
594	Some(Measurement::Classification(ClassificationMeasurement {	1✔
595	measurement: "invalid".into(),	1✔
596	classes: HashMap::new(),	1✔
597	})),	1✔
598	)	1✔
599	.await;	×
600	assert!(result.is_err());	1✔
601	}
602	}

geo-engine / geoengine / 7006568925

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