Fixed tests

pyobs_cloudcover/pipeline/night/altaz_map_generator/altaz_map_generator.py

@@ -0,0 +1,37 @@
+import functools
+from typing import List, Optional, Tuple
+
+import numpy as np
+import numpy.typing as npt
+from cloudmap_rs import AltAzCoord
+
+from pyobs_cloudcover.world_model import WorldModel
+
+
+class AltAzMapGenerator:
+    def __init__(self, model: WorldModel, limiting_altitude: float):
+        self._model = model
+        self._limiting_altitude = limiting_altitude
+
+    @functools.lru_cache(maxsize=None)
+    def __call__(self, image_height: int, image_width: int) -> List[List[Optional[AltAzCoord]]]:
+        x = np.arange(0, image_width)
+        y = np.arange(0, image_height)
+
+        meshed_x, meshed_y = np.meshgrid(x, y)
+
+        alt, az = self._model.pix_to_altaz(meshed_x, meshed_y)
+        alt_az_coords = np.dstack((alt, az))
+
+        return [
+                [
+                    self._alt_az_convert(coord)
+                    for coord in row
+                ] for row in alt_az_coords
+            ]
+
+    def _alt_az_convert(self, alt_az_coords: Tuple[float, float]) -> Optional[AltAzCoord]:
+        if alt_az_coords[0] < np.deg2rad(self._limiting_altitude):
+            return None
+        else:
+            return AltAzCoord(*alt_az_coords)

pyobs_cloudcover/pipeline/night/altaz_map_generator/altaz_map_generator_factory.py

@@ -0,0 +1,12 @@
+from pyobs_cloudcover.pipeline.night.altaz_map_generator.altaz_map_generator import AltAzMapGenerator
+from pyobs_cloudcover.pipeline.night.catalog.catalog_constructor_options import CatalogConstructorOptions
+from pyobs_cloudcover.world_model import WorldModel
+
+
+class AltAzMapGeneratorFactory(object):
+    def __init__(self, options: CatalogConstructorOptions, model: WorldModel):
+        self._limiting_altitude = options.alt_filter
+        self._model = model
+
+    def __call__(self) -> AltAzMapGenerator:
+        return AltAzMapGenerator(self._model, self._limiting_altitude)

pyobs_cloudcover/pipeline/night/catalog/altaz_catalog.py

@@ -1,6 +1,8 @@
 import numpy as np
 import numpy.typing as npt
 
+from pyobs_cloudcover.world_model import WorldModel
+
 
 class AltAzCatalog(object):
     def __init__(self, sao: npt.NDArray[np.int_], alt: npt.NDArray[np.float_], az: npt.NDArray[np.float_],

pyobs_cloudcover/pipeline/night/catalog/pixel_catalog.py

@@ -10,15 +10,19 @@
 
 
 class PixelCatalog(object):
-    def __init__(self, sao: npt.NDArray[np.int_], px: npt.NDArray[np.float_], py: npt.NDArray[np.float_],
+    def __init__(self, sao: npt.NDArray[np.int_], alt: npt.NDArray[np.float_], az: npt.NDArray[np.float_], px: npt.NDArray[np.float_], py: npt.NDArray[np.float_],
                  v_mag: npt.NDArray[np.float_]) -> None:
         self.sao = sao
+        self.alt = alt
+        self.az = az
         self.px = px
         self.py = py
         self.v_mag = v_mag
 
     def _filter(self, condition: npt.NDArray[np.bool_]) -> None:
         self.sao = self.sao[condition]
+        self.alt = self.alt[condition]
+        self.az = self.az[condition]
         self.px = self.px[condition]
         self.py = self.py[condition]
         self.v_mag = self.v_mag[condition]
@@ -48,6 +52,7 @@ def from_altaz_catalog(cls, altaz_catalog: AltAzCatalog, model: WorldModel) -> P
 
         return PixelCatalog(
             altaz_catalog.sao,
+            altaz_catalog.alt, altaz_catalog.az,
             cast(npt.NDArray[np.float_], px), cast(npt.NDArray[np.float_], py),
             altaz_catalog.v_mag
             )

pyobs_cloudcover/pipeline/night/cloud_map_generator/cloud_map_generator.py

@@ -1,8 +1,8 @@
-from typing import List
+from typing import List, Optional
 
 import numpy as np
 import numpy.typing as npt
-from cloudmap_rs import Entry, gen_cloud_map, Average
+from cloudmap_rs import Entry, Average, MagnitudeMapGenerator, AltAzCoord
 
 from pyobs_cloudcover.pipeline.night.catalog.pixel_catalog import PixelCatalog
 
@@ -11,16 +11,24 @@ class CloudMapGenerator(object):
     def __init__(self, radius: float):
         self._radius = radius
 
-    def __call__(self, catalog: PixelCatalog, matches: List[bool], height: int, width: int) -> npt.NDArray[np.float_]:
-        std_entries = [Entry(*entry, True) for entry in zip(catalog.sao, catalog.px, catalog.py, catalog.v_mag)]
-        av_std_vis_map = gen_cloud_map(std_entries, self._radius, width, height)
+    def __call__(self, catalog: PixelCatalog, matches: List[bool], alt_az_image_list: List[List[Optional[AltAzCoord]]]) -> npt.NDArray[np.float_]:
+        star_coords = [AltAzCoord(np.deg2rad(alt), np.deg2rad(az)) for alt, az in zip(catalog.alt, catalog.az)]
+
+        map_generator = MagnitudeMapGenerator(star_coords, alt_az_image_list, np.deg2rad(self._radius))
+
+        std_entries = [Entry(v_mag, True) for v_mag in catalog.v_mag]
+        av_std_vis_map = map_generator.gen_cloud_map(std_entries)
         std_vis_map = self._convert_average_to_value(av_std_vis_map)
 
-        match_entries = [Entry(*entry) for entry in zip(catalog.sao, catalog.px, catalog.py, catalog.v_mag, matches)]
-        av_vis_map = gen_cloud_map(match_entries, self._radius, width, height)
+        no_entries = [Entry(v_mag, False) for v_mag in catalog.v_mag]
+        av_no_vis_map = map_generator.gen_cloud_map(no_entries)
+        no_vis_map = self._convert_average_to_value(av_no_vis_map)
+
+        match_entries = [Entry(v_mag, found) for v_mag, found in zip(catalog.v_mag, matches)]
+        av_vis_map = map_generator.gen_cloud_map(match_entries)
         vis_map = self._convert_average_to_value(av_vis_map)
 
-        return std_vis_map - vis_map
+        return (std_vis_map - vis_map) / (std_vis_map - no_vis_map)
 
     @staticmethod
     def _convert_average_to_value(average_map: List[List[Average]]) -> npt.NDArray[np.float_]:

pyobs_cloudcover/pipeline/night/pipeline.py

@@ -4,6 +4,7 @@
 import numpy.typing as npt
 
 from pyobs_cloudcover.cloud_coverage_info import CloudCoverageInfo
+from pyobs_cloudcover.pipeline.night.altaz_map_generator.altaz_map_generator import AltAzMapGenerator
 from pyobs_cloudcover.pipeline.night.catalog.catalog_constructor import CatalogConstructor
 from pyobs_cloudcover.pipeline.night.cloud_coverage_calculator.coverage_info_calculator import CoverageInfoCalculator
 from pyobs_cloudcover.pipeline.night.cloud_map_generator.cloud_map_generator import CloudMapGenerator
@@ -16,13 +17,15 @@ class NightPipeline(Pipeline):
     def __init__(self,
                  preprocess: Preprocessor,
                  catalog_constructor: CatalogConstructor,
+                 alt_az_list_generator: AltAzMapGenerator,
                  star_reverse_matcher: StarReverseMatcher,
                  cloud_map_generator: CloudMapGenerator,
                  coverage_info_calculator: CoverageInfoCalculator) -> None:
 
         self._preprocess = preprocess
         self._catalog_constructor = catalog_constructor
         self._star_reverse_matcher = star_reverse_matcher
+        self._alt_az_list_generator = alt_az_list_generator
         self._cloud_map_generator = cloud_map_generator
         self._coverage_info_calculator = coverage_info_calculator
 
@@ -33,5 +36,7 @@ def __call__(self, image: npt.NDArray[np.float_], obs_time: datetime.datetime) -
         catalog = self._catalog_constructor(obs_time, img_height, img_width)
         matches = self._star_reverse_matcher(preprocessed_image, catalog)
 
-        cloud_map = self._cloud_map_generator(catalog, matches, img_height, img_width)
+        alt_az_list = self._alt_az_list_generator(img_height, img_width)
+
+        cloud_map = self._cloud_map_generator(catalog, matches, alt_az_list)
         return self._coverage_info_calculator(cloud_map, obs_time)

pyobs_cloudcover/pipeline/night/pipeline_factory.py

@@ -1,5 +1,6 @@
 from astroplan import Observer
 
+from pyobs_cloudcover.pipeline.night.altaz_map_generator.altaz_map_generator_factory import AltAzMapGeneratorFactory
 from pyobs_cloudcover.pipeline.night.catalog.catalog_constructor_factory import CatalogConstructorFactory
 from pyobs_cloudcover.pipeline.night.cloud_coverage_calculator.cloud_info_calculator_factory import \
     CloudInfoCalculatorFactory
@@ -19,19 +20,22 @@ def __init__(self, observer: Observer, model: WorldModel):
     def __call__(self, options: NightPipelineOptions) -> NightPipeline:
         preprocessor_factory = PreprocessorFactory(options.preprocessor_options)
         catalog_constructor_factory = CatalogConstructorFactory(options.catalog_options, self._model, self._observer)
+        altaz_map_generator_factory = AltAzMapGeneratorFactory(options.catalog_options, self._model)
         reverse_matcher_factory = StarReverseMatcherFactory(options.star_matcher_options)
         cloud_map_generator_factory = CloudMapGeneratorFactory(options.cloud_generator_options)
         coverage_info_calculator_factory = CloudInfoCalculatorFactory(options.coverage_info_options, self._model)
 
         preprocessor = preprocessor_factory()
         catalog_constructor = catalog_constructor_factory()
+        altaz_map_generator = altaz_map_generator_factory()
         star_reverse_matcher = reverse_matcher_factory()
         cloud_map_generator = cloud_map_generator_factory()
         coverage_info_calculator = coverage_info_calculator_factory()
 
         pipeline = NightPipeline(
             preprocessor,
             catalog_constructor,
+            altaz_map_generator,
             star_reverse_matcher,
             cloud_map_generator,
             coverage_info_calculator

pyobs_cloudcover/world_model/hour_angle_converter.py

@@ -0,0 +1,15 @@
+import datetime
+
+from astroplan import Observer
+from astropy.coordinates import Angle
+
+
+class HourAngleConverter(object):
+    def __init__(self, observer: Observer):
+        self._observer = observer
+
+    def __call__(self, ra: Angle, obs_time: datetime.datetime) -> float:
+        lst = self._observer.local_sidereal_time(obs_time)
+        self._lst = Angle(lst, unit="hourangle")
+
+        return float((self._lst - ra).rad)

pyobs_cloudcover/world_model/wcs_model.py

@@ -0,0 +1,45 @@
+import datetime
+from typing import Union, Tuple
+
+import numpy as np
+from astroplan import Observer
+from astropy.coordinates import Angle, SkyCoord
+from astropy.wcs import WCS
+from numpy import typing as npt
+
+from pyobs_cloudcover.world_model import WorldModel
+from pyobs_cloudcover.world_model.hour_angle_converter import HourAngleConverter
+
+from astropy import units as u
+
+
+class WCSModel(WorldModel):
+    def __init__(self, wcs: WCS, observer: Observer, config_obs_time: datetime.datetime) -> None:
+        self._wcs = wcs
+        self._observer = observer
+        self._config_obs_time = config_obs_time
+
+        self._ha_conv = HourAngleConverter(observer)
+
+    def pix_to_altaz(self, x: Union[npt.NDArray[np.float_], float], y: Union[npt.NDArray[np.float_], float]) -> \
+            Union[Tuple[npt.NDArray[np.float_], npt.NDArray[np.float_]], Tuple[float, float]]:
+
+        coord = self._wcs.pixel_to_world(x, y)
+        ha = self._ha_conv(Angle(coord.ra.deg, unit="deg"), self._config_obs_time)
+
+        coord = SkyCoord(ra=ha * u.rad, dec=coord.dec, location=self._observer.location, obstime=self._config_obs_time)
+        coord = coord.transform_to("altaz")
+        return coord.alt.rad, coord.az.rad
+
+    def altaz_to_pix(self, alt: Union[npt.NDArray[np.float_], float], az: Union[npt.NDArray[np.float_], float]) -> \
+                Union[Tuple[npt.NDArray[np.float_], npt.NDArray[np.float_]], Tuple[float, float]]:
+
+        coord = SkyCoord(alt=alt * u.rad, az=az * u.rad, frame="altaz",
+                         location=self._observer.location, obstime=self._config_obs_time)
+        coord = coord.transform_to(frame="icrs")
+        ha = self._ha_conv(Angle(coord.ra.deg, unit="deg"), self._config_obs_time)
+
+        coord = SkyCoord(ra=ha * u.rad, dec=coord.dec, location=self._observer.location, obstime=self._config_obs_time)
+
+        pixel = self._wcs.world_to_pixel(coord)
+        return pixel[0], pixel[1]

pyobs_cloudcover/world_model/wcs_model_loader.py

@@ -0,0 +1,28 @@
+import datetime
+
+from astroplan import Observer
+from astropy.wcs import WCS
+from pyobs.images import Image
+
+from astropy import units as u
+from pyobs_cloudcover.world_model.wcs_model import WCSModel
+
+
+class WCSModelLoader(object):
+    def __init__(self, file_path: str) -> None:
+        self._file_path = file_path
+
+    def __call__(self) -> WCSModel:
+        wcs_image = Image.from_file(self._file_path)
+
+        wcs = WCS(header=wcs_image.header)
+        config_obs_time = datetime.datetime.fromisoformat(wcs_image.header["DATE-OBS"])
+
+        observer = Observer(
+            latitude=wcs_image.header["LATITUDE"] * u.deg,
+            longitude=wcs_image.header["LONGITUD"] * u.deg,
+            elevation=wcs_image.header["HEIGHT"] * u.m,
+        )
+
+        return WCSModel(wcs=wcs, observer=observer, config_obs_time=config_obs_time)
+

pyproject.toml

@@ -8,8 +8,8 @@ packages = [{include = "pyobs_cloudcover"}]
 
 [tool.poetry.dependencies]
 python = ">=3.11,<3.12"
-pyobs-core = {version = ">=1.13.0", extras = ["full"]}
-influxdb-client = "^1.42.0"
+#pyobs-core = {version = ">=1.13.0", extras = ["full"]}
+#influxdb-client = "^1.42.0"
 
 [tool.poetry.group.dev.dependencies]
 pytest = "^7.4.0"

src/alt_az_coord.rs

@@ -129,6 +129,16 @@ mod tests
         assert!(to.distance(&from) - PI/2.0 < 1e-10);
     }
 
+    #[test]
+    fn test_distance_alt_eq()
+    {
+        let from = AltAzCoord::new(PI/2.0, 0.0);
+        let to1 = AltAzCoord::new(PI/4.0, 0.0);
+        let to2 = AltAzCoord::new(PI/4.0, PI);
+
+        assert!(from.distance(&to2) - from.distance(&to1) < 1e-10);
+    }
+
     #[test]
     fn test_distance_az()
     {

src/average.rs

@@ -1,6 +1,7 @@
 use pyo3::prelude::*;
 
 #[pyclass]
+#[derive(Debug)]
 pub struct Average
 {
     #[pyo3(get)]
@@ -30,16 +31,16 @@ impl Average
 
     pub fn calc_weighted(values: &Vec<f64>, weights: &Vec<f64>) -> Option<Average>
     {
-        let average = Self::calc(values)?;
 
         let norm: f64 = weights.iter().fold(0.0, |acc, x| acc + x);
 
         let sum: f64 = values.iter().zip(weights).fold(0.0, |acc, (x, w)| acc + x * w);
         let weighted_average: f64 = sum/norm;
 
-        let weighted_std: f64 = average.get_std() * weights.iter().fold(0.0, |acc, x| acc + x.powf(2.0)).sqrt();
+        let weighted_std: f64 = values.iter().zip(weights).map(|(x, w)| (weighted_average - x).powi(2)/w.powi(2)).sum();
+        let std: f64 = (weighted_std/((values.len() as f64 - 1.0) * weights.iter().map(|x| 1.0/x.powi(2)).sum::<f64>())).sqrt();
 
-        Some(Average {value: weighted_average, std: weighted_std})
+        Some(Average {value: weighted_average, std})
     }
 
     pub fn get_value(&self) -> f64
@@ -53,6 +54,13 @@ impl Average
     }
 }
 
+impl PartialEq for Average
+{
+    fn eq(&self, other: &Self) -> bool {
+        self.value == other.value && self.std == other.std
+    }
+}
+
 #[cfg(test)]
 mod tests
 {