oscarperpinan / solar Goto Github PK

View Code? Open in Web Editor NEW

33.0 33.0 13.0 9.86 MB

Solar Radiation and Photovoltaic Systems with R

Home Page: http://oscarperpinan.github.io/solar/

License: GNU General Public License v3.0

R 100.00%

solar's People

Contributors

Stargazers

Watchers

Forkers

parrad cosuroca stanislawswierc aabercrombie0492 fxcebx rangeo frimane hossein-madadi drroad stiloenergy nlebang solarization

solar's Issues

calcSol not using correct UTC time stamps

Given

library(solaR)

test <- calcSol(48.27, fBTd(mode = "serie",
                            start = as.POSIXct("2015-01-01"),
                            end = as.POSIXct("2015-02-01")),
                sample = "15 min")

head(test@solI$Bo0, n=40)

the resulting test object has the first non-zero value of solI$Bo0 at 8 a.m. UTC. According to another source of sunrise time, the local time (UTC+1 in this case) of sunrise was at 7:45 a.m. So the series of solI$Bo0 seems to be off by one hour.

time(test@solI$Bo0[1])

gives

: [1] "2015-01-01 UTC"

and

Sys.getenv("TZ")

gives "UTC".

sessionInfo()

R version 3.2.0 (2015-04-16)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

locale:
[1] LC_COLLATE=English_United States.1252 
[2] LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] solaR_0.41          latticeExtra_0.6-26 RColorBrewer_1.1-2 
[4] lattice_0.20-33     zoo_1.7-12         

loaded via a namespace (and not attached):
[1] compiler_3.2.0 tools_3.2.0    grid_3.2.0

NA values in Gef

Fix fInclin so NA values in irradiation on the inclined plane are set to 0.

Define * operator

Define * operator and additional classes (generator, inverter) to operate (for example) G0 * generator * inverter

Import packages other than base

http://developer.r-project.org/blosxom.cgi/R-devel/NEWS/2015/06/29#n2015-06-29

Please check docu of "fSolD"

under "value", it says:
ws:
Sunset angle (in radians) for each day of year. Due to the convention which con-
siders that the solar hour angle is negative before midday, this angle is negative.

This is confusing: If the function really returns sunSET angle, and the sun sets AFTER midday, the angle should be POSITIVE, shouldn't it?

Problem with 'NA' values in 'Gef' when using 'calcGef' (package version 0.41)

Hi Oscar,
Thank you for creating this package! By its description it seems as if it can solve a lot of the problems I will need to deal with in my research on PV systems. I installed your package version 0.41 from CRAN. I wanted to use calcGef to calculate global irradiation on the generator plane from the values of intradaily global and diffuse irradiation on the horizontal plane. When I started using the package with my data I ran into some problems with NA values for Gef.

Could you help me understand if this is a problem in the package itself, or the way I use it?

Please notice that I do not have ambient temperature, but from what I know global and diffuse irradiation should be sufficient for this calculation. That's why I ignored the following warning:

In checkG0Ta(file, maxmin = TRUE) :
  Ambient temperature data is not available. A new column with a constant value
  has been added.

Here is the code I used:

library("solaR")
library("dplyr")

# Read solar data
channels <- read.csv("channels_solar_30_fixed.csv",
  col.names = c("time", "glob0", "indr0", "glob30", "indr30"))

# Convert radiation from [kW/m2] to [W/m2]
channels_solar[, -1] <- channels_solar[, -1] * 1000

# Make sure global radiation is greater or equal to indirect
channels_solar <- channels_solar %>%
  mutate(
    glob0 = max(glob0, indr0),
    glob30 = max(glob30, indr30))

# Set location for Dublin IE
lat <- 53.347778,
lon <- -6.259722

# Prepare zoo object with the right columns
meteo_zoo <- zoo(
  data.frame(
    G0 = channels_solar$glob0,
    D0 = channels_solar$indr0,
    B0 = channels_solar$glob0 - channels_solar$indr0),
  local2Solar(channels_solar$time, lon=lon))

gef <- calcGef(
  lat,
  dataRad = meteo_zoo,
  modeRad = 'bdI',
  corr = 'none',
  beta = 30,
  alfa = 0)

print(gef)
xyplot(gef)

gef_df <- as.data.frameI(gef)
head(gef_df, 48)

Here is the result for a single day which shows NA values.

                    Gef      Bef Def day month year
2008-12-31 23:34:57   0        0   0 366    12 2008
2009-01-01 00:04:57   0        0   0   1     1 2009
2009-01-01 00:34:57   0        0   0   1     1 2009
2009-01-01 01:04:57   0        0   0   1     1 2009
2009-01-01 01:34:57   0        0   0   1     1 2009
2009-01-01 02:04:57   0        0   0   1     1 2009
2009-01-01 02:34:57   0        0   0   1     1 2009
2009-01-01 03:04:57   0        0   0   1     1 2009
2009-01-01 03:34:57   0        0   0   1     1 2009
2009-01-01 04:04:57   0        0   0   1     1 2009
2009-01-01 04:34:57   0        0   0   1     1 2009
2009-01-01 05:04:57   0        0   0   1     1 2009
2009-01-01 05:34:57   0        0   0   1     1 2009
2009-01-01 06:04:57   0        0   0   1     1 2009
2009-01-01 06:34:57   0        0   0   1     1 2009
2009-01-01 07:04:57   0        0   0   1     1 2009
2009-01-01 07:34:57   0        0   0   1     1 2009
2009-01-01 08:04:57   0        0   0   1     1 2009
2009-01-01 08:34:57  NA 13465037  NA   1     1 2009
2009-01-01 09:04:57  NA  5412793  NA   1     1 2009
2009-01-01 09:34:57  NA  4026320  NA   1     1 2009
2009-01-01 10:04:57  NA  3446681  NA   1     1 2009
2009-01-01 10:34:57  NA  3167176  NA   1     1 2009
2009-01-01 11:04:57  NA  2988023  NA   1     1 2009
2009-01-01 11:34:57  NA  2898145  NA   1     1 2009
2009-01-01 12:04:57  NA  2850976  NA   1     1 2009
2009-01-01 12:34:57  NA  2913485  NA   1     1 2009
2009-01-01 13:04:57  NA  2985269  NA   1     1 2009
2009-01-01 13:34:57  NA  3153930  NA   1     1 2009
2009-01-01 14:04:57  NA  3496556  NA   1     1 2009
2009-01-01 14:34:57  NA  4140284  NA   1     1 2009
2009-01-01 15:04:57  NA  5636209  NA   1     1 2009
2009-01-01 15:34:57  NA 16391656  NA   1     1 2009
2009-01-01 16:04:57   0        0   0   1     1 2009
2009-01-01 16:34:57   0        0   0   1     1 2009
2009-01-01 17:04:57   0        0   0   1     1 2009
2009-01-01 17:34:57   0        0   0   1     1 2009
2009-01-01 18:04:57   0        0   0   1     1 2009
2009-01-01 18:34:57   0        0   0   1     1 2009
2009-01-01 19:04:57   0        0   0   1     1 2009
2009-01-01 19:34:57   0        0   0   1     1 2009
2009-01-01 20:04:57   0        0   0   1     1 2009
2009-01-01 20:34:57   0        0   0   1     1 2009
2009-01-01 21:04:57   0        0   0   1     1 2009
2009-01-01 21:34:57   0        0   0   1     1 2009
2009-01-01 22:04:57   0        0   0   1     1 2009
2009-01-01 22:34:57   0        0   0   1     1 2009
2009-01-01 23:04:57   0        0   0   1     1 2009

Kind regards,
Stan

Surface Orientation Factor issue (Azimuth)

First of all, you will allow me to give thanks because your work is so useful for my thesis!

My problem: Using SIAR data and prodGCPV for multiple orientations (beta, alfa) i get a nice contour plot of the Surface Orientation Factor. However, it is always moved 30 degrees (two hours) to west:

I attach my project (1MB) here: Controller (Main), data (SIAR), example plot and results in rData.

I suppose the error was in time data (solar, local, utc) but the hourly sumatory (for a year) is symetric and centered on 12:00.
I used several SIAR stations and years taking into account:

Siar. Descripción de datos publicados en la web: "Los datos se muestran en formato de hora UTC desde el día 1 de abril de 2014, para datos anteriores los datos se muestran en hora solar."
Conclusión: solaR trabaja con hora solar - UTC: Paso a solar time los datos posteriores de 1 Abril 2014.

Controller code:

library('data.table')
library('lubridate')
library('solaR')

##########################################################################
# constants & functions
##########################################################################

# Constantes del plot
C.PLOT = list(
  PALETTE = colorRampPalette(c("dodgerblue4","deepskyblue", "gold", "darkorange2", "firebrick"), space = "Lab"),
  LEVELS = c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.85,0.9,0.92,0.94,0.96,0.98,0.99,1),
  PALETTE.BW = colorRampPalette(c('#151515', '#F0F0F0'), space = "Lab")
)

readSiar = function(path, pattern, lon = NULL){
  # Lectura de datos
  files = list.files(path = path, pattern = pattern, full.names = T)
  message('Reading files...')
  print(files)
  csv= lapply(files, function(f){
    data.table(read.csv(f, header = T, sep = ';', dec = ',', fileEncoding = 'UTF-16LE'))})
  dt = rbindlist(csv)
  message('Data columns')
  print(names(dt))
  
  # Defino index
  dt[, time := as.POSIXct(paste0(Fecha, ' ', paste(HoraMin%/%100, HoraMin%%100, '00', sep = ':')), format = '%d/%m/%Y %H:%M:%S', tz = 'UTC')]
  # Fuerzo TZ
  attr(dt$time, 'tzone') = 'UTC'
  
  ##########################################################################
  # Siar - Datos publicados en web. [http://www.mapama.gob.es/es/desarrollo-rural/temas/gestion-sostenible-regadios/Datos_Publicados_tcm7-325744.pdf]
  # Los datos se muestran en formato de hora UTC desde el dÃa 1 de abril de 2014, para datos anteriores los datos se muestran en hora solar. 
  # ConclusiÃ³n: solaR trabaja con hora solar: Paso a  solar los datos posteriores al 1 Abril 2014
  dt[year(time) > 2014 | (year(time) == 2014 & month(time) > 3), time := local2Solar(time, lon = lon)]
  # dt[year(time) < 2014 | (year(time) == 2014 & month(time) < 4), time := solar2local(time, tz = 'UTC', lon = lon)]
  ##########################################################################

  # Subset solaR data: Temperatura en ÂºC y RadiaciÃ³n media en W/m2
  tempname = names(dt) %like% 'Temp'
  radname = names(dt) %like% 'Rad'
  res = dt[, .SD, .SDcols = c('time', names(dt)[names(dt) %like% 'Temp'], names(dt)[names(dt) %like% 'Rad'])]
  names(res)[2:3] = c('Ta', 'G0')
  
  res
}


##########################################################################
# inputs and read data (Siar)
##########################################################################

# Variables a analizar
vars = c('Eac', 'Edc', 'Yf', 'Gefd', 'Gd')

# Azimut (Âº). SÃ³lo coordenadas positivas (Oeste). Luego se aplica simÃ©tria
az.sim = c(2,4,6,8,10,30,45,60,75,90)
# InclinaciÃ³n (Âº)
alt = c(0,5,10,20,30,40,50,60,70,80,90)

# Datos EstaciÃ³n MeteorolÃ³gica
num = 17
prov = 46
stations = fread(file.path(getwd(), 'datos', 'stations.csv'), header = T, encoding = 'UTF-8')
station = stations[N_Estacion == num & N_Provincia == prov]
message(paste0('Station selected: ', station$Estacion))
print(station)
latitud = station$lat
longitud = station$lon

# Hora solar de acuerdo con la bibliografia
pica02 = readSiar(file.path(getwd(), 'datos'), 'V17_Picassent_01_01_2002_31_12_2002', lon = longitud)
# Hora UTC de acuerdo con la bibliografia
pica15 = readSiar(file.path(getwd(), 'datos'), 'V17_Picassent_01_01_2015_31_12_2015', lon = longitud)

# Otros datos
# pica02 = readSiar(file.path(getwd(), 'datos'), 'V107_Manises_01_01_2010_31_12_2010', lon = longitud)
# pica15 = readSiar(file.path(getwd(), 'datos'), 'V26_Betera_01_01_2015_31_12_2015', lon = longitud)

# Compruebo TZ y DST
attr(pica02$time, 'tzone')
attr(pica15$time, 'tzone')
sum(as.POSIXlt(pica02$time)$isdst)
sum(as.POSIXlt(pica15$time)$isdst)
print(pica02)
print(pica15)


##########################################################################
# solaR procedure for multiple alpha,beta
##########################################################################

# Constructor Meteo
meteo02 = dfI2Meteo(file = pica02, lat = latitud, format = format(pica02$time), time.col = 'time', source = 'Picassent')
meteo15 = dfI2Meteo(file = pica15, lat = latitud, format = format(pica15$time), time.col = 'time', source = 'Picassent')

# Lista donde se guarda los resultados de prodGCPV en cada orientaciÃ³n
li.prod02 = list()
li.prod15 = list()

if (length(az.sim) > 1) {
  az = c(sort(-1*az.sim), 0, az.sim)
} else az = az.sim

# # Mejora el rendimiento en el bucle pero no se asigna la clase de cada item
# length(li.prod02) = length(az)*length(alt)
# length(li.prod15) = length(az)*length(alt)

for(i in az) {
  for (j in alt){
    id = paste0('a', i, '-b', j)
    message(paste('alpha =', i, '; alt =', j))
    li.prod02[[id]] = prodGCPV(lat = latitud, modeTrk = 'fixed', 
                               modeRad = 'bdI', dataRad = meteo02, 
                               keep.night = T, sunGeometry = 'michalsky', 
                               corr = 'BRL', iS = 2, horizBright = F, 
                               HCPV = F, alb = 0.2, beta = j, alfa = i)
    li.prod15[[id]] = prodGCPV(lat = latitud, modeTrk = 'fixed', 
                               modeRad = 'bdI', dataRad = meteo15, 
                               keep.night = T, sunGeometry = 'michalsky', 
                               corr = 'BRL', iS = 2, horizBright = F, 
                               HCPV = F, alb = 0.2, beta = j, alfa = i)
  }
}


##########################################################################
# Anual results 
##########################################################################

# Resultados anuales para cada orientaciÃ³n 
dt02 = data.table()
dt15 = data.table()

for (item in li.prod02) {
  dt02 = rbind(dt02, data.table(Eac = sum(item@prody$Eac), Edc = sum(item@prody$Edc), 
                                Yf = sum(item@prody$Yf), Gefd = sum(item@Gefy$Gefd), 
                                Gd = sum(item@Gefy$Gd), alpha = item@angGen$alfa, beta = item@angGen$beta), 
               use.names = T, fill = T)
}
for (item in li.prod15) {
  dt15 = rbind(dt15, data.table(Eac = sum(item@prody$Eac), Edc = sum(item@prody$Edc), 
                                Yf = sum(item@prody$Yf), Gefd = sum(item@Gefy$Gefd), 
                                Gd = sum(item@Gefy$Gd), alpha = item@angGen$alfa, beta = item@angGen$beta), 
               use.names = T, fill = T)
}

##########################################################################
# Plot results 
##########################################################################

# Defino x, y para contour plot fuera del bucle
x = az
y = alt

#Recorro todas las variables que analizo
for (var in vars){
  
  # Resultados anuales para cada orientaciÃ³n ordenados en una matriz
  mat02 = matrix(0, nrow = length(az), ncol = length(alt))
  mat15 = matrix(0, nrow = length(az), ncol = length(alt))
  for(i in c(1:length(az))){
    for(j in c(1:length(alt))){
      mat02[i,j] = dt02[alpha == az[i] & beta == alt[j]][[var]]
      mat15[i,j] = dt15[alpha == az[i] & beta == alt[j]][[var]]
    }
  }
  # Plot: Se desviado al oeste !!!!
  mat = list(mat02, mat15)
  ye = c(2002,2015)
  for(i in 1:length(mat)) {
    z = mat[[i]]
    wh = which(z == max(z), arr.ind = T) + c(x[1], y[1])
    message(paste('OrientaciÃ³n Ã³ptima: alpha =', wh[1],', beta =', wh[2]))
    filled.contour(x,y,z,col=C.PLOT$PALETTE(17),
                   levels = C.PLOT$LEVELS*max(z),
                   key.axes = axis(3,NULL),
                   xlab = expression(paste(alpha," (Âº)")),
                   ylab = expression(paste(beta," (Âº)")),
                   plot.title = {main =paste('year =', ye[i], '; var =', var)},
                   plot.axes = {
                     axis(1, at = az)
                     axis(2, at = alt)
                     contour(x,y,z, level=C.PLOT$LEVELS*max(z), add=T, lwd=1,drawlabels = T,labels = C.PLOT$LEVELS,
                             labcex =1.1,col = grey(0,0.8))})
    
  }
}

##########################################################################
# Daily distribution of production
##########################################################################

# Pero sÃ que estÃ¡ centrada la producciÃ³n y la irradiaciÃ³n al mediodÃa solar, 
# incluso estÃ¡ desplazada a la maÃ±ana: Â¿quÃ© pasa?
# Nota: prod02[[1]] es de inclinaciÃ³n = beta = 0
test = as.data.table(li.prod02[[1]]@prodI)
test[, time := index(li.prod02[[1]]@prodI)]
test[, sum(Pac, na.rm=T), by = hour(time)]

save(dt02, dt15, file = 'res.RData')

Extra-atmospheric upper bound in fCompI is not documented

NA's in Gef at sunrise and sunset

Oscar, buenas tardes.
En primer lugar, felicitaciones y gracias por solaR, es una gran herramienta.
Soy Daniel Paniagua de Argentina, por lo que aclaro que escribiré en español para poder expresarme mejor, sabiendo que puedes entenderme. Si debo hacerlo en inglés por algún standard de Github por favor no dejes de avisarme.

Estoy teniendo dificultades para utilizar la función fInclin, ya que al utilizar G0 con datos cada una hora (intradaily), obtengo NA y NaN en Gef en los horarios cercanos al amanecer y al atardecer. Revisé los distintos temas abiertos (y cerrados) aquí y no conseguí resolverlo. Verfiqué principalemnte que G0 sea mayor a Bo0 en esos casos, además de otras cuestiones particulares.
Consideré reemplazarlos por cero, pero veo que en algunos casos el valor de radiación de G0 en esos horarios no es despreciable.

Dejo el código abajo y un archivo para importar G0 ("G0.txt", lo paso como txt para poder adjuntar, se debe cambiar extensión a csv).

Muchas gracias!
Saludos cordiales.
Daniel Paniagua

library(solaR)
library(dplyr)

#Importar archivo csv de G0
G0 <- read.csv("G0.csv")
G0 <- select(G0, date, G0, Ta)
G0$date <- as.POSIXlt(G0$date)

#Crear zoo de G0
zoo_G0 <- zoo(data.frame(G0=G0$G0, Ta=G0$Ta), order.by = G0$date)

#Crear objeto Meteo
G0_meteo <- zoo2Meteo(zoo_G0, lat=-32.7)

#Crear objeto G0
G0_G0 <- calcG0(lat=-32.7, modeRad='bdI', dataRad=G0_meteo)

#Crear objeto Sol
sol <- calcSol(lat=-32.7, BTi=index(zoo_G0))

#Cálculo de ángulos
theta <- fTheta(sol=sol, beta=30, alfa = 180, modeTrk = "fixed")

#Cálculo de radiación en el plano inclinado (Gef)
Gef <- fInclin(compI=G0_G0, angGen=theta)

#Pasa Gef a data.frame
df.Gef <- as.data.frame(Gef)

G0.txt

seq.default error

Dear Oscar,

I found this package short time ago and I'm trying to get data for my location.

I tried calcG0 with the values from the documentation and and mode aguilar. It was running without any problem.

Then I changed the G0dm values

G0dm=c(0.805,1.56,3.04,4.490,5.17,5.7,5.52,4.67,3.37,1.96,1.0,0.674)*1000

Then I got an error

> calcG0(lat=48.7, modeRad='aguiar',dataRad=G0dm) Error in seq.default(rng[1], rng[2], length = 11) : 'from' cannot be NA, NaN or infinite
Is this an error or do I not understand the logic of the library? I tested it with R 3.3.0 and 3.2.4 same result.

Regards
Midon

Use units package

https://r-quantities.github.io/units/

calcGef with beta=0 is not equal to g0

Why is the calculation of calcGef with beta=0 not equal to g0? Thanks for your help.

My test Code:

lat = 34.8363472
date = c("2017-12-14")
G0 = c(3600)
TempMax = c(30)
TempMin = c(0)
df = data.frame(date, G0, TempMax, TempMin)
bd=df2Meteo(df, dates.col='date', lat=lat, format='%Y-%m-%d')
g0 <- calcG0(lat=lat, modeRad='bd', dataRad=bd)
print(as.zooD(g0))

gef <- calcGef(lat = lat, modeRad = "prev", dataRad = g0, beta =0, alfa = 0)
print(as.zooD(gef))

solaR sometimes gives NaNs for Pdc and NREL input data

NREL data from http://rredc.nrel.gov/solar/old_data/nsrdb/1991-2010/ provides for global radiation G0 zero values for night hours. From the other hand the expression
fd <- D0/G0 is used in the fCompI() function. This in some cases finally produces NaN in output results for some (not all) hours where G0=0.

Problem disappears if to add to G0 very small constant

Use the units package

This: http://r-spatial.org//r/2016/08/16/units2.html

Adopt maptools::sun-methods?

Oscar,

Is it possible that this package could receive the sun-methods from maptools https://maptools.r-forge.r-project.org/reference/sun-methods.html? Some 6 CRAN packages use the methods, and as.maptools retires would see changes in results using for example solartime. Sebastian Luque, the author of the methods, would be happy to see them moved from maptools.

Roger