Atmospheric correction
774161
def atmos(
    ctx,
    atmo,
    contrast,
    bias,
    jobs,
    out_dtype,
    src_path,
    dst_path,
    creation_options,
    as_color,
):
    
    if as_color:
        click.echo(
            "rio color {} {} {}".format(
                src_path, dst_path, simple_atmo_opstring(atmo, contrast, bias)
            )
        )
        exit(0)

    with rasterio.open(src_path) as src:
        opts = src.profile.copy()
        windows = [(window, ij) for ij, window in src.block_windows()]

    opts.update(**creation_options)
    opts["transform"] = guard_transform(opts["transform"])

    out_dtype = out_dtype if out_dtype else opts["dtype"]
    opts["dtype"] = out_dtype

    args = {"atmo": atmo, "contrast": contrast, "bias": bias, "out_dtype": out_dtype}

    jobs = check_jobs(jobs)

    if jobs > 1:
        with riomucho.RioMucho(
            [src_path],
            dst_path,
            atmos_worker,
            windows=windows,
            options=opts,
            global_args=args,
            mode="manual_read",
        ) as mucho:
            mucho.run(jobs)
    else:
        with rasterio.open(dst_path, "w", **opts) as dest:
            with rasterio.open(src_path) as src:
                rasters = [src]
                for window, ij in windows:
                    arr = atmos_worker(rasters, window, ij, args)
                    dest.write(arr, window=window)
427828
def re_normalize_flux(self, kwargs_ps, norm_factor):
        
        for i, model in enumerate(self.point_source_type_list):
            if model == 'UNLENSED':
                kwargs_ps[i]['point_amp'] *= norm_factor
            elif model in ['LENSED_POSITION', 'SOURCE_POSITION']:
                if self._fixed_magnification_list[i] is True:
                    kwargs_ps[i]['source_amp'] *= norm_factor
                else:
                    kwargs_ps[i]['point_amp'] *= norm_factor
        return kwargs_ps
1009132
def simple_atmo_opstring(haze, contrast, bias):
    """"""
    gamma_b = 1 - haze
    gamma_g = 1 - (haze / 3.0)
    ops = (
        "gamma g {gamma_g}, " "gamma b {gamma_b}, " "sigmoidal rgb {contrast} {bias}"
    ).format(gamma_g=gamma_g, gamma_b=gamma_b, contrast=contrast, bias=bias)
    return ops
1216669
function pixelToLatitude (y, zoom) {
  var latRad = Math.atan(Math.sinh(Math.PI * (1 - 2 * y / zScale(zoom))))
  return toDegrees(latRad)
}
377116
def refract(alt_degrees, temperature_C, pressure_mbar):
    """"""
    alt = alt_degrees
    while True:
        alt1 = alt
        alt = alt_degrees + refraction(alt, temperature_C, pressure_mbar)
        converged = abs(alt - alt1) <= 3.0e-5
        if converged.all():
            break
    return alt
775392
def update_view_bounds(self):
        """"""
        self.view_bounds.left = self.pan.X - self.world_center.X
        self.view_bounds.top = self.pan.Y - self.world_center.Y
        self.view_bounds.width = self.world_center.X * 2
        self.view_bounds.height = self.world_center.Y * 2
1122479
def enhance_dark_holes(image, min_radius, max_radius, mask=None):
    '''
    '''
    #
    # Do 4-connected erosion
    #
    se = np.array([[False, True, False],
                   [True, True, True],
                   [False, True, False]])
    #
    # Invert the intensities
    #
    inverted_image = image.max() - image
    previous_reconstructed_image = inverted_image
    eroded_image = inverted_image
    smoothed_image = np.zeros(image.shape)
    for i in range(max_radius+1):
        eroded_image = grey_erosion(eroded_image, mask=mask, footprint = se)
        reconstructed_image = grey_reconstruction(eroded_image, inverted_image,
                                                  footprint = se)
        output_image = previous_reconstructed_image - reconstructed_image
        if i >= min_radius:
            smoothed_image = np.maximum(smoothed_image,output_image)
        previous_reconstructed_image = reconstructed_image
    return smoothed_image
720290
def enforce_vertical_symmetry(pixmap):
    ''''''
    mirror = []
    for item in pixmap:
        y = item[0]
        x = item[1]

        if x <= IMAGE_APEX:
            diff_x = diff(x, IMAGE_APEX)
            mirror.append((y, x + (2 * diff_x) - 1))

        if x > IMAGE_APEX:
            diff_x = diff(x, IMAGE_APEX)
            mirror.append((y, x - (2 * diff_x) - 1))

    return mirror + pixmap
1169120
def _calibrate(self, radiance, coefs, channel, calibration):
        """"""
        if self._is_vis(channel):
            if not calibration == 'reflectance':
                raise ValueError('Cannot calibrate VIS channel to '
                                 '{}'.format(calibration))
            return self._calibrate_vis(radiance=radiance, k=coefs['k'])
        else:
            if not calibration == 'brightness_temperature':
                raise ValueError('Cannot calibrate IR channel to '
                                 '{}'.format(calibration))

            # Since the scanline-detector assignment is unknown, use the average
            # coefficients for all scanlines.
            mean_coefs = {'a': np.array(coefs['a']).mean(),
                          'b': np.array(coefs['b']).mean(),
                          'n': np.array(coefs['n']).mean(),
                          'btmin': coefs['btmin'],
                          'btmax': coefs['btmax']}
            return self._calibrate_ir(radiance=radiance, coefs=mean_coefs)
262365
def _gti_dirint_lt_90(poa_global, aoi, aoi_lt_90, solar_zenith, solar_azimuth,
                      times, surface_tilt, surface_azimuth, pressure=101325.,
                      use_delta_kt_prime=True, temp_dew=None, albedo=.25,
                      model='perez', model_perez='allsitescomposite1990',
                      max_iterations=30):
    
    I0 = get_extra_radiation(times, 1370, 'spencer')
    cos_zenith = tools.cosd(solar_zenith)
    # I0h as in Marion 2015 eqns 1, 3
    I0h = I0 * np.maximum(0.065, cos_zenith)

    airmass = atmosphere.get_relative_airmass(solar_zenith, model='kasten1966')
    airmass = atmosphere.get_absolute_airmass(airmass, pressure)

    # these coeffs and diff variables and the loop below
    # implement figure 1 of Marion 2015

    # make coeffs that is at least 30 elements long so that all
    # coeffs can be assigned as specified in Marion 2015.
    # slice below will limit iterations if necessary
    coeffs = np.empty(max(30, max_iterations))
    coeffs[0:3] = 1
    coeffs[3:10] = 0.5
    coeffs[10:20] = 0.25
    coeffs[20:] = 0.125
    coeffs = coeffs[:max_iterations]  # covers case where max_iterations < 30

    # initialize diff
    diff = pd.Series(9999, index=times)
    best_diff = diff

    # initialize poa_global_i
    poa_global_i = poa_global

    for iteration, coeff in enumerate(coeffs):

        # test if difference between modeled GTI and
        # measured GTI (poa_global) is less than 1 W/m^2
        # only test for aoi less than 90 deg
        best_diff_lte_1 = best_diff <= 1
        best_diff_lte_1_lt_90 = best_diff_lte_1[aoi_lt_90]
        if best_diff_lte_1_lt_90.all():
            # all aoi < 90 points have a difference <= 1, so break loop
            break

        # calculate kt and DNI from GTI
        kt = clearness_index(poa_global_i, aoi, I0)  # kt from Marion eqn 2
        disc_dni = np.maximum(_disc_kn(kt, airmass)[0] * I0, 0)
        kt_prime = clearness_index_zenith_independent(kt, airmass)
        # dirint DNI in Marion eqn 3
        dni = _dirint_from_dni_ktprime(disc_dni, kt_prime, solar_zenith,
                                       use_delta_kt_prime, temp_dew)

        # calculate DHI using Marion eqn 3 (identify 1st term on RHS as GHI)
        # I0h has a minimum zenith projection, but multiplier of DNI does not
        ghi = kt * I0h                  # Kt * I0 * max(0.065, cos(zen))
        dhi = ghi - dni * cos_zenith    # no cos(zen) restriction here

        # following SSC code
        dni = np.maximum(dni, 0)
        ghi = np.maximum(ghi, 0)
        dhi = np.maximum(dhi, 0)

        # use DNI and DHI to model GTI
        # GTI-DIRINT uses perez transposition model, but we allow for
        # any model here
        all_irrad = get_total_irradiance(
            surface_tilt, surface_azimuth, solar_zenith, solar_azimuth,
            dni, ghi, dhi, dni_extra=I0, airmass=airmass,
            albedo=albedo, model=model, model_perez=model_perez)

        gti_model = all_irrad['poa_global']

        # calculate new diff
        diff = gti_model - poa_global

        # determine if the new diff is smaller in magnitude
        # than the old diff
        diff_abs = diff.abs()
        smallest_diff = diff_abs < best_diff

        # save the best differences
        best_diff = diff_abs.where(smallest_diff, best_diff)

        # on first iteration, the best values are the only values
        if iteration == 0:
            best_ghi = ghi
            best_dni = dni
            best_dhi = dhi
            best_kt_prime = kt_prime
        else:
            # save new DNI, DHI, DHI if they provide the best consistency
            # otherwise use the older values.
            best_ghi = ghi.where(smallest_diff, best_ghi)
            best_dni = dni.where(smallest_diff, best_dni)
            best_dhi = dhi.where(smallest_diff, best_dhi)
            best_kt_prime = kt_prime.where(smallest_diff, best_kt_prime)

        # calculate adjusted inputs for next iteration. Marion eqn 4
        poa_global_i = np.maximum(1.0, poa_global_i - coeff * diff)
    else:
        # we are here because we ran out of coeffs to loop over and
        # therefore we have exceeded max_iterations
        import warnings
        failed_points = best_diff[aoi_lt_90][~best_diff_lte_1_lt_90]
        warnings.warn(
            ('%s points failed to converge after %s iterations. best_diff:\n%s'
             % (len(failed_points), max_iterations, failed_points)),
            RuntimeWarning)

    # return the best data, whether or not the solution converged
    return best_ghi, best_dni, best_dhi, best_kt_prime
258837
def offset_to_sky(skydir, offset_lon, offset_lat,
                  coordsys='CEL', projection='AIT'):
    """"""

    offset_lon = np.array(offset_lon, ndmin=1)
    offset_lat = np.array(offset_lat, ndmin=1)

    w = create_wcs(skydir, coordsys, projection)
    pixcrd = np.vstack((offset_lon, offset_lat)).T

    return w.wcs_pix2world(pixcrd, 0)
454065
private static void adjustUncalibrated(DMatrixRMaj rectifyLeft, DMatrixRMaj rectifyRight,
										   RectangleLength2D_F64 bound, double scale) {
		// translation
		double deltaX = -bound.x0*scale;
		double deltaY = -bound.y0*scale;

		// adjustment matrix
		SimpleMatrix A = new SimpleMatrix(3,3,true,new double[]{scale,0,deltaX,0,scale,deltaY,0,0,1});
		SimpleMatrix rL = SimpleMatrix.wrap(rectifyLeft);
		SimpleMatrix rR = SimpleMatrix.wrap(rectifyRight);

		rectifyLeft.set(A.mult(rL).getDDRM());
		rectifyRight.set(A.mult(rR).getDDRM());
	}
157569
def _get_aeff_corrections(intensity_ratio, mask):
        
        nebins = len(intensity_ratio.data)
        aeff_corrections = np.zeros((nebins))
        for i in range(nebins):
            bright_pixels_intensity = intensity_ratio.data[i][mask.data[i]]
            mean_bright_pixel = bright_pixels_intensity.mean()
            aeff_corrections[i] = 1. / mean_bright_pixel

        print("Aeff correction: ", aeff_corrections)
        return aeff_corrections
1261755
private boolean refineWithBundleAdjustment(SceneObservations observations) {
		if( scaleSBA ) {
			scaler.applyScale(structure,observations);
		}

		sba.setVerbose(verbose,verboseLevel);
		sba.setParameters(structure,observations);
		sba.configure(converge.ftol,converge.gtol,converge.maxIterations);

		if( !sba.optimize(structure) ) {
			return false;
		}

		if( scaleSBA ) {
			// only undo scaling on camera matrices since observations are discarded
			for (int i = 0; i < structure.views.length; i++) {
				DMatrixRMaj P = structure.views[i].worldToView;
				scaler.pixelScaling.get(i).remove(P,P);
			}
			scaler.undoScale(structure,observations);
		}
		return true;
	}
974240
def calibrate(self,
                  dataset_id,
                  pre_launch_coeffs=False,
                  calib_coeffs=None):
        
        tic = datetime.now()

        if calib_coeffs is None:
            calib_coeffs = {}

        units = {'reflectance': '%',
                 'brightness_temperature': 'K',
                 'counts': '',
                 'radiance': 'W*m-2*sr-1*cm ?'}

        if dataset_id.name in ("3a", "3b") and self._is3b is None:
            # Is it 3a or 3b:
            is3b = np.expand_dims(
                np.bitwise_and(
                    np.right_shift(self._data['scnlinbit'], 0), 1) == 1, 1)
            self._is3b = np.repeat(is3b,
                                   self._data['hrpt'][0].shape[0], axis=1)

        try:
            vis_idx = ['1', '2', '3a'].index(dataset_id.name)
            ir_idx = None
        except ValueError:
            vis_idx = None
            ir_idx = ['3b', '4', '5'].index(dataset_id.name)

        if vis_idx is not None:
            coeffs = calib_coeffs.get('ch' + dataset_id.name)
            ds = create_xarray(
                _vis_calibrate(self._data,
                               vis_idx,
                               dataset_id.calibration,
                               pre_launch_coeffs,
                               coeffs,
                               mask=(dataset_id.name == '3a' and self._is3b)))
        else:
            ds = create_xarray(
                _ir_calibrate(self._header,
                              self._data,
                              ir_idx,
                              dataset_id.calibration,
                              mask=(dataset_id.name == '3b' and
                                    np.logical_not(self._is3b))))

        if dataset_id.name == '3a' and np.all(np.isnan(ds)):
            raise ValueError("Empty dataset for channel 3A")
        if dataset_id.name == '3b' and np.all(np.isnan(ds)):
            raise ValueError("Empty dataset for channel 3B")

        ds.attrs['units'] = units[dataset_id.calibration]
        ds.attrs.update(dataset_id._asdict())

        logger.debug("Calibration time %s", str(datetime.now() - tic))

        return ds
1051217
function(){
            var x = this.x,
                y = this.y,
                fieldRange = this.fovFieldRange,
                direction = this.fovDirection,
                maxViewDistance = this.fovMaxViewDistance;
            this.fov.update(x, y, fieldRange, direction, maxViewDistance, this);
        }
883989
function boxBlurH(scl, tcl, w, h, r) {
    var iarr = 1 / (r+r+1); // averaging adjustment parameter
    for (var i = 0; i < h; i++) {
        var ti = i * w, // current target index
            li = ti, // current left side of the examined range
            ri = ti + r, // current right side of the examined range
            fv = scl[ti], // first value in the row
            lv = scl[ti + w - 1], // last value in the row
            val = (r+1) * fv, // target value, accumulated over examined points
            j;
        // Sum the source values in the box
        for (j = 0; j < r; j++) { val += scl[ti + j]; }
        // Compute the target value by taking the average of the surrounding
        // values. This is done by adding the deviations so far and adjusting,
        // accounting for the edges by extending the first and last values.
        for (j = 0  ; j <= r ; j++) { val += scl[ri++] - fv       ; tcl[ti++] = val*iarr; }
        for (j = r+1; j < w-r; j++) { val += scl[ri++] - scl[li++]; tcl[ti++] = val*iarr; }
        for (j = w-r; j < w  ; j++) { val += lv        - scl[li++]; tcl[ti++] = val*iarr; }
    }
}