fit_cie_sky_model {rcaiman} | R Documentation |
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 |
z |
SpatRaster built with |
a |
SpatRaster built with |
sky_points |
The data.frame returned by |
zenith_dn |
Numeric vector of length one. Zenith digital number, see
|
sun_coord |
An object of class list. The result of a call to
|
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 |
rmse |
Logical vector of length one. If it is |
method |
Optimization method to use. See |
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)