How to use the pyorbital.astronomy function in pyorbital

def __call__(self, projectables, *args, **kwargs):
        (area, lat, lon) = self._get_area_lat_lon(projectables)

        sza = pyorbital.astronomy.sun_zenith_angle(
                projectables[0].start_time, lon, lat)

        maskproj = self._convert_xr_to_ma(projectables)

        flsinput = {'ir108': maskproj[1],
                    'ir039': maskproj[0],
                    'sza': sza,
                    'lat': lat,
                    'lon': lon,
                    'time': projectables[0].start_time
                    }

        # Compute fog mask
        flsalgo = NightFogLowStratusAlgorithm(**flsinput)
        fls, mask = flsalgo.run()

self.info['sun_zen_corrected'] to the original channel, so
        this can be used as a pointer to the corrected data.
        '''

        try:
            from pyorbital import astronomy
        except ImportError:
            LOG.warning("Could not load pyorbital.astronomy")
            return None

        if lons is None or lats is None:
            # Read coordinates
            lons, lats = self.area.get_lonlats()
    
        # Calculate Sun zenith angles and the cosine
        zen_angles = astronomy.sun_zenith_angle(time_slot, 
                                                lons, lats)
        cos_zen = astronomy.cos_zen(time_slot, lons, lats)

        # Copy the channel
        new_ch = copy.deepcopy(self)
        # Update the name
        if name is None:
            new_ch.name += '_SZC'
        else:
            new_ch.name = name

        if mode == 'cos':
            # Cosine correction
            lim_y, lim_x = np.where(zen_angles < limit)
            new_ch.data[lim_y, lim_x] /= cos_zen[lim_y, lim_x]
            # Use constant value (the limit) for larger zenith

fog_mask = fog_mask | ice_mask
    filters['ice'] = np.nansum(fog_mask) - prev
    prev += filters['ice']

    # Thin cirrus is detected by means of the split-window IR channel
    # brightness temperature difference (T10.8 –T12.0 ). This difference is
    # compared to a threshold dynamically interpolated from a lookup table
    # based on satellite zenith angle and brightness temperature at 10.8 μm
    # (Saunders and Kriebel, 1988)
    chn120_ma = np.ma.masked_where(cloud_mask | snow_mask | ice_mask, chn120)
    chn108_ma = np.ma.masked_where(cloud_mask | snow_mask | ice_mask, chn108)

    bt_diff = chn108 - chn120

    # Calculate sun zenith angles
    sza = astronomy.sun_zenith_angle(time, lon, lat)

    minsza = np.min(sza)
    maxsza = np.max(sza)
    logger.debug("Found solar zenith angles from %s to %s°" % (minsza,
                                                               maxsza))

    # Calculate secant of sza
    # secsza = np.ma.masked_where(cloud_mask | snow_mask | ice_mask,
    #                             (1 / np.cos(np.deg2rad(sza))))
    secsza = 1 / np.cos(np.deg2rad(sza))

    # Lookup table for BT difference thresholds at certain sec(sun zenith
    # angles) and 10.8 μm BT
    lut = {260: {1.0: 0.55, 1.25: 0.60, 1.50: 0.65, 1.75: 0.90, 2.0: 1.10},
           270: {1.0: 0.58, 1.25: 0.63, 1.50: 0.81, 1.75: 1.03, 2.0: 1.13},
           280: {1.0: 1.30, 1.25: 1.61, 1.50: 1.88, 1.75: 2.14, 2.0: 2.30},

fog_mask = fog_mask | ice_mask
    filters['ice'] = np.nansum(fog_mask) - prev
    prev += filters['ice']

    # Thin cirrus is detected by means of the split-window IR channel
    # brightness temperature difference (T10.8 –T12.0 ). This difference is
    # compared to a threshold dynamically interpolated from a lookup table
    # based on satellite zenith angle and brightness temperature at 10.8 μm
    # (Saunders and Kriebel, 1988)
    chn120_ma = np.ma.masked_where(cloud_mask | snow_mask | ice_mask, chn120)
    chn108_ma = np.ma.masked_where(cloud_mask | snow_mask | ice_mask, chn108)

    bt_diff = chn108 - chn120

    # Calculate sun zenith angles
    sza = astronomy.sun_zenith_angle(time, lon, lat)

    minsza = np.min(sza)
    maxsza = np.max(sza)
    logger.debug("Found solar zenith angles from %s to %s°" % (minsza,
                                                               maxsza))

    # Calculate secant of sza
    # secsza = np.ma.masked_where(cloud_mask | snow_mask | ice_mask,
    #                             (1 / np.cos(np.deg2rad(sza))))
    secsza = 1 / np.cos(np.deg2rad(sza))

    # Lookup table for BT difference thresholds at certain sec(sun zenith
    # angles) and 10.8 μm BT
    lut = {260: {1.0: 0.55, 1.25: 0.60, 1.50: 0.65, 1.75: 0.90, 2.0: 1.10},
           270: {1.0: 0.58, 1.25: 0.63, 1.50: 0.81, 1.75: 1.03, 2.0: 1.13},
           280: {1.0: 1.30, 1.25: 1.61, 1.50: 1.88, 1.75: 2.14, 2.0: 2.30},

Args:
            | ir108 (:obj:`ndarray`): Array for the 10.8 μm channel.
            | ir087 (:obj:`ndarray`): Array for the 8.7 μm channel.
            | ir120 (:obj:`ndarray`): Array for the 12.0 μm channel.
            | time (:obj:`datetime`): Datetime object for the satellite scence.
            | lat (:obj:`ndarray`): Array of latitude values.
            | lon (:obj:`ndarray`): Array of longitude values.

        Returns:
            Filter image and filter mask.
        """
        logger.info("Applying Cirrus Filter")
        # Get infrared channel difference
        self.bt_diff = self.ir108 - self.ir120
        # Calculate sun zenith angles
        sza = astronomy.sun_zenith_angle(self.time, self.lon, self.lat)
        minsza = np.min(sza)
        maxsza = np.max(sza)
        logger.debug("Found solar zenith angles from %s to %s°" % (minsza,
                                                                   maxsza))
        # Calculate secant of sza
        secsza = 1 / np.cos(np.deg2rad(sza))

        # Apply lut to BT and sza values
        # Vectorize LUT functions for numpy arrays
        vfind_nearest_lut_sza = np.vectorize(self.find_nearest_lut_sza)
        vfind_nearest_lut_bt = np.vectorize(self.find_nearest_lut_bt)
        vapply_lut = np.vectorize(self.apply_lut)

        secsza_lut = vfind_nearest_lut_sza(secsza)
        chn108_ma_lut = vfind_nearest_lut_bt(self.ir108)

def _get_day_percentage(self, refl_swath):
        swath_name = refl_swath["swath_definition"]["swath_name"]
        if swath_name not in self._day_percentage:
            from pyorbital import astronomy
            lons = refl_swath["swath_definition"].get_longitude_array()
            lats = refl_swath["swath_definition"].get_latitude_array()
            invalid_mask = refl_swath.get_data_mask()
            sza_data = astronomy.sun_zenith_angle(refl_swath["begin_time"], lons, lats)
            valid_day_mask = (sza_data < self.sza_threshold) & ~invalid_mask
            fraction_day = numpy.count_nonzero(valid_day_mask) / (float(sza_data.size) - numpy.count_nonzero(invalid_mask))
            self._day_percentage[swath_name] = fraction_day * 100.0
        else:
            LOG.debug("Using cached day percentage")
        return self._day_percentage[swath_name]