Added plottin test script

.gitignore

@@ -70,3 +70,5 @@ docs/_build/
 
 # Pyenv
 .python-version
+
+*.sh

plotting.py

@@ -0,0 +1,93 @@
+import sys
+import datetime
+
+import numpy as np
+from astroplan import Observer
+from astropy.visualization import PercentileInterval
+import astropy.units as u
+from pyobs.images import Image
+
+from pyobs_cloudcover.pipeline.night.altaz_map_generator.altaz_map_generator import AltAzMapGenerator
+from pyobs_cloudcover.pipeline.night.catalog.altaz_catalog_loader import AltAzCatalogLoader
+from pyobs_cloudcover.pipeline.night.catalog.catalog_constructor import CatalogConstructor
+from pyobs_cloudcover.pipeline.night.cloud_coverage_calculator.coverage_calculator import CoverageCalculator
+from pyobs_cloudcover.pipeline.night.cloud_coverage_calculator.coverage_change_calculator import \
+    CoverageChangeCalculator
+from pyobs_cloudcover.pipeline.night.cloud_coverage_calculator.coverage_info_calculator import CoverageInfoCalculator
+from pyobs_cloudcover.pipeline.night.cloud_coverage_calculator.zenith_masker import ZenithMasker
+from pyobs_cloudcover.pipeline.night.cloud_map_generator.cloud_map_generator import CloudMapGenerator
+from pyobs_cloudcover.pipeline.night.pipeline import NightPipeline
+from pyobs_cloudcover.pipeline.night.preprocessor.background_remover import BackgroundRemover
+from pyobs_cloudcover.pipeline.night.preprocessor.image_binner import ImageBinner
+from pyobs_cloudcover.pipeline.night.preprocessor.image_masker import ImageMasker
+from pyobs_cloudcover.pipeline.night.preprocessor.preprocessor import Preprocessor
+from pyobs_cloudcover.pipeline.night.star_reverse_matcher.detector.sigma_treshhold_detector import \
+    SigmaThresholdDetector
+from pyobs_cloudcover.pipeline.night.star_reverse_matcher.star_reverse_matcher import StarReverseMatcher
+from pyobs_cloudcover.pipeline.night.star_reverse_matcher.window import ImageWindow
+from pyobs_cloudcover.world_model.wcs_model_loader import WCSModelLoader
+
+
+def main() -> None:
+    import matplotlib.pyplot as plt
+
+    wcs_file, catalog_file, good_image_file, bad_image_file = sys.argv[1:]
+
+    good_image = Image.from_file(good_image_file)
+    good_obs_time = datetime.datetime.fromisoformat(good_image.header["DATE-OBS"])
+
+    pipeline = build_pipeline(wcs_file, catalog_file, good_image.data.shape)
+    good_coverage_info = pipeline(good_image.data, good_obs_time)
+
+    bad_image = Image.from_file(bad_image_file)
+    bad_obs_time = datetime.datetime.fromisoformat(bad_image.header["DATE-OBS"])
+    bad_coverage_info = pipeline(bad_image.data, bad_obs_time)
+
+    cutoff = 5.8
+
+    plt.figure(figsize=(20, 10))
+    plt.subplot(121)
+    plt.imshow(PercentileInterval(99)(ImageBinner(2)(good_image.data)), cmap="gray")
+    plt.imshow(good_coverage_info.cloud_cover_map, alpha=(good_coverage_info.cloud_cover_map < cutoff).astype(np.float_))
+    plt.colorbar()
+
+    plt.subplot(122)
+    plt.imshow(PercentileInterval(99)(ImageBinner(2)(bad_image.data)), cmap="gray")
+    plt.imshow(bad_coverage_info.cloud_cover_map, alpha=(bad_coverage_info.cloud_cover_map < cutoff).astype(np.float_))
+    plt.colorbar()
+
+    plt.show()
+
+def build_pipeline(wcs_file, catalog_file, image_shape) -> NightPipeline:
+    observer = Observer(latitude=51.559299 * u.deg, longitude=9.945472 * u.deg, elevation=201 * u.m)
+
+    wcs_model_factory = WCSModelLoader(wcs_file)
+    model = wcs_model_factory()
+
+    mask = ImageMasker(np.ones(image_shape).astype(np.bool_))
+    binner = ImageBinner(2)
+    background_remover = BackgroundRemover(3.0, (5, 5))
+    preprocessor = Preprocessor(mask, binner, background_remover)
+
+    altaz_catalog_loader = AltAzCatalogLoader.from_csv(catalog_file)
+    catalog_constructor = CatalogConstructor(altaz_catalog_loader, model, observer, 30.0, 7.0, 4.0)
+
+    altaz_list_generator = AltAzMapGenerator(model, 30.0)
+
+    reverse_matcher = StarReverseMatcher(SigmaThresholdDetector(3.0), ImageWindow(6.0))
+
+    cloud_map_gem = CloudMapGenerator(7.0)
+
+    coverage_calculator = CoverageCalculator(0.5)
+    coverage_change_calculator = CoverageChangeCalculator()
+    zenith_masker = ZenithMasker(80)
+    cloud_coverage_info_calculator = CoverageInfoCalculator(coverage_calculator, coverage_change_calculator,
+                                                            zenith_masker)
+
+    pipeline = NightPipeline(preprocessor, catalog_constructor, altaz_list_generator, reverse_matcher, cloud_map_gem,
+                             cloud_coverage_info_calculator)
+
+    return pipeline
+
+if __name__ == '__main__':
+    main()

pyobs_cloudcover/pipeline/night/altaz_map_generator/altaz_map_generator.py

@@ -2,7 +2,6 @@
 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
@@ -34,4 +33,4 @@ def _alt_az_convert(self, alt_az_coords: Tuple[float, float]) -> Optional[AltAzC
         if alt_az_coords[0] < np.deg2rad(self._limiting_altitude):
             return None
         else:
-            return AltAzCoord(*alt_az_coords)
+            return AltAzCoord(*alt_az_coords)

pyobs_cloudcover/pipeline/night/catalog/catalog_constructor.py

@@ -17,8 +17,8 @@ def __init__(self, catalog_loader: AltAzCatalogLoader, model: WorldModel, observ
         self._v_mag_limit = v_mag_limit
         self._distance_limit = distance_limit
 
-    def __call__(self, time: datetime.datetime, height: int, width: int) -> PixelCatalog:
-        altaz_catalog = self._catalog_loader(self._observer, time)
+    def __call__(self, obs_time: datetime.datetime, height: int, width: int) -> PixelCatalog:
+        altaz_catalog = self._catalog_loader(self._observer, obs_time)
         altaz_catalog.filter_alt(self._alt_limit)
         altaz_catalog.filter_v_mag(self._v_mag_limit)
 

pyobs_cloudcover/pipeline/night/catalog/pixel_catalog.py

@@ -1,5 +1,6 @@
 from __future__ import annotations
 
+import datetime
 from typing import cast
 
 import numpy as np

pyobs_cloudcover/pipeline/night/cloud_coverage_calculator/cloud_info_calculator_factory.py

@@ -16,7 +16,7 @@ def __init__(self, options: CloudInfoCalculatorOptions, model: WorldModel):
     def __call__(self) -> CoverageInfoCalculator:
         coverage_calculator = CoverageCalculator(self._options.cloud_threshold)
         coverage_change_calculator = CoverageChangeCalculator()
-        zenith_masker = ZenithMasker(self._options.altitude_limit, self._model)
+        zenith_masker = ZenithMasker(self._options.altitude_limit)
         cloud_coverage_info_calculator = CoverageInfoCalculator(coverage_calculator, coverage_change_calculator,
                                                                 zenith_masker)
 

pyobs_cloudcover/pipeline/night/cloud_coverage_calculator/coverage_info_calculator.py

@@ -1,8 +1,10 @@
 import datetime
 from copy import copy
+from typing import List, Optional
 
 import numpy as np
 import numpy.typing as npt
+from cloudmap_rs import AltAzCoord
 
 from pyobs_cloudcover.cloud_coverage_info import CloudCoverageInfo
 from pyobs_cloudcover.pipeline.night.cloud_coverage_calculator.coverage_calculator import CoverageCalculator
@@ -17,9 +19,9 @@ def __init__(self, coverage_calculator: CoverageCalculator, coverage_change_calc
         self._coverage_change_calculator = coverage_change_calculator
         self._zenith_masker = zenith_masker
 
-    def __call__(self, cloud_map: npt.NDArray[np.float_], obs_time: datetime.datetime) -> CloudCoverageInfo:
+    def __call__(self, cloud_map: npt.NDArray[np.float_], obs_time: datetime.datetime, alt_az_list: List[List[Optional[AltAzCoord]]]) -> CloudCoverageInfo:
         change = self._coverage_change_calculator(cloud_map)
         coverage = self._coverage_calculator(cloud_map)
-        zenith_coverage = self._coverage_calculator(self._zenith_masker(cloud_map))
+        zenith_coverage = self._coverage_calculator(self._zenith_masker(cloud_map, alt_az_list))
 
         return CloudCoverageInfo(copy(cloud_map), coverage, zenith_coverage, change, obs_time)

pyobs_cloudcover/pipeline/night/cloud_coverage_calculator/zenith_masker.py

@@ -1,25 +1,19 @@
 from copy import copy
+from typing import List, Optional
 
 import numpy as np
 import numpy.typing as npt
-
-from pyobs_cloudcover.world_model import WorldModel
+from cloudmap_rs import AltAzCoord
 
 
 class ZenithMasker(object):
-    def __init__(self, altitude: float, model: WorldModel) -> None:
-        self._zenith = np.array(model.altaz_to_pix(0.0, 0))
-        offset = np.array(model.altaz_to_pix(np.deg2rad(altitude), 0.0))
-
-        self._radius = np.sqrt(np.sum(np.square(self._zenith - offset)))
+    def __init__(self, altitude: float) -> None:
+        self._altitude = altitude
 
-    def __call__(self, image: npt.NDArray[np.float_]) -> npt.NDArray[np.float_]:
-        ny, nx = image.shape
-        x, y = np.arange(0, nx), np.arange(0, ny)
-        x_coordinates, y_coordinates = np.meshgrid(x, y)
+    def __call__(self, image: npt.NDArray[np.float_], alt_az_list: List[List[Optional[AltAzCoord]]]) -> npt.NDArray[np.float_]:
+        mask = [[entry.alt > self._altitude if entry is not None else False for entry in row] for row in alt_az_list]
 
-        circ_mask = (x_coordinates - self._zenith[0])**2 + (y_coordinates - self._zenith[1])**2 <= self._radius**2
         masked_image = copy(image)
-        masked_image[~circ_mask] = None
+        masked_image[mask] = None
 
         return masked_image

pyobs_cloudcover/pipeline/night/cloud_map_generator/cloud_map_generator.py

@@ -16,19 +16,11 @@ def __call__(self, catalog: PixelCatalog, matches: List[bool], alt_az_image_list
 
         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)
-
-        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) / (std_vis_map - no_vis_map)
+        return vis_map
 
     @staticmethod
     def _convert_average_to_value(average_map: List[List[Average]]) -> npt.NDArray[np.float_]:

pyobs_cloudcover/pipeline/night/pipeline.py

@@ -36,7 +36,16 @@ 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)
 
+        '''
+        import matplotlib.pyplot as plt
+        from astropy.visualization import PercentileInterval
+        plt.imshow(PercentileInterval(99)(preprocessed_image), cmap="gray")
+        plt.scatter(catalog.px[matches], catalog.py[matches], alpha=0.2, c="green", s=6)
+        plt.scatter(catalog.px[~np.array(matches)], catalog.py[~np.array(matches)], alpha=0.2, c="red", s=6)
+        plt.show()
+        '''
+
         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)
+        return self._coverage_info_calculator(cloud_map, obs_time, alt_az_list)

pyobs_cloudcover/pipeline/night/preprocessor/preprocessor.py

@@ -13,8 +13,9 @@ def __init__(self, masker: ImageMasker, binner: ImageBinner, background_remover:
         self._background_remover = background_remover
 
     def __call__(self, image: npt.NDArray[np.float_]) -> npt.NDArray[np.float_]:
-        masked_image = self._masker(image)
+        flipped_image = np.flip(image, axis=0)
+        masked_image = self._masker(flipped_image)
         binned_image = self._binner(masked_image)
         processed_image = self._background_remover(binned_image)
 
-        return processed_image
+        return binned_image

pyobs_cloudcover/world_model/hour_angle_converter.py

@@ -1,5 +1,8 @@
 import datetime
+from typing import Union, cast
 
+import numpy as np
+import numpy.typing as npt
 from astroplan import Observer
 from astropy.coordinates import Angle
 
@@ -8,8 +11,8 @@ class HourAngleConverter(object):
     def __init__(self, observer: Observer):
         self._observer = observer
 
-    def __call__(self, ra: Angle, obs_time: datetime.datetime) -> float:
+    def __call__(self, ra: Angle, obs_time: datetime.datetime) -> Union[float, npt.NDArray[np.float_]]:
         lst = self._observer.local_sidereal_time(obs_time)
         self._lst = Angle(lst, unit="hourangle")
 
-        return float((self._lst - ra).rad)
+        return cast(Union[float, npt.NDArray[np.float_]], (self._lst - ra).rad)

pyobs_cloudcover/world_model/wcs_model_loader.py

@@ -1,10 +1,10 @@
 import datetime
 
 from astroplan import Observer
+from astropy import units as u
+from astropy.io import fits
 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
 
 
@@ -13,7 +13,7 @@ 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_image = fits.open(self._file_path)[0]
 
         wcs = WCS(header=wcs_image.header)
         config_obs_time = datetime.datetime.fromisoformat(wcs_image.header["DATE-OBS"])

pyproject.toml

@@ -48,7 +48,8 @@ test = [
     "pytest-cov>=4.1.0",
     "pytest-mock>=3.11.1",
     "pytest-asyncio>=0.21.1",
-    "mypy>=1.9.0"
+    "mypy>=1.9.0",
+    "matplotlib"
 ]
 
 [tool.maturin]