fit_cie_sky_model {rcaiman}

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 Li et al. (2016)'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 Li et al. (2016) 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 Lang et al. (2010). 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 (Lang et al. 2013) 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 Li et al. (2016) to know more about these coefficients.

Note

If you use this function in your research, please cite Lang et al. (2010) in addition to this package (⁠citation("rcaiman"⁠).

References

Lang M, Kodar A, Arumäe T (2013). “Restoration of above canopy reference hemispherical image from below canopy measurements for plant area index estimation in forests.” Forestry Studies, 59(1), 13–27. doi:10.2478/fsmu-2013-0008.

Lang M, Kuusk A, M~ottus M, Rautiainen M, Nilson T (2010). “Canopy gap fraction estimation from digital hemispherical images using sky radiance models and a linear conversion method.” Agricultural and Forest Meteorology, 150(1), 20–29. doi:10.1016/j.agrformet.2009.08.001.

Li DH, Lou S, Lam JC, Wu RH (2016). “Determining solar irradiance on inclined planes from classified CIE (International Commission on Illumination) standard skies.” Energy, 101, 462–470. doi:10.1016/j.energy.2016.02.054.

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)