oscarperpinan / solar

Solar Radiation and Photovoltaic Systems with R
http://oscarperpinan.github.io/solar/
GNU General Public License v3.0
35 stars 13 forks source link

Surface Orientation Factor issue (Azimuth) #12

Closed JoanViana closed 6 years ago

JoanViana commented 7 years ago

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: rplot 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:

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')
oscarperpinan commented 7 years ago

Hi,

I have reviewed your code and the data, but unfortunately I cannot find the cause of the problem. I suspect there is a problem with the time index of the data, but nowadays I have no time to prove it.

I have prepared code (extracted from help(calcG0)) with a different data to show that solaR is performing correctly. You may want to adapt it to your data in order to find the problem.

##NREL-MIDC
##La Ola, Lanai
##Latitude: 20.76685o North
##Longitude: 156.92291o West
##Elevation: 381 meters AMSL
##Time Zone: -10.0
latitud <- 20.77
longitud <- -156.94

NRELurl <- 'http://goo.gl/fFEBN'

dat <- read.table(NRELurl, header=TRUE, sep=',')
names(dat) <- c('date', 'hour', 'G0', 'B', 'D0', 'Ta')

##B is direct normal. We need direct horizontal.
dat$B0 <- dat$G0-dat$D0

##http://www.nrel.gov/midc/la_ola_lanai/instruments.html:
##The datalogger program runs using Greenwich Mean Time (GMT),
##data is converted to Hawaiin Standard Time (HST) after data collection
idxLocal <- with(dat, as.POSIXct(paste(date, hour), format='%m/%d/%Y %H:%M', tz='HST'))
idx <- local2Solar(idxLocal, lon=-156.9339)

NRELMeteo <- zoo(dat[,c('G0', 'D0', 'B0', 'Ta')], idx)

## Alpha and Beta combinations
M <- expand.grid(alfa = seq(-75, 75, 25),
                 beta = seq(0, 60, 15))

## Auxiliary function, a wrapper around calcGef
myFoo <- function(ab, datos)
{
    cat(ab, '\n')
    p <- calcGef(lat = latitud,
                modeRad = 'bdI', dataRad = datos, 
                corr = 'BRL',
                horizBright = FALSE,
                alfa = ab[1], beta = ab[2])
    ## Return the aggregated value of calcGef
    sum(as.data.frameD(p, complete = TRUE)$Gefd)/1000
}

## Apply the function to each combination of alpha/beta
M$Gef <- apply(M, 1, myFoo, datos = NRELMeteo)

levelplot(Gef/max(Gef) ~ alfa * beta, data = M,
          contour = TRUE,
          par.settings = custom.theme.2)

x

On the other hand, you should consider this simpler version of the funcion to read your data:

library("solaR")

readSIAR <- function(x)
{
    datos <- read.csv2(x, fileEncoding = 'UTF-16LE')
    datos <- datos[, c("Fecha", "HoraMin",
                           "Radiación..W.m2.", "Temp.Media..ºC.")]
    tt <- with(datos,
               as.POSIXct(paste0(
                   Fecha, ' ',
                   paste(HoraMin%/%100, HoraMin%%100, '00', sep = ':')),
                   format = '%d/%m/%Y %H:%M:%S'))
    datos$time <- local2Solar(tt, longitud)
    datos[, c("Fecha", "HoraMin")] <- NULL
    names(datos) <- c("G0", "Ta", "time")
    datos
}

pica02 <- readSIAR('/tmp/datos/V17_Picassent_01_01_2002_31_12_2002.csv')

xyplot(G0 ~ as.numeric(format(time, "%H%M"))/100 - 12, data = pica02, grid = TRUE)

meteo02 = dfI2Meteo(file = pica02,
                    lat = latitud,
                    time.col = 'time',
                    source = 'Picassent')
xyplot(meteo02, grid = TRUE)
oscarperpinan commented 6 years ago

I close this issue because there has been no answer after several months.