document()
library(microclima)
tme <- as.POSIXlt(c(1:740) * 24 * 3600, origin = "2015-01-15 00:00", tz = 'GMT')
hourlydata <- hourlyNCEP(ncepdata = NA, lat = 50, long = -5, tme)
dem <- get_dem(lat = 50, long = -5.2, resolution = 3)
rw <- .pointradwind(hourlydata, dem, lat = 50, long = -5.2, l = 1, x = 1)
lat=50
long = -5.2
x=1
l=1
albr = 0.15
zmin = 0
slope = NA
aspect = NA
horizon = NA
difani = TRUE
m <- is_raster(dem)
m[is.na(m)] <- zmin
m[m < zmin] <- zmin
dem <- if_raster(m, dem)
xy <- data.frame(x = long, y = lat)
coordinates(xy) = ~x + y
proj4string(xy) = "+init=epsg:4326"
xy <- as.data.frame(spTransform(xy, crs(dem)))
reso <- res(dem)[1]
wsc36 <- 0
wsc36atground <- 0
for (i in 0:35) {
wscr <- windcoef(dem, i*10, hgt = 2, res = reso)
wscr2 <- windcoef(dem, i*10, hgt = 0, res = reso)
wsc36[i + 1] <- extract(wscr, xy)
wsc36atground[i + 1] <- extract(wscr2, xy)
}
wshelt <- 0
wsheltatground <- 0
hourlydata$winddir <- hourlydata$winddir%%360
for (i in 1:length(hourlydata$winddir)) {
daz <- round(hourlydata$winddir[i] / 10, 0) + 1
daz[daz > 36] <- daz[daz > 36] - 36
wshelt[i] <- wsc36[daz]
wsheltatground[i] <- wsc36atground[daz]
}
if (class(slope) == "logical") {
slope <- terrain(dem, unit = 'degrees')
slope <- extract(slope, xy)
}
if (class(aspect) == "logical") {
aspect <- terrain(dem, opt = 'aspect', unit = 'degrees')
aspect <- extract(aspect, xy)
}
svf <- skyviewveg(dem, array(l, dim = dim(dem)[1:2]),
array(x, dim = dim(dem)[1:2]), res = reso)
fr <- canopy(array(l, dim = dim(dem)[1:2]), array(x, dim = dim(dem)[1:2]))
svf <- extract(svf, xy)
if (class(horizon) == "logical") {
ha36 <- 0
for (i in 0:35) {
har <- horizonangle(dem, i*10, reso)
ha36[i + 1] <- atan(extract(har, xy)) * (180/pi)
}
} else ha36 <- rep(horizon, 36)
tme <- as.POSIXlt(hourlydata$obs_time)
tme <- as.POSIXlt(tme + 0, tz = 'UTC')
ha <- 0
jd <- julday(tme$year + 1900, tme$mon + 1, tme$mday)
for (i in 1:length(tme)) {
saz <- solazi(tme$hour[i], lat, long, jd[i], merid = 0)
saz <- round(saz / 10, 0) + 1
saz <- ifelse(saz > 36, 1, saz)
ha[i] <- ha36[saz]
}
swrad <-.shortwave.ts(hourlydata$rad_dni, hourlydata$rad_dif, jd, tme$hour,
lat, long, slope, aspect, ha, svf, x, l, albr, merid = 0,
difani = difani)
dni=hourlydata$rad_dni
dif=hourlydata$rad_dif
localtime=tme$hour
slope
aspect
ha[1:25]
svv=svf
merid=0
dst=0
difani=T
sa <- solalt(locatime, lat, long, jd, merid)
document()
rw <- .pointradwind(hourlydata, dem, lat = 50, long = -5.2, l = 1, x = 1)
m <- is_raster(dem)
m[is.na(m)] <- zmin
m[m < zmin] <- zmin
dem <- if_raster(m, dem)
xy <- data.frame(x = long, y = lat)
coordinates(xy) = ~x + y
proj4string(xy) = "+init=epsg:4326"
xy <- as.data.frame(spTransform(xy, crs(dem)))
reso <- res(dem)[1]
wsc36 <- 0
wsc36atground <- 0
for (i in 0:35) {
wscr <- windcoef(dem, i*10, hgt = 2, res = reso)
wscr2 <- windcoef(dem, i*10, hgt = 0, res = reso)
wsc36[i + 1] <- extract(wscr, xy)
wsc36atground[i + 1] <- extract(wscr2, xy)
}
wshelt <- 0
wsheltatground <- 0
hourlydata$winddir <- hourlydata$winddir%%360
for (i in 1:length(hourlydata$winddir)) {
daz <- round(hourlydata$winddir[i] / 10, 0) + 1
daz[daz > 36] <- daz[daz > 36] - 36
wshelt[i] <- wsc36[daz]
wsheltatground[i] <- wsc36atground[daz]
}
if (class(slope) == "logical") {
slope <- terrain(dem, unit = 'degrees')
slope <- extract(slope, xy)
}
if (class(aspect) == "logical") {
aspect <- terrain(dem, opt = 'aspect', unit = 'degrees')
aspect <- extract(aspect, xy)
}
svf <- skyviewveg(dem, array(l, dim = dim(dem)[1:2]),
array(x, dim = dim(dem)[1:2]), res = reso)
fr <- canopy(array(l, dim = dim(dem)[1:2]), array(x, dim = dim(dem)[1:2]))
svf <- extract(svf, xy)
if (class(horizon) == "logical") {
ha36 <- 0
for (i in 0:35) {
har <- horizonangle(dem, i*10, reso)
ha36[i + 1] <- atan(extract(har, xy)) * (180/pi)
}
} else ha36 <- rep(horizon, 36)
tme <- as.POSIXlt(hourlydata$obs_time)
tme <- as.POSIXlt(tme + 0, tz = 'UTC')
ha <- 0
jd <- julday(tme$year + 1900, tme$mon + 1, tme$mday)
for (i in 1:length(tme)) {
saz <- solazi(tme$hour[i], lat, long, jd[i], merid = 0)
saz <- round(saz / 10, 0) + 1
saz <- ifelse(saz > 36, 1, saz)
ha[i] <- ha36[saz]
}
swrad <-.shortwave.ts(hourlydata$rad_dni, hourlydata$rad_dif, jd, tme$hour,
lat, long, slope, aspect, ha, svf, x, l, albr, merid = 0,
difani = difani)
sa <- solalt(localtime, lat, long, jd, merid)
alt <- sa * (pi / 180)
zen <- pi / 2 - alt
saz <- solazi(localtime, lat, long, jd, merid)
azi <- saz * (pi / 180)
sl <- slope * (pi / 180)
aspe <- aspect * (pi / 180)
si <- cos(zen) * cos(sl) + sin(zen) * sin(sl) * cos(azi - aspe)
si[sa < ha] <- 0
si[si < 0] <- 0
dirr <- si * dni
document()
sa <- solalt(localtime, lat, long, jd, merid)
alt <- sa * (pi / 180)
zen <- pi / 2 - alt
saz <- solazi(localtime, lat, long, jd, merid)
azi <- saz * (pi / 180)
sl <- slope * (pi / 180)
aspe <- aspect * (pi / 180)
si <- cos(zen) * cos(sl) + sin(zen) * sin(sl) * cos(azi - aspe)
si[sa < ha] <- 0
si[si < 0] <- 0
dirr <- si * dni
sa <- solalt(localtime, lat, long, jd, merid)
alt <- sa * (pi / 180)
zen <- pi / 2 - alt
saz <- solazi(localtime, lat, long, jd, merid)
azi <- saz * (pi / 180)
sl <- slope * (pi / 180)
aspe <- aspect * (pi / 180)
si <- cos(zen) * cos(sl) + sin(zen) * sin(sl) * cos(azi - aspe)
si[sa < ha] <- 0
si[si < 0] <- 0
dirr <- si * dni
rw <- .pointradwind(hourlydata, dem, lat = 50, long = -5.2, l = 1, x = 1)
m <- is_raster(dem)
m[is.na(m)] <- zmin
m[m < zmin] <- zmin
dem <- if_raster(m, dem)
xy <- data.frame(x = long, y = lat)
coordinates(xy) = ~x + y
proj4string(xy) = "+init=epsg:4326"
xy <- as.data.frame(spTransform(xy, crs(dem)))
reso <- res(dem)[1]
wsc36 <- 0
wsc36atground <- 0
for (i in 0:35) {
wscr <- windcoef(dem, i*10, hgt = 2, res = reso)
wscr2 <- windcoef(dem, i*10, hgt = 0, res = reso)
wsc36[i + 1] <- extract(wscr, xy)
wsc36atground[i + 1] <- extract(wscr2, xy)
}
wshelt <- 0
wsheltatground <- 0
hourlydata$winddir <- hourlydata$winddir%%360
for (i in 1:length(hourlydata$winddir)) {
daz <- round(hourlydata$winddir[i] / 10, 0) + 1
daz[daz > 36] <- daz[daz > 36] - 36
wshelt[i] <- wsc36[daz]
wsheltatground[i] <- wsc36atground[daz]
}
if (class(slope) == "logical") {
slope <- terrain(dem, unit = 'degrees')
slope <- extract(slope, xy)
}
if (class(aspect) == "logical") {
aspect <- terrain(dem, opt = 'aspect', unit = 'degrees')
aspect <- extract(aspect, xy)
}
svf <- skyviewveg(dem, array(l, dim = dim(dem)[1:2]),
array(x, dim = dim(dem)[1:2]), res = reso)
fr <- canopy(array(l, dim = dim(dem)[1:2]), array(x, dim = dim(dem)[1:2]))
svf <- extract(svf, xy)
if (class(horizon) == "logical") {
ha36 <- 0
for (i in 0:35) {
har <- horizonangle(dem, i*10, reso)
ha36[i + 1] <- atan(extract(har, xy)) * (180/pi)
}
} else ha36 <- rep(horizon, 36)
tme <- as.POSIXlt(hourlydata$obs_time)
tme <- as.POSIXlt(tme + 0, tz = 'UTC')
ha <- 0
jd <- julday(tme$year + 1900, tme$mon + 1, tme$mday)
for (i in 1:length(tme)) {
saz <- solazi(tme$hour[i], lat, long, jd[i], merid = 0)
saz <- round(saz / 10, 0) + 1
saz <- ifelse(saz > 36, 1, saz)
ha[i] <- ha36[saz]
}
swrad <-.shortwave.ts(hourlydata$rad_dni, hourlydata$rad_dif, jd, tme$hour,
lat, long, slope, aspect, ha, svf, x, l, albr, merid = 0,
difani = difani)
fr <- as.vector(canopy(array(l, dim = c(1,1)), array(x, dim = c(1,1))))
sa <- solalt(tme$hour, lat, long, jd, merid = 0)
kk <- ((x^2 + 1/(tan(sa * (pi/180))^2))^0.5)/(x + 1.774 * (x + 1.182)^(-0.733))
trd <- exp(-kk * l)
trf <- (1 - fr)
fgd <- fd * trd
document()
m <- is_raster(dem)
m[is.na(m)] <- zmin
m[m < zmin] <- zmin
dem <- if_raster(m, dem)
xy <- data.frame(x = long, y = lat)
coordinates(xy) = ~x + y
proj4string(xy) = "+init=epsg:4326"
xy <- as.data.frame(spTransform(xy, crs(dem)))
reso <- res(dem)[1]
wsc36 <- 0
wsc36atground <- 0
for (i in 0:35) {
wscr <- windcoef(dem, i*10, hgt = 2, res = reso)
wscr2 <- windcoef(dem, i*10, hgt = 0, res = reso)
wsc36[i + 1] <- extract(wscr, xy)
wsc36atground[i + 1] <- extract(wscr2, xy)
}
wshelt <- 0
wsheltatground <- 0
hourlydata$winddir <- hourlydata$winddir%%360
for (i in 1:length(hourlydata$winddir)) {
daz <- round(hourlydata$winddir[i] / 10, 0) + 1
daz[daz > 36] <- daz[daz > 36] - 36
wshelt[i] <- wsc36[daz]
wsheltatground[i] <- wsc36atground[daz]
}
if (class(slope) == "logical") {
slope <- terrain(dem, unit = 'degrees')
slope <- extract(slope, xy)
}
if (class(aspect) == "logical") {
aspect <- terrain(dem, opt = 'aspect', unit = 'degrees')
aspect <- extract(aspect, xy)
}
svf <- skyviewveg(dem, array(l, dim = dim(dem)[1:2]),
array(x, dim = dim(dem)[1:2]), res = reso)
fr <- canopy(array(l, dim = dim(dem)[1:2]), array(x, dim = dim(dem)[1:2]))
svf <- extract(svf, xy)
if (class(horizon) == "logical") {
ha36 <- 0
for (i in 0:35) {
har <- horizonangle(dem, i*10, reso)
ha36[i + 1] <- atan(extract(har, xy)) * (180/pi)
}
} else ha36 <- rep(horizon, 36)
tme <- as.POSIXlt(hourlydata$obs_time)
tme <- as.POSIXlt(tme + 0, tz = 'UTC')
ha <- 0
jd <- julday(tme$year + 1900, tme$mon + 1, tme$mday)
for (i in 1:length(tme)) {
saz <- solazi(tme$hour[i], lat, long, jd[i], merid = 0)
saz <- round(saz / 10, 0) + 1
saz <- ifelse(saz > 36, 1, saz)
ha[i] <- ha36[saz]
}
sw <-.shortwave.ts(hourlydata$rad_dni, hourlydata$rad_dif, jd, tme$hour,
lat, long, slope, aspect, ha, svf, x, l, albr, merid = 0,
difani = difani)
hourlyrad <- data.frame(swrad = sw$fgc, skyviewfact = svf, canopyfact = sw$cfc,
whselt = wsheltatground, windspeed = wshelt *  windsp,
slope = slope, aspect = aspect)
windsp <- windheight(hourlydata$windspeed, 10, 1)
hourlyrad <- data.frame(swrad = sw$fgc, skyviewfact = svf, canopyfact = sw$cfc,
whselt = wsheltatground, windspeed = wshelt *  windsp,
slope = slope, aspect = aspect)
xxx<-sw$fgc
length(xxx)
sa <- solalt(localtime, lat, long, jd, merid)
alt <- sa * (pi / 180)
zen <- pi / 2 - alt
saz <- solazi(localtime, lat, long, jd, merid)
azi <- saz * (pi / 180)
sl <- slope * (pi / 180)
aspe <- aspect * (pi / 180)
si <- cos(zen) * cos(sl) + sin(zen) * sin(sl) * cos(azi - aspe)
si[sa < ha] <- 0
si[si < 0] <- 0
dirr <- si * dni
if (difani) {
k <- dni / 4.87
} else k <- 0
isor <- 0.5 * dif * (1 + cos(sl)) * (1 - k)
cisr <- k * dif * si
sdi <- (slope + ha) * (pi / 180)
refr <- 0.5 * albr * (1 - cos(sdi)) * dif
fd <- dirr + cisr
fdf <- isor + refr
kk <- ((x ^ 2 + 1 / (tan(sa * (pi / 180)) ^ 2)) ^ 0.5) /
(x + 1.774 * (x + 1.182) ^ (-0.733))
trd <- exp(-kk * l)
fr <- as.vector(canopy(array(l, dim = c(1,1)), array(x, dim = c(1,1))))
trf <- (1 - fr)
fgd <- fd * trd * (1 - albr)
fged <- fdf * trf * (1 - albr) * svv
fgc <- fgd + fged
length(fgc)
fgc <- fgd + fged
cfc <- ((1 - trd) * fd + fr * fdf) / (fd + fdf)
cfc[is.na(cfc)] <- ((1 - trd[is.na(cfc)]) * 0.5 + fr * 0.5) / (0.5 + 0.5)
xxx=list(fgc, cfc)
length(xxx$fgc)
head(xxx)
document()
m <- is_raster(dem)
m[is.na(m)] <- zmin
m[m < zmin] <- zmin
dem <- if_raster(m, dem)
xy <- data.frame(x = long, y = lat)
coordinates(xy) = ~x + y
proj4string(xy) = "+init=epsg:4326"
xy <- as.data.frame(spTransform(xy, crs(dem)))
reso <- res(dem)[1]
wsc36 <- 0
wsc36atground <- 0
for (i in 0:35) {
wscr <- windcoef(dem, i*10, hgt = 2, res = reso)
wscr2 <- windcoef(dem, i*10, hgt = 0, res = reso)
wsc36[i + 1] <- extract(wscr, xy)
wsc36atground[i + 1] <- extract(wscr2, xy)
}
wshelt <- 0
wsheltatground <- 0
hourlydata$winddir <- hourlydata$winddir%%360
for (i in 1:length(hourlydata$winddir)) {
daz <- round(hourlydata$winddir[i] / 10, 0) + 1
daz[daz > 36] <- daz[daz > 36] - 36
wshelt[i] <- wsc36[daz]
wsheltatground[i] <- wsc36atground[daz]
}
if (class(slope) == "logical") {
slope <- terrain(dem, unit = 'degrees')
slope <- extract(slope, xy)
}
if (class(aspect) == "logical") {
aspect <- terrain(dem, opt = 'aspect', unit = 'degrees')
aspect <- extract(aspect, xy)
}
svf <- skyviewveg(dem, array(l, dim = dim(dem)[1:2]),
array(x, dim = dim(dem)[1:2]), res = reso)
fr <- canopy(array(l, dim = dim(dem)[1:2]), array(x, dim = dim(dem)[1:2]))
svf <- extract(svf, xy)
if (class(horizon) == "logical") {
ha36 <- 0
for (i in 0:35) {
har <- horizonangle(dem, i*10, reso)
ha36[i + 1] <- atan(extract(har, xy)) * (180/pi)
}
} else ha36 <- rep(horizon, 36)
tme <- as.POSIXlt(hourlydata$obs_time)
tme <- as.POSIXlt(tme + 0, tz = 'UTC')
ha <- 0
jd <- julday(tme$year + 1900, tme$mon + 1, tme$mday)
for (i in 1:length(tme)) {
saz <- solazi(tme$hour[i], lat, long, jd[i], merid = 0)
saz <- round(saz / 10, 0) + 1
saz <- ifelse(saz > 36, 1, saz)
ha[i] <- ha36[saz]
}
sw <-.shortwave.ts(hourlydata$rad_dni, hourlydata$rad_dif, jd, tme$hour,
lat, long, slope, aspect, ha, svf, x, l, albr, merid = 0,
difani = difani)
windsp <- windheight(hourlydata$windspeed, 10, 1)
hourlyrad <- data.frame(swrad = sw$fgc, skyviewfact = svf, canopyfact = sw$cfc,
whselt = wsheltatground, windspeed = wshelt *  windsp,
slope = slope, aspect = aspect)
head(hourlyradwind, 20)
head(hourlyrad, 20)
library(RNCEP)
tmeout <- tme
tme <- .tme.sort(tme)
long <- ifelse(long > 180, long - 360, long)
int <- as.numeric(tme[2]) - as.numeric(tme[1])
lgth <- (length(tme) * int) / (24 * 3600)
tme <- as.POSIXlt(c(0:(lgth - 1)) * 3600 * 24, origin = min(tme), tz = 'UTC')
yrs <- unique(tme$year + 1900)
mths <- unique(tme$mon + 1)
ll <- data.frame(x = long, y = lat)
pre <- NCEP.gather('prate.sfc', level = 'gaussian',
years.minmax = c(min(yrs),max(yrs)),
months.minmax = c(min(mths):max(mths)),
lat.southnorth = c(ll$y,ll$y), lon.westeast = c(ll$x,ll$x),
return.units = FALSE, status.bar = FALSE, reanalysis2 = reanalysis2)
reanalysis2 = TRUE
head(tme)
library(devtools)
document()
devtools::document()
document()
devtools::document()
devtools::document()
demworld
library(microclima)
plot(demworld)
library(raster)
plot(demworld)
devtools::document()
library(raster)
# Get global dem
dem <- raster("C:/NOAA_NCEP/demtiles/dem.tif")
# NCEP layer reanalyses 1
tmean <- raster("C:/NicheMapR/NOAAdem/r1_air.2m.gauss.2018.nc", band = 1)
e1 <- extent(c(-0.9375,180,-89.49406, 89.49406))
e2 <- extent(c(180,360,-89.49406, 89.49406))
tmean1 <- crop(tmean, e1)
tmean2 <- crop(tmean, e2)
tmean2 <- shift(tmean2, x = -360)
tmean <- mosaic(tmean1, tmean2, fun = mean)
demworld <- resample(dem, tmean)
demworld[is.na(demworld)] <- 0
devtools::use_data(demworld, overwrite=T)
devtools::document()
demworld
library(devtools)
document()
document()
document()
document()
devtools::document()
document()
devtools::document()
?hourlytemp
devtools::document()
?hourlytemp
hist(microvars$dni[1:48])
document()
devtools::document()
devtools::document()
devtools::document()
library(devtools)
install_packages('backports')
install_package('backports')
install.packages('backports')
library(devtools)
document()
document()
library(devtools)
document()
library(devtools)
document()
?hourlyNCEP
tme <- as.POSIXlt(c(0:30) * 24 * 3600, origin ="2015-01-15 00:00", tz = "UTC")
# NB takes a while to download data
hdata<- hourlyNCEP(NA, 50, -5, tme)
head(hdata)
library(devtools)
document()
document()
devtools::document()
devtools::document()
document()
devtools::document()
?runauto_ncep
library(microclima)
?runauto.ncep
library(NicheMapR)
library(raster)
require(NicheMapR)
# Get dem for Pico, Azores
r <- get_dem(lat = 38.467429, long = -28.398995, resolution = 30)
plot(r)
# Takes ~ c. 15 minutes to run
temps <- runauto.ncep(r, "10/06/2010", "15/06/2010", hgt = 0.1, l = NA, x = NA,
habitat = "Barren or sparsely vegetated")
library(devtools)
document()
library(devtools)
document()
