fit_cie_sky_model: Fit CIE sky model

In rcaiman: CAnopy IMage ANalysis

Fit CIE sky model

Description

Use maximum likelihood to estimate the coefficients of the CIE sky model that best fit to data sampled from a canopy photograph.

Usage

fit_cie_sky_model(
  r,
  z,
  a,
  sky_points,
  zenith_dn,
  sun_coord,
  custom_sky_coef = NULL,
  std_sky_no = NULL,
  general_sky_type = NULL,
  twilight = TRUE,
  rmse = FALSE,
  method = "BFGS"
)

Arguments

r	SpatRaster. A normalized greyscale image. Typically, the blue channel extracted from a canopy photograph. Please see read_caim() and normalize().
z	SpatRaster built with zenith_image().
a	SpatRaster built with azimuth_image().
sky_points	The data.frame returned by extract_rl() or a data.frame with same structure and names.
zenith_dn	Numeric vector of length one. Zenith digital number, see extract_rl() for how to obtain it.
sun_coord	An object of class list. The result of a call to extract_sun_coord() or an object with same structure and names. See also row_col_from_zenith_azimuth() in case you want to provide values based on date and time of acquisition and the suncalc package.
custom_sky_coef	Numeric vector of length five. Custom starting coefficients of the sky model. By default, they are drawn from standard skies.
std_sky_no	Numeric vector. Standard sky number from \insertCiteLi2016;textualrcaiman's Table 1.
general_sky_type	Character vector of length one. It could be any of these: "Overcast", "Clear", or "Partly cloudy". See Table 1 from \insertCiteLi2016;textualrcaiman for additional details.
twilight	Logical vector of length one. If it is TRUE and the initial standard parameters belong to the "Clear" general sky type, sun zenith angles from 90 to 96 degrees will be tested (civic twilight). This is necessary since extract_sun_coord() would mistakenly recognize the center of what can be seen of the solar corona as the solar disk.
rmse	Logical vector of length one. If it is TRUE, the criteria for selecting the best sky model is to choose the one with the lowest root mean square error (RMSE) calculated by using the sky_points argument as the source of reference values. Otherwise, the criteria is to evaluate the whole image by calculating the out-of-range index as \sum_{i = 1}^{N}(r_i/sky_i)^2, where r is the r argument, sky is the raster obtained from the fitted model with cie_sky_model_raster() and zenith_dn, i is the index that represents the position of a given pixel on the raster grid, and N is the total number of pixels that satisfy either of these inequalities: r_i/sky_i<0 and r_i/sky_i>1.
method	Optimization method to use. See optim.

Details

This function is based on \insertCiteLang2010;textualrcaiman. In theory, the best result would be obtained with data showing a linear relation between digital numbers and the amount of light reaching the sensor. See extract_radiometry() and read_caim_raw() for further details. As a compromise solution, gbc() can be used.

The following code exemplifies how this package can be used to compare the manually-guided fitting provided by HSP \insertCiteLang2013rcaiman against the automatic fitting provided by this package. The code assumes that the user is working within an RStudio project located in the HSP project folder.

r <- read_caim("manipulate/IMG_1013.pgm") %>% normalize()
z <- zenith_image(ncol(r), lens())
a <- azimuth_image(z)
manual_input <- read_manual_input(".", "IMG_1013" )
sun_coord <- manual_input$sun_coord$row_col
sun_coord <- zenith_azimuth_from_row_col(z, sun_coord, lens())
sky_points <- manual_input$sky_points
rl <- extract_rl(r, z, a, sky_points)
model <- fit_cie_sky_model(r, z, a, rl$sky_points, rl$zenith_dn, sun_coord)
cie_sky <- model$relative_luminance * model$zenith_dn
plot(r/cie_sky)

r <- read_caim("manipulate/IMG_1013.pgm")
sky_coef <- read_opt_sky_coef(".", "IMG_1013")
cie_sky_manual <- cie_sky_model_raster(z, a, sun_coord$zenith_azimuth, sky_coef)
cie_sky_manual <- cie_sky_manual * manual_input$zenith_dn
plot(r/cie_sky_manual)

Value

object from the class list. The result includes the following: (1) the output produced by bbmle::mle2(), (2) the 5 coefficients, (3 and 4) observed and predicted values, (5) the digital number at the zenith, (6) the sun coordinates –zenith and azimuth angle in degrees–, and (7) the description of the standard sky from which the initial coefficients were drawn. See \insertCiteLi2016;textualrcaiman to know more about these coefficients.

Note

If you use this function in your research, please cite \insertCiteLang2010;textualrcaiman in addition to this package (⁠citation("rcaiman"⁠).

References

\insertAllCited

See Also

Other Sky Reconstruction Functions: cie_sky_model_raster(), fit_coneshaped_model(), fit_trend_surface(), fix_reconstructed_sky(), interpolate_sky_points(), ootb_sky_reconstruction()

Examples

## Not run: 
caim <- read_caim() %>% normalize()
z <- zenith_image(ncol(caim), lens())
a <- azimuth_image(z)

# Manual method after Lang et al. (2010)
# ImageJ can be used to digitize points
path <- system.file("external/sky_points.csv",
                    package = "rcaiman")
sky_points <- read.csv(path)
sky_points <- sky_points[c("Y", "X")]
colnames(sky_points) <- c("row", "col")
head(sky_points)
plot(caim$Blue)
points(sky_points$col, nrow(caim) - sky_points$row, col = 2, pch = 10)

xy <- c(210, 451) #originally captured with click() after x11()
sun_coord <- zenith_azimuth_from_row_col(z, z, c(nrow(z) - xy[2],xy[1]))
points(sun_coord$row_col[2], nrow(caim) - sun_coord$row_col[1],
       col = 3, pch = 1)

rl <- extract_rl(caim$Blue, z, a, sky_points)

set.seed(7)
model <- fit_cie_sky_model(caim$Blue, z, a, rl$sky_points,
                           rl$zenith_dn, sun_coord,
                           general_sky_type = "Clear",
                           rmse = FALSE,
                           twilight = FALSE,
                           method = "SANN")
summary(model$mle2_output)
plot(model$obs, model$pred)
abline(0,1)
r2 <- lm(model$pred~model$obs) %>% summary(.) %>% .$r.squared
r2
sky_cie <- cie_sky_model_raster(z, a,
                                model$sun_coord$zenith_azimuth,
                                model$coef) * model$zenith_dn
plot(sky_cie)
plot(caim$Blue/sky_cie)

# A quick demonstration of how to use interpolation to improve sky modelling
# after Lang et al. (2010)
sky <- interpolate_sky_points(rl$sky_points, caim$Blue, rmax = ncol(caim)/7)
plot(sky)
sky <- sky * rl$zenith_dn * (1 - r2) + sky_cie * r2
sky <- terra::cover(sky, sky_cie)
plot(sky)
plot(caim$Blue/sky)

# how to provide a custom starting coefficient
model <- fit_cie_sky_model(caim$Blue, z, a, rl$sky_points,
                           rl$zenith_dn, sun_coord,
                           custom_sky_coef = model$coef,
                           method = "SANN")
plot(model$obs, model$pred, ylim = range(model$obs))
abline(0,1)

## End(Not run)

rcaiman documentation built on Nov. 15, 2023, 1:08 a.m.