fit_cie_sky_model: Fit CIE sky model

View source: R/fit_cie_sky_model.R

fit_cie_sky_modelR Documentation

Fit CIE sky model


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


  custom_sky_coef = NULL,
  std_sky_no = NULL,
  general_sky_type = NULL,
  twilight = TRUE,
  rmse = FALSE,
  method = "BFGS"



SpatRaster. A normalized greyscale image. Typically, the blue channel extracted from a canopy photograph. Please see read_caim and normalize.


SpatRaster built with zenith_image.


SpatRaster built with azimuth_image.


The data.frame returned by extract_rl or a data.frame with same structure and names.


Numeric vector of length 1. Zenith digital number, see extract_rl for how to obtain it.


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 R package 'suncalc'.


Numeric vector of length five. Custom starting coefficients of the sky model. By default, they are drawn from standard skies.


Numeric vector. Standard sky number from Table 1 from \insertCiteLi2016;textualrcaiman.


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.


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.


Logical vector of length one. If it is TRUE, the criteria for selecting the best sky model is to choose the one with less root mean square error calculated by using sky_points as reference values. Otherwise, the criteria is to evaluate the whole hemisphere by calculating the product between the square ratio of r to the sky model and the fraction of pixels from this new layer that are above one or below zero, and selecting the sky model that produce the least value.


Optimization method to use. See optim.


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. However, because the CIE sky model is indeed the adjoin of two mathematical model, it is capable of handling any non-linearity since it is not a physical model with strict assumptions.

Ultimately, when the goal is to calculate the ratio of canopy to sky digital numbers, if the latter is accurately constructed, any non-linearity will be canceled. Please, see interpolate_sky_points for further considerations.

Nevertheless, the recommended input for this function is data pre-processed with the HSP software package \insertCiteLang2013rcaiman. Please, refer to write_sky_points for additional details about HSP.

The following code exemplifies how this package can be used to compare the manually-guided fitting provided by HSP 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") 
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

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

If you use this function in your research, please cite \insertCiteLang2010;textualrcaiman in addition to this package.


The result includes the following: (1) the output produced by mle2, (2) the 5 coefficients, (3) observed and predicted values, the sun coordinates –zenith and azimuth angle in degrees–, (4) the relative luminance image calculated for every pixel using the estimated coefficients and corresponding sun coordinates, (4) the digital number at the zenith, and (5) the description of the standard sky from which the initial coefficients were drawn. See \insertCiteLi2016;textualrcaiman to know more about these coefficients.



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()


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

bin <- ootb_obia(caim, z, a)
bin <- bin & mask_hs(z, 0, 80)

r <- gbc(caim$Blue*255)
g <- sky_grid_segmentation(z, a, 10)
sun_coord <- extract_sun_coord(r, z, a, bin, g)
sky_points <- extract_sky_points(r, bin, g)
rl <- extract_rl(r, z, a, sky_points)
model <- fit_cie_sky_model(r, z, a, rl$sky_points,
                           rl$zenith_dn, sun_coord,
                           rmse = TRUE,
                           general_sky_type = "Partly cloudy")
sky_cie <- model$relative_luminance * model$zenith_dn
sky_cie <- normalize(sky_cie, 0, 1, TRUE)

#to provide custom starting coefficient
path <- system.file("external", package = "rcaiman")
skies <- utils::read.csv(file.path(path, "15_CIE_standard_skies.csv"))
custom_sky_coef <- skies[9, 1:5] %>% as.numeric()
fit_cie_sky_model(r, z, a, rl$sky_points,
                 rl$zenith_dn, sun_coord,
                 rmse = TRUE,
                 custom_sky_coef = custom_sky_coef)

#restricted view canopy photo
path <- system.file("external/APC_0020.jpg", package = "rcaiman")
caim <- read_caim(path)
caim <- normalize(caim)
diameter <- calc_diameter(lens(), sqrt(nrow(caim)^2 + ncol(caim)^2)/2, 90)
z <- zenith_image(diameter, lens())
caim <- expand_noncircular(caim, z, c(ncol(caim)/2, nrow(caim)/2))
m <- !$Red)
a <- azimuth_image(z)
caim[!m] <- 0
z <- normalize(z, 0, 90) * 20 # a diagonal FOV of 40 degrees, a rough guess

bin <- ootb_obia(caim)

g <- sky_grid_segmentation(z, a, 5, sequential = TRUE)
col <- terra::unique(g) %>% nrow() %>% rainbow() %>% sample()
plot(g, col = col)
r <- gbc(caim$Blue*255)
sun_coord <- extract_sun_coord(r, z, a, bin, g)
sky_points <- extract_sky_points(r, bin, g)
rl <- extract_rl(r, z, a, sky_points)
model <- fit_cie_sky_model(r, z, a, rl$sky_points,
                           rl$zenith_dn, sun_coord, twilight = FALSE)
sky_cie <- model$relative_luminance * model$zenith_dn

## End(Not run)

rcaiman documentation built on Sept. 20, 2022, 1:05 a.m.