mc_saf: Monte Carlo spatially-resolved spherical albedo

In AlexCast/apsfs: Monte Carlo code for simulation of atmospheric Point Spread Functions and spatially-resolved Spherical Albedo Functions

Monte Carlo spatially-resolved spherical albedo

Description

Backward Monte Carlo Code to calculate the elastic scattering spatially-resolved spherical albedo function (integrated in bin area) of a vertically stratified but horizontally homogeneous atmosphere at monocromatic regime. The spatially-resolved spherical albedo is the angular integral of the bi-directional scattering-surface reflectance distribution function (BSSRDF). The spherical albedo is defined as the normalized irradiance backscattered by the atmosphere when the input irradiance at the bottom is isotropic (Vermote et al., 1997).

Usage

mc_saf(atm, geom, res, ext, np, mnw, cdf_aer, cdf_ray = NULL, press = NULL)

Arguments

atm	Atmosphere model, see Details.
geom	One of "annular", "sectorial" or "grid", see Details.
res	Resolution (km), [0,10]. See Details.
ext	Maximum spatial extent at geom resolution (km), see Details.
np	Number of photons to trace, [0,2^64 - 1].
mnw	Minimal weight to keep tracing photon package (unitless), [0,1].
cdf_aer	Aerosol scattering CDF LUT (unitless).
press	Atmospheric pressure at surface (mbar). Can be omitted if present on the attributes of atm. See Details.

Details

Vermote, E. F., Tanre, D., Deuze, J. L., Herman, M., Morcette, J.-J. (1997). Second Simulation of the Satellite Signal in the Solar Spectrum, 6S: an overview. IEEE Transactions on Geoscience and Remote Sensing, 35(3), 675–686.

This backward Monte Carlo radiative transfer code calculates the area integrals of the spatially-resolved spherical albedo of the atmosphere due to elastic scattering in a layered atmosphere. Turbulence effects and inelastic scattering are not included. Additional simplifications of the current surface reflectance (Bidirectional Reflectance Distribution Function effects not included). The surface is considered totally absorbing, such that multiple interaction between surface and atmosphere are not included.

The sensor position is at the origin, at the bottom of the atmosphere, and the zenith view angle is 180 degrees (looking at the zenith). The FOV is hemispherical.

The atmosphere profile atm should be a named data frame as the one created by function get_opt_atm_prfl. If the profile is create manually or with user functions, it must contain at least the same variables and column names. It may also include pressure in the attributes (attr(atm, "press") <- value), and in this case the parameter press can be omitted.

The initial direction of the photon package is sampled assuming equal sensibility of the sensor to all directions within its FOV (hemispherical in the case of irradiance).

At each iteration, a free optical path is sampled and the photon package is moved in metric scale. A check is made to verify if photon reached the surface or if it escaped the system. If it reached the surface, its intensity ("photon package weight") is added to the appropriate accumulator bin. The accumulator geometry can be: (1) 'annular', which only keep track of the annular bin in which the photon reached the surface; (2) 'sectorial', wich divide the annuli in sectors for azimuthal assymetry; or (3) 'grid', which will keep track of the x,y position. The 'annular' accumulator is appropriate for symmetrical conditions, as nadir view over a Lambertian surface. Note that regardless of extent of the accumulator, a final bin will be added between ext and Inf. The resolution of the accumulator will have an impact on the performance, as too high resolution, specially in 'grid' geometry will require more photons to be traced to reduce statistical sampling variation. Note that when using 'sectorial' geometry, angular resolution in azimuth is fixed to 1 degree and break points returned by the function refer only to the radius break points (angular break points will be 0:360 always).

The reduced performance with high resolution grid and the need for high resolution close to the target pixel while low resolution is sufficient for larger distances suggests two separate runs. Those grids can be combined later with the function fit_grid_dpsf to produce a Zernike polynomial fit to the density PSF. This model can then be used to generate a discrete PSF kernel at any sensor resolution.

To reduce the random sampling fluctuation in the sectorial and grid geometries, the grid is 'folded' on its axis of symmetry (sensor azimuth line), averaged and then 'unfolded'.

Note that the returned values are fractional contribution per accumulation bin (integrals), and so are not normalized per area. Conversion to per area contribution can be made with function an_psf.

Value

A list with elements 'bin_phtw' containing the fractional contribution of downward reflected photons per accumulation bin (depends on geom), 'bin_brks' containing the break points of the accumulation bins, 'bin_mid' with the mid point of the bins and 'metadata' containing the input parameters.

Examples

# Using the C version of the code:
data(us76)
atm <- rayleigh_od(us76, lambda = 550) %>%
       get_opt_atm_prfl(atm = us76, tau_aer = 0.5, H_aer = 2, w0_aer = 0.89, 
         tau_ray_z = ., a_mol_z = rep(0, nrow(us76)))
cdf_aer <- cdf_lut(continental_ph_6sv[, c(1, 8)])
psf_int <- mc_saf(atm     = atm, 
                  geom    = 'annular', 
                  res     = 0.03, 
                  ext     = 10, 
                  np      = 1E5, 
                  mnw     = 1E-6,
                  cdf_aer = cdf_aer)
par(mar = c(5, 6, 3, 2))
x <- psf_int$bin_brks
plot(x, cum_psf(psf_int), xaxs = "i",  yaxs = "i", xlab = "Radius (km)", 
ylab = expression(2*pi*integral("r'"*f("r'")*d*"r'", 0, r)), xlim = c(0, 10))

AlexCast/apsfs documentation built on Feb. 1, 2024, 9:48 p.m.