fit_cie_sky_model: Fit CIE sky model
In GastonMauroDiaz/rcaiman: CAnopy IMage ANalysis

fit_cie_sky_model

R Documentation

Fit CIE sky model

Description

Use general-purpuse optimization to find the parameters of the CIE sky model that best fit data sampled from a canopy photograph.

Usage

fit_cie_sky_model(
  rr,
  sun_coord,
  custom_sky_coef = NULL,
  std_sky_no = NULL,
  general_sky_type = NULL,
  twilight = 60,
  method = "BFGS",
  loss = "MAE"
)

Arguments

`rr`	An object of class list. The output of `extract_rel_radiance()` or an object with same structure and names.
`sun_coord`	An object of class list. The output of `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 or a numeric matrix with five columns. 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`	Numeric vector of length one. Sun zenith angle (in degrees). If the sun zenith angle provided through the `sun_coord` argument is below this value, sun zenith angles from 85 to 96 degrees will be tested when selecting initial optimization parameters from the general sky types Clear or Partially Cloudy (specifically, for standard sky numbers 7 to 15). This adjustment is necessary because `extract_sun_coord()` can mistakenly identify the visible center of the solar corona as the solar disk. Since `extract_sun_coord()` cannot output a zenith angle below 90 degrees, setting this value to 90 is equivalent to disabling this step. Actually, `extract_sun_coord()` cannot output a value very close to 90, therefore the testing start at 85, civic twilight is from 90 to 96 degrees.
`method`	The method to be used. See ‘Details’. Can be abbreviated.
`loss`	Character vector of length one. Specifies the error metric to use. Options are `"MAE"` (Mean Absolute Error) or `"RMSE"` (Root Mean Squared Error). Defaults to `"MAE"`. `"MAE"` is more robust to inaccurate sky points and outliers, such as those found in cloudy skies.

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_minmax()
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
rr <- extract_rel_radiance(r, z, a, sky_points)
model <- fit_cie_sky_model(rr, 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_m <- cie_sky_image(z, a, sun_coord$zenith_azimuth, sky_coef)
cie_sky_m <- cie_sky_manual * manual_input$zenith_dn
plot(r/cie_sky_m)

Value

A list with the following components:

The output produced by stats::optim
The 5 coefficients of the CIE model
Observed values
Predicted values
The estimated digital number at the zenith
The sun coordinates
The method used for optimization
The starting parameters

Note

The point selection tool of ‘ImageJ’ software can be used to manually digitize points and create a CSV file from which to read coordinates (see Examples). After digitizing the points on the image, use the dropdown menu Analyze>Measure to open the Results window. To obtain the CSV file, use File>Save As...

The QGIS software can also be used to manually digitize points. In order to do that, drag and drop the image in an empty project, create an new vector layer, digitize points manually, save the editions, and close the project. To create the new vector layer go to the dropdown menu Layer>Create Layer>New Geopackage Layer...

Choose "point" in the Geometry type dropdown list and make sure the CRS is EPSG:7589. To be able to input the points, remember to click first on the Toogle Editing icon, and then on the Add Points Feature icon.

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

References

\insertAllCited

Examples

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

# Manual method following 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)

## Idem for QGIS
## These points were automatically generated with ootb_fit_cie_sky_model(),
## but same method apply to manually generated points.
path <- system.file("external/ootb_sky_sky_points.gpkg",
                    package = "rcaiman")
sky_points <- terra::vect(path)
sky_points <- terra::extract(caim, sky_points, cells = TRUE)
sky_points <- terra::rowColFromCell(caim, sky_points$cell) %>% as.data.frame()
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) #taken with click() after x11(), then hardcoded here
sun_coord <- zenith_azimuth_from_row_col(z, a, c(nrow(z) - xy[2],xy[1]))
points(sun_coord$row_col[2], nrow(caim) - sun_coord$row_col[1],
       col = 3, pch = 1)

rr <- extract_rel_radiance(caim$Blue, z, a, sky_points)

set.seed(7)
model <- fit_cie_sky_model(rr, sun_coord,
                           general_sky_type = "Clear",
                           twilight = 90,
                           method = "CG")
summary(model$mle2_output)
plot(model$obs, model$pred)
abline(0,1)
lm(model$pred~model$obs) %>% summary()

sky_cie <- cie_sky_image(z, a,
                                model$sun_coord$zenith_azimuth,
                                model$coef) * model$zenith_dn
plot(sky_cie)
plot(caim$Blue/sky_cie)

## End(Not run)

GastonMauroDiaz/rcaiman documentation built on April 14, 2025, 9:39 a.m.