Skip to content

Commit

Permalink
remove deleted insol dependency
Browse files Browse the repository at this point in the history
  • Loading branch information
fpirotti committed Dec 10, 2023
1 parent b4a9eae commit afbe7a6
Show file tree
Hide file tree
Showing 4 changed files with 126 additions and 11 deletions.
3 changes: 1 addition & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -34,8 +34,7 @@ Suggests:
testthat,
rmarkdown,
shinyWidgets,
DBI,
insol,
DBI,
comf,
leafem,
shinyjs,
Expand Down
117 changes: 117 additions & 0 deletions R/sunpos.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,117 @@

#' sunpos
#' @description Returns a matrix of azimuth and zenith angles of the sun given the
#' unit vectors from the observer to the direction of the sun.
#'
#' @param sunv coordinates x, y, z of the unit vector in the direction
#' of the sun.
#'
#' @author Javier G. Corripio
#'
#' @return A matrix of azimuth and zenith angles.
#' @export
#'
#' @examples
#' sunpos(matrix( (1:6)/6, nrow=2, ncol=3))
sunpos <- function (sunv)
{
if (nargs() < 1) {
cat("USAGE: sunpos(sunvector)\n 3D vector\n")
return()
}
azimuth = degrees(pi - atan2(sunv[, 1], sunv[, 2]))
zenith = degrees(acos(sunv[, 3]))
return(cbind(azimuth, zenith))
}

#' radians
#' @description Converts angles to radians from degrees
#'
#' @param degree angle in degrees
#'
#' @author Javier G. Corripio
#'
#' @return angle in radians.
#' @export
#'
#' @examples
#' sunpos(matrix( (1:6)/6, nrow=2, ncol=3))
radians <- function (degree)
{
radian = degree * (pi/180)
return(radian)
}


#' sunvector
#' @description Calculates a unit vector in the direction of the sun from
#' the observer position.
#'
#'
#' @param jd Julian Day and decimal fraction.
#' @param latitude Latitude of observer in degrees and decimal fraction.
#' @param longitude Longitude of observer in degrees and decimal fraction.
#' @param timezone Time zone, west is negative.
#'
#' @details To calculate the sunvector to the nearest hour, give the Julian Day
#' with a precission better than 1/24; to approximate it to the nearest minute
#' use decimal fractions smaller than 1/(24*60), and so on.
#'
#' @author Javier G. Corripio
#'
#' @return 3 column matrix with the x, y , z coordinates of the sun vector.
#' @export
#'
#' @examples
#' sunvector(JD(Sys.time()),51.4778,-0.0017,0)
sunvector <- function (jd, latitude, longitude, timezone)
{
if (nargs() < 4) {
cat("USAGE: sunvector(jd,latitude,longitude,timezone)\n values in jd, degrees, hours\n")
return()
}
omegar = hourangle(jd, longitude, timezone)
deltar = radians(declination(jd))
lambdar = radians(latitude)
svx = -sin(omegar) * cos(deltar)
svy = sin(lambdar) * cos(omegar) * cos(deltar) - cos(lambdar) *
sin(deltar)
svz = cos(lambdar) * cos(omegar) * cos(deltar) + sin(lambdar) *
sin(deltar)
return(cbind(svx, svy, svz))
}



#' radians
#' @description Computes Julian Day from dates as POSIXct object.
#'
#' @usage JD(x, inverse=FALSE)
#'
#' @param x POSIXct object.
#' @param inverse Logical. If false (default) returns the Julian Days
#' corresponding to given dates. If TRUE returns the date corresponding
#' to input Julian days
#'
#' @details Class "POSIXct" represents the (signed) number of seconds since the
#' beginning of 1970 (in the UTC timezone) as a numeric vector, and Julian Day is
#' the number of days since January 1, 4713 BCE at noon UTC, so the Julian Day
#' is calculated as numeric(POSIXct)+2440587.5 days.
#'
#' @author Javier G. Corripio
#'
#' @return Julian Day
#' @export
#'
#' @examples
#' JD(Sys.time())
JD <- function (x, inverse = FALSE)
{
if (inverse) {
return(as.POSIXct((x - 2440587.5) * 86400, origin = ISOdate(1970,
1, 1, 0, 0, 0), format = "%Y-%m-%d %H:%M:%S"))
}
else {
return(as.numeric(x)/86400 + 2440587.5)
}
}
14 changes: 7 additions & 7 deletions data-raw/app/app.R
Original file line number Diff line number Diff line change
Expand Up @@ -383,7 +383,7 @@ solarApp <- function() {
lsCond$clo<<-input$clo
shinyWidgets::updateAirDateInput(session = shiny::getDefaultReactiveDomain(),
inputId = "bins", value=Sys.time())
sunposition<<-insol::sunpos(insol::sunvector(insol::JD(Sys.time()), input$lat, input$long, 1))[1,]
sunposition<<-sunpos(sunvector(JD(Sys.time()), input$lat, input$long, 1))[1,]

shiny::updateNumericInput(inputId = "forceSunElev",value = round(sunposition[[2]]) )
shiny::updateNumericInput(inputId = "forceSunAngle",value = round(sunposition[[1]]) )
Expand Down Expand Up @@ -416,7 +416,7 @@ solarApp <- function() {

# browser()
req(input$bins)
dd<-insol::sunpos(insol::sunvector(insol::JD(input$bins), input$lat, input$long, 1))
dd<-sunpos(sunvector(JD(input$bins), input$lat, input$long, 1))

sunposition <<- dd[1,]
lsCond <<- getMeteoStation(input$bins )
Expand Down Expand Up @@ -589,15 +589,15 @@ solarApp <- function() {
length = 2,
by = '1 day')[2]

days = insol::JD(seq(as.POSIXct(as.Date(input$bins)),
days = JD(seq(as.POSIXct(as.Date(input$bins)),
as.POSIXct(as.Date(dayafter)), by = 'min'))
sh<-seq(as.POSIXct(as.Date(input$bins)), as.POSIXct(as.Date(dayafter)), by = 'hour')
sh.lab <- format(sh, "%H:00")
hours <- insol::JD(sh)
hours <- JD(sh)

sp <- sunposition
sps <- insol::sunpos(insol::sunvector(days, input$lat, input$long, 1))
sps.h <- insol::sunpos(insol::sunvector(hours, input$lat, input$long, 1))
sps <- sunpos(sunvector(days, input$lat, input$long, 1))
sps.h <- sunpos(sunvector(hours, input$lat, input$long, 1))
sps <- sps[ sps[, "zenith"]>0 & sps[, "zenith"] < 90, ]

sps.h.w <- sps.h[, "zenith"]>0 & sps.h[, "zenith"] < 90
Expand All @@ -612,7 +612,7 @@ solarApp <- function() {
radial.lim = seq(0, 90, 30),
radial.labels = c("", "30°", "60°", "90°"),
clockwise = TRUE,
label.pos= insol::radians(seq(0, 360, 45)) ,
label.pos= radians(seq(0, 360, 45)) ,
labels = c("N", "NE", "E", "SE", "S", "SW", "W", "NW", "N"),
rp.type = 's',
cex=0.1
Expand Down
3 changes: 1 addition & 2 deletions data-raw/prepareData.R
Original file line number Diff line number Diff line change
@@ -1,5 +1,4 @@
library(terra)
library(insol)
library(lidR)
library(rPET)
library(data.table)
Expand Down Expand Up @@ -63,7 +62,7 @@ makeplot <- function(datec="2022-08-03", cutoffZenith=80){
}
}

dd <- data.frame(insol::sunpos(insol::sunvector(insol::JD(ss), 46, 11.9701, 0)))
dd <- data.frame(sunpos(sunvector(JD(ss), 46, 11.9701, 0)))
dd$date <- ss
dd2<- data.table(dd[dd$zenith<cutoffZenith,])
dd3 <- dd2[as.integer(rownames(dd2))%%5==0, ]
Expand Down

0 comments on commit afbe7a6

Please sign in to comment.