#' Create \emph{skyscapeR.horizon} object from Az/Alt data
#'
#' This function creates a \emph{skyscapeR.horizon} object from measurements of
#' azimuth and altitude.
#' @param az Array of azimuth values
#' @param alt Array of altitude values.
#' @param alt.unc (Optional) Either a single value or an array of altitude
#' uncertainty.
#' @param loc Location, a vector containing the latitude, longitude and elevation of
#' the location, in this order.
#' @param name Name of site.
#' @param smooth Boolean to control whether to smooth horizon profile using rolling mean.
#' Defaults to FALSE
#' @param .scale Rolling mean window for smoothing. See \code{\link{createHWT}}
#' @seealso \code{\link{plot.skyscapeR.horizon}}
#' @import zoo
#' @export
#' @seealso \code{\link{createHWT}}, \code{\link{downloadHWT}}
#' @examples
#' # Create a skyscapeR.horizon from 5 measurements:
#' az <- c(0,90,180,270,360)
#' alt <- c(0,5,5,0,0)
#' hor <- createHor(az, alt, 0.1, c(40.1,-8), 'Test')
#' plot(hor)
createHor <- function(az, alt, alt.unc=0.5, loc, name='', smooth=F, .scale=1000) {
# return result
hor <- c()
hor$metadata$name <- name
if (length(loc)<3) { cat('No site elevation given. Assuming 0 metres above sea level.'); loc <- c(loc,0) }
if (length(loc)<2) { stop('Site latitude, longitude and elevation needed.') }
hor$metadata$georef <- loc; names(hor$metadata$georef) <- c('Lat','Lon','Elev'); dim(hor$metadata$georef) <- c(1,3)
if (length(alt.unc) > 1 & length(alt.unc) < length(az)) { stop('Altitude uncertainty must be either a single value or have the same length as "az"') }
if (length(alt.unc) == 1 ) { alt.unc <- rep(alt.unc, length(az)) }
az <- c(az-360, az, az+360); alt <- rep(alt, 3); alt.unc <- rep(alt.unc, 3)
ind <- which(duplicated(az)); if (length(ind)>0) { az <- az[-ind]; alt <- alt[-ind]; alt.unc <- alt.unc[-ind] }
xx <- seq(0-90, 360+90, 0.01)
alt <- approx(az, alt, xx)$y
alt.unc <- approx(az, alt.unc, xx)$y
hor$data <- data.frame(az = xx, alt = alt, alt.unc = alt.unc)
if (smooth) {
yy <- zoo::rollmean(alt, .scale, fill=0)
hor$data$alt <- yy
yy <- zoo::rollmean(alt.unc, .scale, fill=0)
hor$data$alt.unc <- yy
}
class(hor) <- "skyscapeR.horizon"
return(hor)
}
#' Exports a \emph{skyscapeR.horizon} object into \emph{Stellarium} format
#'
#' This function exports any \emph{skyscapeR.horizon} object into the landscape
#' format of \emph{Stellarium}, ready to be imported.
#' @param hor Horizon data in \emph{skyscapeR.horizon} format.
#' @param name Horizon name to be displayed in \emph{Stellarium}, if different
#' from one in \emph{skyscapeR.horizon} object.
#' @param author (Optional) Author, to be included in \emph{landscape.ini} file.
#' @param description (Optional) Description, to be included in \emph{landscape.ini} file.
#' @param ground_col Color of ground. Defaults to \emph{Stellarium}'s default.
#' @param hor_col Color of horizon line. Defaults to \emph{Stellarium}'s default.
#' @seealso \code{\link{createHor}}, \code{\link{downloadHWT}}, \code{\link{plot.skyscapeR.horizon}}
#' @references \href{https://stellarium.org/}{Stellarium: a free open source planetarium}
#' @export
#' @import utils
#' @examples
#' # Downloads horizon data from HeyWhatsThat and exports it into Stellarium:
#' \dontrun{
#' hor <- downloadHWT('HIFVTBGK')
#' exportHor(hor, name='Test', description='Test horizon export to Stellarium')
#' }
exportHor = function(hor, name, author="skyscapeR", description, ground_col, hor_col) {
if (!is(hor,'skyscapeR.horizon')) { stop('No skyscapeR.horizon object found.') }
if (missing(name)) { name = hor$name }
if (missing(description)) { description <- paste0("Horizon created using skyscapeR ", packageVersion('skyscapeR'), ".") }
if (substr(description,nchar(description),nchar(description)) != ".") { description <- paste0(description,". Horizon created using skyscapeR ", packageVersion('skyscapeR'), ".") }
if (missing(ground_col)) {ground_col = ".15,.45,.15"}
if (missing(hor_col)) {hor_col = ".75,.45,.45"}
# Horizon data
data <- data.frame(x = hor$data$az, y = hor$data$alt)
write.table(sprintf("%3.2f %4.2f", data$x, data$y), file="horizon.txt", quote=FALSE, sep = " ", row.names=F, col.names=F, fileEncoding = "UTF-8")
# Landscape.ini file
fileConn<-file("landscape.ini", encoding='UTF-8')
string.text <- c("[landscape]", paste0("name = ",name), "type = polygonal", paste0("author = ",author),
paste0("description = ", description), "polygonal_horizon_list = horizon.txt",
"polygonal_angle_rotatez = 0.00001", paste0("ground_color = ",ground_col),
paste0("horizon_line_color = ",hor_col), "",
"[location]", "planet = Earth", paste0("latitude = ",hor$metadata$georef[1],"d"),
paste0("longitude = ",hor$metadata$georef[2],"d"), paste0("altitude = ",hor$metadata$georef[3]))
writeLines(string.text, fileConn)
close(fileConn)
# Save as zip file
zip(zipfile=paste0(name,'-horizon.zip'), files=c("landscape.ini", "horizon.txt"))
file.remove(c('landscape.ini', 'horizon.txt'))
}
#' Download horizon data from \emph{HeyWhatsThat}
#'
#' This function downloads previously created horizon data
#' from \emph{HeyWhatsThat}, given its ID, and saves it as
#' a \emph{skyscapeR.horizon} object.
#' @param HWTID This is the 8 character ID attributed by
#' \emph{HeyWhatsThat.com}
#' @export
#' @import utils
#' @references \href{http://heywhatsthat.com/}{HeyWhatsThat.com}
#' @seealso \code{\link{createHWT}}
#' @examples
#' \dontrun{
#' # Retrieve horizon data for \href{https://www.heywhatsthat.com/?view=HIFVTBGK}{Liverpool Cathedral}:
#' hor <- downloadHWT('HIFVTBGK')
#' }
downloadHWT <- function(HWTID) {
if (nchar(HWTID) != 8) { stop('Incorrect HeyWhatsThat ID.') }
## Horizon metadata
test <- readLines(paste0("http://www.heywhatsthat.com/iphone/pan.cgi?id=",HWTID))
hor <- c()
# Lat/Lon/Elev
mypattern = '<div class=\"details_data\">([^<]*)</div>'
datalines = grep(mypattern,test,value=TRUE)
getexpr = function(s,g)substring(s,g,g+attr(g,'match.length')-1)
gg = gregexpr(mypattern,datalines)
matches = mapply(getexpr,datalines,gg)
result = gsub(mypattern,'\\1',matches)
names(result) = NULL
result[c(1,2,4)]
Lat <- as.numeric(strtrim(result[1],regexpr("°", result[1])[1]-1))
if (substr(result[1],nchar(result[1]),nchar(result[2])) == "S") { Lat <- -Lat }
Lon <- as.numeric(strtrim(result[2],regexpr("°", result[2])[1]-1))
if (substr(result[2],nchar(result[2]),nchar(result[2])) == "W") { Lon <- -Lon }
aux <- strtrim(result[4],regexpr(" ", result[4])[1]-1)
Elev <- as.numeric(substr(aux,2,nchar(aux)))
# Site Name
grep("pan_top_title", test)
mypattern = '<div id=\"pan_top_title\" class=\"ellipsis\" style=\"position: absolute; top: 46px; width: 296px\">([^<]*)</div>'
datalines = grep(mypattern,test,value=TRUE)
getexpr = function(s,g)substring(s,g,g+attr(g,'match.length')-1)
gg = gregexpr(mypattern,datalines)
matches = mapply(getexpr,datalines,gg)
result = gsub(mypattern,'\\1',matches)
names(result) = NULL
Name <- result
## Horizon data
hor.ex <- unique(substr(list.files(tempdir()),1,8)) # check if already downloaded
if (sum(HWTID == hor.ex) == 0) {
if (NROW(hor.ex) > 500) {
# delete oldest
details <- file.info(file.path(tempdir(),list.files(file.path(tempdir()))))
details <- details[with(details, order(as.POSIXct(ctime))), ]
files <- unique(substr(rownames(details),10,17))[1]
files <- paste0(file.path(tempdir()),"/",c(files, paste0(files,"-0-1")), ".png")
file.remove(files)
}
# download new one
curdir <- getwd()
setwd(tempdir())
download.file(paste0('http://www.heywhatsthat.com/api/horizon.csv?id=',HWTID,'&resolution=.125'), mode ='wb', destfile=paste0(HWTID,'.csv'), quiet=T)
setwd(curdir)
}
horizon <- read.csv(file.path(tempdir(), paste0(HWTID, '.csv')))
## Altitude Error
delta <- 9.73 # 9.73m for SRTM data
horizon$error <- rep(NA,NROW(data))
for (i in 1:NROW(horizon)) {
hor.alt <- horizon$altitude[i]
delta.elev <- horizon$distance..m.[i] * tan( hor.alt * pi/180)
aux0 <- atan( delta.elev / horizon$distance..m.[i] ) * 180/pi
aux1 <- atan( (delta.elev + delta) / horizon$distance..m.[i] ) * 180/pi
aux2 <- atan( (delta.elev - delta) / horizon$distance..m.[i] ) * 180/pi
horizon$error[i] <- mean(abs(c(aux1,aux2)-aux0))
}
# return result
hor$metadata <- c()
hor$metadata$ID <- HWTID
hor$metadata$name <- Name
hor$metadata$georef <- c(Lat, Lon, Elev); names(hor$metadata$georef) <- c('Lat','Lon', 'Elev'); dim(hor$metadata$georef) <- c(1,3)
hor$metadata$elevation <- Elev
hor$data <- data.frame(az = horizon$bin.bottom, alt = horizon$altitude, alt.unc = horizon$error)
class(hor) <- "skyscapeR.horizon"
return(hor)
}
#' Create and download horizon data from \emph{HeyWhatsThat}
#'
#' This function send a data request to \emph{HeyWhatsThat}, for the
#' creation of a horizon profile for a give Lat/Lon and elevation. It
#' then downloads the data and saves it as a \emph{skyscapeR.horizon}
#' object.
#' @param lat The latitude of the location.
#' @param lon The longitude of the location.
#' @param elevation (Optional) The elevation of the observer above
#' ground level in meters. Default is 1.6 meters (eye level).
#' @param name (Optional) Name for horizon.
#' @param src (Optional) Request source ID for \emph{HeyWhatsThat}. Default is
#' 'skyscapeR'. Only change this if you have been given a source ID by the
#' creator of \emph{HeyWhatsThat}.
#' @param verbose (Optional) Boolean switch to control output. Default is \emph{TRUE}.
#' @export
#' @import utils httr
#' @references \href{http://heywhatsthat.com/}{HeyWhatsThat.com}
#' @seealso \code{\link{downloadHWT}}
#' @examples
#' \dontrun{
#' # Create and retrieve horizon data for the London Mithraeum:
#' hor <- createHWT(lat=ten(51,30,45), lon=ten(0,5,26.1), name='London Mithraeum')
#' }
createHWT <- function(lat, lon, elevation=1.6, name, src='skyscapeR', verbose=T) {
if (verbose) { cat('Sending request to HeyWhatsThat servers...') }
test <- httr::GET(url=paste0('http://www.heywhatsthat.com/api/query?src=',src,'&lat=',lat,'&lon=',lon,'&elev_agl=',elevation))
if (test$headers$`x-app-status` == '200 OK') {
if (verbose) { cat('Done.\nWaiting for computation to finish...') }
} else { stop(test$headers$`x-app-status`) }
HWT.ID <- substr(httr::content(test),1,8)
Sys.sleep(30)
t <- 0
while (t==0) {
test <- httr::GET(url=paste0('http://www.heywhatsthat.com/api/ready?src=',src,'&id=',HWT.ID))
t <- as.numeric(substr(httr::content(test),1,1))
Sys.sleep(60)
}
if (verbose) { cat('Done.\nDownloading data...') }
hor <- downloadHWT(HWT.ID)
hor$metadata$name <- name
if (verbose) { cat('Done.\n') }
return(hor)
}
#' Retrieves horizon altitude for a given azimuth from a given horizon profile
#'
#' This function retrieves the horizon altitude for a given azimuth from
#' a previously created \emph{skyscapeR.horizon} object via spline interpolation.
#' @param hor A \emph{skyscapeR.horizon} object from which to retrieve horizon altitude.
#' @param az Array of azimuth(s) for which to retrieve horizon altitude(s).
#' @param return.unc (Optional) Boolean switch control where to output altitude uncertainty.
#' Default is \emph{FALSE}.
#' @export
#' @import stats
#' @seealso \code{\link{createHor}}, \code{\link{downloadHWT}}
#' @examples
#' hor <- downloadHWT('HIFVTBGK')
#' hor2alt(hor, 90)
hor2alt = function(hor, az, return.unc=F) {
hh <- splinefun(hor$data$az, hor$data$alt)
alt <- round(hh(az), 2)
if (return.unc) {
hh <- splinefun(hor$data$az, hor$data$alt.unc)
unc <- round(hh(az), 2)
alt[2] <- unc
}
return(alt)
}
#' Retrieves maximum declination for a given horizon
#'
#' @param hor Either a \emph{skyscapeR.horizon} object or the latitude of the place.
#' @param alt (Optional) Altitude of the horizon at northernmost point. Only necessary if no
#' \emph{skyscapeR.horizon} object is provided. Default is 0.
#' @param range (Optional) Range of azimuths to check for maximum declination. Default is 25
#' degrees either side of North.
#' @export
#' @examples
#' hor <- downloadHWT('HIFVTBGK')
#' hor2max.dec(hor)
hor2max.dec = function(hor, alt=0, range=25) {
xx <- seq(-range, range, 0.01)
if (is(hor,'skyscapeR.horizon')) { decs <- az2dec(xx, hor) } else { decs <- az2dec(xx, hor[1], alt) }
return(max(decs))
}
#' Retrieves minimum declination for a given horizon
#'
#' @param hor Either a \emph{skyscapeR.horizon} object or the latitude of the place.
#' @param alt (Optional) Altitude of the horizon at southernmost point. Only necessary if no
#' \emph{skyscapeR.horizon} object is provided. Default is 0.
#' @param range (Optional) Range of azimuths to check for maximum declination. Default is 25
#' degrees either side of South
#' @export
#' @examples
#' hor <- downloadHWT('HIFVTBGK')
#' hor2min.dec(hor)
hor2min.dec = function(hor, alt=0, range=25) {
xx <- seq(-range, range, 0.01)+180
if (is(hor,'skyscapeR.horizon')) { decs <- az2dec(xx, hor) } else { decs <- az2dec(xx, hor[1], alt) }
return(min(decs))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.