12466060820

Committed 23 Dec 2024 11:26AM UTC coverage: 90.353% (-0.2%) from 90.512%

Build # 12466060820

Build Type

Pull #998

github

Committed by

web-flow

Commit Message

Merge 66ab0655c into 34e12969f

Pull Request Pull Request #998: Quota and Data usage Logging

Run Details

834 of 1211 new or added lines in 66 files covered. (68.87%)

222 existing lines in 18 files now uncovered.

133834 of 148123 relevant lines covered (90.35%)

54353.18 hits per line

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

96.01

/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 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,
    path: WorkflowOperatorPath,
    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.clone(), context)
            .await?;
        let input = initialized_source.raster;

        let in_desc = input.result_descriptor();

        for band in in_desc.bands.iter() {
            match &band.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(
                in_desc
                    .bands
                    .iter()
                    .map(|b| RasterBandDescriptor {
                        name: b.name.clone(),
                        measurement: Measurement::Continuous(ContinuousMeasurement {
                            measurement: "reflectance".into(),
                            unit: Some("fraction".into()),
                        }),
                    })
                    .collect::<Vec<_>>(),
            )?,
        };

        let initialized_operator = InitializedReflectance {
            name,
            path,
            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.result_descriptor.clone(),
                self.params,
            ))))
        } else {
            Err(Error::InvalidRasterDataType)
        }
    }

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

    fn name(&self) -> &'static str {
        Reflectance::TYPE_NAME
    }

    fn path(&self) -> WorkflowOperatorPath {
        self.path.clone()
    }
}

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

impl<Q> ReflectanceProcessor<Q>
where
    Q: RasterQueryProcessor<RasterType = PixelOut>,
{
    pub fn new(
        source: Q,
        result_descriptor: RasterResultDescriptor,
        params: ReflectanceParams,
    ) -> Self {
        Self {
            source,
            result_descriptor,
            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,
        ResultDescription = RasterResultDescriptor,
    >,
{
    type Output = RasterTile2D<PixelOut>;
    type SpatialBounds = SpatialPartition2D;
    type Selection = BandSelection;
    type ResultDescription = RasterResultDescriptor;

    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())
    }

    fn result_descriptor(&self) -> &Self::ResultDescription {
        &self.result_descriptor
    }
}

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

geo-engine / geoengine / 12466060820

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