Recurso eólico. Cubos espacio-temporales

Librerias y Datos

library(raster) # rasterbrick (datos raster)
library(ncdf4) # datos netCDF
library(rasterVis) # visualización multitemporal
library(sf) # datos vectoriales
library(latticeExtra) # visualización multitemporal

# Visualización de flechas
library(rWind)
library(fields)
library(shape)

# Rosa de los vientos
library(tibble) # tipo de dato tibble
library(clifro)

# Histogramas mensuales
library(reshape)
library(ggplot2)

# Visualización de series temporales
library(highcharter)

# Agregación temporal
library(xts)

# Distribución de Weibull
library(fitdistrplus)

Área de Estudio

Subiremos dos conjuntos de datos de polígonos. Uno en formato RData con la forma de Ecuador, y otro en formato shp con la forma de la cuenca del río Paute.

Ecuador <- load("shp/Ecuador_shape_4326.RData")
Ecuador <- ecuador
crs(Ecuador) <- crs("+proj=longlat +datum=WGS84 +no_defs")

Paute <- st_read("shp/RioPaute.shp", crs = st_crs("+proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs")) # SRS 32717 # libreria sf

## Reading layer `RioPaute' from data source 
##   `C:\Users\HP_USER\Desktop\Desktop\CURSO UDA (EOLICO)\11_Cubos espacio-temporales para el recurso eólico, gestión espacial y temporal\cubos_espacio-temporales\cubos_espacio-temporales\shp\RioPaute.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 681758 ymin: 9637966 xmax: 803433.7 ymax: 9745363
## Projected CRS: +proj=utm +zone=17 +south +datum=WGS84 +units=m +no_defs

Paute <- st_transform(Paute, crs = st_crs("+proj=longlat +datum=WGS84 +no_defs")) # transformar el sistema de referencia espacial
Paute <- as(Paute, "Spatial") # libreria sp

ERA5 Land.c componentes u y v

Resolución espacial de 0.1 grados aproximadamente, equivalente a unos 9 km
Componentes u y v a una altura de 10 m
Datos disponibles a nivel mensual (también disponibles a nivel horario y promedio mensual por hora)
ERA5 a una resolución de aproximadamente 0.25 grados, equivalente a unos 28 km, disponible en 37 niveles de presión.

V Meridional (N a S)

Componente hacia el Norte del viento a 10 metros. Es la velocidad horizontal del aire moviéndose hacia el norte, a una altura de diez metros sobre la superficie de la Tierra, en metros por segundo. Se debe tener cuidado al comparar esta variable con observaciones, ya que las observaciones del viento varían en pequeñas escalas espaciales y temporales y se ven afectadas por el terreno local, la vegetación y los edificios que solo se representan en promedio en el Sistema Integrado de Pronóstico de ECMWF. Esta variable se puede combinar con la componente U del viento a 10 metros para obtener la velocidad y dirección del viento horizontal a 10 metros.

U Zonal (W a E, O a E)

Componente hacia el Este del viento a 10 metros. Es la velocidad horizontal del aire moviéndose hacia el este, a una altura de diez metros sobre la superficie de la Tierra, en metros por segundo. Se debe tener cuidado al comparar esta variable con observaciones, ya que las observaciones del viento varían en pequeñas escalas espaciales y temporales y se ven afectadas por el terreno local, la vegetación y los edificios que solo se representan en promedio en el Sistema Integrado de Pronóstico de ECMWF. Esta variable se puede combinar con la componente V del viento a 10 metros para obtener la velocidad y dirección del viento horizontal a 10 metros.

Haremos uso de dos tipos de datos ERA5. Primero, un archivo .tif con las componentes u y v para enero de 2005. Esto se utilizará para trabajar con un mes específico. En segundo lugar, un archivo .nc se utilizará para trabajar con una serie temporal de imágenes. Estos datos constituyen un cubo de datos espacio-temporales. Exploraremos primero la dimensión espacial trabajando con un mes y una serie temporal de meses, y más tarde la dimensión temporal trabajando con un sitio.

Spatio-temporal datacube

1. Trabajando con un mes (1 mapa)

Leer los datos

v<-raster("data_demo/v_wind_ERA5LAND201501.tif")
u<-raster("data_demo/u_wind_ERA5LAND201501.tif")

u

## class      : RasterLayer 
## dimensions : 101, 71, 7171  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : -82.05, -74.95, -7.05, 3.05  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : u_wind_ERA5LAND201501.tif 
## names      : u_component_of_wind_10m

## class      : RasterLayer 
## dimensions : 101, 71, 7171  (nrow, ncol, ncell)
## resolution : 0.1, 0.1  (x, y)
## extent     : -82.05, -74.95, -7.05, 3.05  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : v_wind_ERA5LAND201501.tif 
## names      : v_component_of_wind_10m

Recortar y enmascarar las imágenes con la forma de Ecuador. Visualizar las componentes u y v del viento.

Puedes probar diferentes paletas de colores: pasmaTheme, YlOrRdTheme, BuRdTheme, RdBuTheme, GrTheme, and BTCTheme.

v <-mask(crop(v, Ecuador), Ecuador)
u <-mask(crop(u, Ecuador), Ecuador)

levelplot(u, par.settings = infernoTheme, margin= FALSE, main="U component at 10m") + latticeExtra::layer(sp.polygons(Ecuador, col='black'))

levelplot(v, par.settings = infernoTheme, margin= FALSE, main="V component at 10m") + latticeExtra::layer(sp.polygons(Ecuador, col='black'))

Velocidad y dirección del viento

Las componentes U y V son coordenadas cartesianas, por lo tanto, necesitamos combinarlas en coordenadas polares, es decir, longitud del vector (velocidad) y dirección (desde el norte 0°).

https://confluence.ecmwf.int/pages/viewpage.action?pageId=133262398

Para esto, utilizaremos las siguientes fórmulas:

      wind_speed=sqrt(u**2+v**2)
      
      wind_direction=180+180/pi atan2(v,u)

speed <- sqrt(v**2+u**2)
direction<- 180+(180/pi) * atan2(v,u)

levelplot(speed, par.settings = YlOrRdTheme, margin= FALSE, main="Speed at 10m m/s") + latticeExtra::layer(sp.polygons(Ecuador, col='black'))

levelplot(direction, par.settings = BuRdTheme, margin= FALSE, main="Direction 0-360  at 10m") + latticeExtra::layer(sp.polygons(Ecuador, col='black'))

También exploremos una forma más tradicional de visualizar datos de viento con flechas. Para esto, superpondremos las flechas de direction sobre la rasterización de speed. Para hacerlo, direction debe ser un data.frame (es decir, una tabla).

dir<- as.data.frame(direction, xy=TRUE)
names(dir)<- c("x", "y", "dir")

image.plot(speed, col=bpy.colors(1000, alpha=0.8),zlim=c(0,7),xlab="longitude", ylab="latitude")
Arrowhead(dir$x, dir$y, angle=arrowDir(dir), arr.length=0.1, arr.type="curved")

dir_u<- as.data.frame(u, xy=TRUE)
names(dir_u)<- c("x", "y", "u")

dir_v<- as.data.frame(v, xy=TRUE)
names(dir_v)<- c("x", "y", "v")
      
speed_df <- as.data.frame(speed, xy=TRUE)
names(speed_df)<- c("x", "y", "speed")

dir$u <- dir_u$u
dir$v <- dir_v$v


ggplot() +
  geom_raster(data = speed_df, aes(x = x, y = y, fill = speed)) +
  geom_segment(data = dir, aes(x = x, y = y, xend = x + u*0.05, yend = y + v*0.05, color="white"),
               arrow = arrow(type = "closed", length = unit(1.5, "mm")),
               linewidth = 0.1, lineend = "round", linejoin = "round") +
  
  labs(title = "Visualización tradicional de datos de viento con flechas en Ecuador") +
  theme_minimal()

## Warning: Removed 1837 rows containing missing values (`geom_segment()`).

No estoy satisfecha con este mapa. Las flechas son grandes y no permiten ver claramente el mapa de velocidad. Intentemos hacer este mismo mapa para un área más pequeña, como la cuenca de Paute. Para esto, nuestras capas raster se recortarán con la forma de Paute.

speed_paute <-mask(crop(speed, Paute), Paute)
dir_paute <-mask(crop(direction, Paute), Paute)

dir.df.paute<- as.data.frame(dir_paute, xy=TRUE)
names(dir.df.paute)<- c("x", "y", "dir")

image.plot(speed_paute, col=bpy.colors(1000, alpha=0.8),zlim=c(0,2),xlab="longitude", ylab="latitude")
Arrowhead(dir.df.paute$x, dir.df.paute$y, angle=
arrowDir(dir.df.paute), arr.length=0.2, arr.type="curved")

speed_paute <-mask(crop(speed, Paute), Paute)
dir_u_paute <-mask(crop(u, Paute), Paute)
dir_v_paute <-mask(crop(v, Paute), Paute)

dir_u.df.paute<- as.data.frame(dir_u_paute, xy=TRUE)
names(dir_u.df.paute)<- c("x", "y", "u")

dir_v.df.paute<- as.data.frame(dir_v_paute, xy=TRUE)
names(dir_v.df.paute)<- c("x", "y", "v")

dir_u.df.paute$u <- dir_u.df.paute$u
dir_u.df.paute$v <- dir_v.df.paute$v

speed_paute.df<- as.data.frame(speed_paute, xy=TRUE)
names(speed_paute.df)<- c("x", "y", "speed")

ggplot() +
  geom_raster(data = speed_paute.df, aes(x = x, y = y, fill = speed)) +
  geom_segment(data = dir_u.df.paute, aes(x = x, y = y, xend = x + u*0.03, yend = y + v*0.03),
               arrow = arrow(type = "closed", length = unit(1, "mm")),
               linewidth = 0.1) +
  
  labs(title = "Visualización tradicional de datos de viento con flechas en Ecuador") +
  theme_minimal()

## Warning: Removed 50 rows containing missing values (`geom_segment()`).

Densidad de energía eólica

Medida de la energía eólica disponible en un sitio específico (W/m2)

Adaptado de <https://conwx.com/why-invest-in-high-quality-power-production-forecasts/>
Explicación de la Formula https://www.youtube.com/watch?v=gHm7nPWYveo

Densidad de energía eólica se calcula a partir de la velocidad del viento y densidad del aire. Ya tenemos la velocidad del viento. En cuanto a la densidad del aire, existen varios enfoques para calcularla. Aquí, con el fin de simplificar el guion, utilizaremos un valor constante. Algunas alternativas son:

Para aproximar la densidad del aire ρ, se podría utilizar la fórmula de altitud barométrica para temperatura constante. ρ=ρ0exp(−z/Hρ) donde ρ0=1.225kg/m3 es la densidad atmosférica estándar al nivel del mar y temperatura estándar, z es la altitud y Hρ=8.55km es la altitud de escala para la densidad. Tutorial 5 - Desarrollar una aplicación basada en el ejemplo de un indicador climático — Documentación de la Caja de herramientas del Almacén de Datos Climáticos 1.1.1”

Otro enfoque. Recomendado si tienes los datos.
Atlas Eólico Global. Capa de densidad del aire a diferentes altitudes

Wind Power wp = windspeed ** 3 * air_density * 0.5Units W/m**2

air_density = 1.225  # NOTA importante. Esto es simplemente para facilitar el análisis en este paso. 1.225 kg/m³ es la densidad atmosférica estándar al nivel del mar y a una temperatura estándar. En Ecuador, la densidad del aire cambia con la altitud, por lo que debería utilizarse un mapa de densidad del aire.

wp = speed ** 3 * air_density * 0.5

levelplot(wp, par.settings = YlOrRdTheme, margin= FALSE, main="Wind Power at 10m W/m**2") + latticeExtra::layer(sp.polygons(Ecuador, col='black'))

Rosa de los Vientos

Primero se crea un data.frame (tabla) llamado df para almacenar los valores de velocidad y dirección.

df<-NULL
df$s <- values(speed)
df$d <-  values(direction)


windrose(speed = df$s,
                 direction = df$d,
                 n_directions = 12,
                 speed_cuts = seq(0,4,1),
                 ggtheme='minimal',
                 col_pal = 'YlGnBu')

2.Series temporales de imágenes

Leer datos

Lee los datos en formato nc. Las imágenes están apiladas (brick) en una única variable. Se grafican las primeras 4 imágenes del componente v.

wind_v<-brick("data_demo/ERA_viento_u_2015_2019.nc")
wind_u<-brick("data_demo/ERA_viento_v_2015_2019.nc")

wind_v

## class      : RasterBrick 
## dimensions : 40, 28, 1120, 60  (nrow, ncol, ncell, nlayers)
## resolution : 0.25, 0.25  (x, y)
## extent     : -82, -75, -7, 3  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : ERA_viento_u_2015_2019.nc 
## names      :        X1,        X2,        X3,        X4,        X5,        X6,        X7,        X8,        X9,       X10,       X11,       X12,       X13,       X14,       X15, ... 
## min values : -2.166206, -1.798296, -2.613769, -3.230329, -3.508509, -3.534022, -3.161539, -3.269866, -2.992448, -2.507188, -2.538931, -2.845659, -2.119399, -3.636128, -2.674843, ... 
## max values :  4.080397,  3.026091,  2.083539,  2.411764,  2.930612,  3.899724,  4.376393,  5.102653,  4.972281,  4.894050,  4.815600,  4.574602,  3.557852,  2.461947,  1.925570, ... 
## z (unknown): 1, 60 (min, max)
## varname    : variable

levelplot(wind_v[[1:4]], main="V component at 10m")+ latticeExtra::layer(sp.polygons(Ecuador, col='black'))

Mejoremos algunas cosas. Primero, renombraremos cada imagen según la fecha. Segundo, enmascararemos la imagen con la forma de Ecuador.

date<- seq.Date(as.Date("2015-01-01"), as.Date("2019-12-01"), "months")
names(wind_v)<- date
names(wind_u)<- date

levelplot(wind_v[[1:4]],main="V component at 10m")+ latticeExtra::layer(sp.polygons(Ecuador, col='black'))

levelplot(wind_v[[c(1,13,25,37,49)]],main="V component at 10m")+ latticeExtra::layer(sp.polygons(Ecuador, col='black'))

wind_v <-crop(wind_v, Ecuador)
wind_u <-crop(wind_u, Ecuador)

levelplot(wind_v[[1:4]], main="Wind-v")+ latticeExtra::layer(sp.polygons(Ecuador, col='black'))

levelplot(wind_u[[1:4]], main="Wind-u")+ latticeExtra::layer(sp.polygons(Ecuador, col='black'))

Velocidad y dirección del viento

La misma fórmula utilizada para una sola imagen se puede aplicar al bloque raster

speed <- sqrt(wind_v**2+wind_u**2)
direction<- 180+(180/pi) * atan2(wind_v,wind_u)

levelplot(speed[[1:4]], main="Wind speed at 10m m/s", par.settings = YlOrRdTheme, margin= FALSE) + latticeExtra::layer(sp.polygons(Ecuador, col='black'))

levelplot(direction[[1:4]], par.settings = BuRdTheme, margin= FALSE, main="Direction 0-360°") + latticeExtra::layer(sp.polygons(Ecuador, col='black'))

Guardar archivo en el disco

# writeRaster(speed[[1:12]], filename = "wind_2015.nc", format = "CDF", overwrite = TRUE)

# brick("wind_2015.nc")

# writeRaster(speed[[1]], filename = "ERA5wind_2015_01.nc", format = "CDF", overwrite = TRUE)
# writeRaster(speed[[2]], filename = "ERA5wind_2015_02.nc", format = "CDF", overwrite = TRUE)
# writeRaster(speed[[3]], filename = "ERA5wind_2015_03.nc", format = "CDF", overwrite = TRUE)
# writeRaster(speed[[4]], filename = "ERA5wind_2015_04.nc", format = "CDF", overwrite = TRUE)
# writeRaster(speed[[5]], filename = "ERA5wind_2015_05.nc", format = "CDF", overwrite = TRUE)
# writeRaster(speed[[6]], filename = "ERA5wind_2015_06.nc", format = "CDF", overwrite = TRUE)
# 
# 
# writeRaster(speed[[7]], filename = "ERA5wind_2015_07.nc", format = "CDF", overwrite = TRUE)
# writeRaster(speed[[8]], filename = "ERA5wind_2015_08.nc", format = "CDF", overwrite = TRUE)
# writeRaster(speed[[9]], filename = "ERA5wind_2015_09.nc", format = "CDF", overwrite = TRUE)
# writeRaster(speed[[10]], filename = "ERA5wind_2015_10.nc", format = "CDF", overwrite = TRUE)
# writeRaster(speed[[11]], filename = "ERA5wind_2015_11.nc", format = "CDF", overwrite = TRUE)
# writeRaster(speed[[12]], filename = "ERA5wind_2015_12.nc", format = "CDF", overwrite = TRUE)

Explora el histograma de la velocidad del viento de cada mes. Primero se realizará esto con la función hist de la biblioteca raster (¡No homogeneizó el eje y!!), y luego con la biblioteca ggplot. Se utilizarán los primeros 12 meses.

#Option 1
hist(speed[[1:12]])

#Option 2
df <- as.data.frame(speed[[1:12]], xy = TRUE) %>%
    melt(id.vars = c('x','y'))

df

ggplot(df) +
  geom_histogram(aes(value)) +
    facet_wrap(~variable)

Climatologia

months <-  as.integer(format(date, "%m"))
speed_climatoloty <- stackApply(speed, months, fun=mean)
direction_climatoloty <- stackApply(direction, months, fun=mean)

names(speed_climatoloty)<- month.abb
names(direction_climatoloty)<- month.abb

levelplot(speed_climatoloty, par.settings = YlOrRdTheme, margin= FALSE, main="Climatology speed m/s") + latticeExtra::layer(sp.polygons(Ecuador, col='black'))

levelplot(direction_climatoloty, par.settings = BuRdTheme, margin= FALSE, main="Climatology direction 0-360°") + latticeExtra::layer(sp.polygons(Ecuador, col='black'))

Climatologia de rosa de los vientos

df <- as.data.frame(speed_climatoloty, xy = TRUE) %>%
    melt(id.vars = c('x','y')) 

df.d <- as.data.frame(direction_climatoloty, xy = TRUE) %>%
    melt(id.vars = c('x','y'))

df$direction<-df.d$value
df

windrose(speed = df$value,
                 direction = df$direction,
                 n_directions = 12,
                 speed_cuts = seq(0,4,1),
                  facet= df$variable,
                 ggtheme='minimal',
                 col_pal = 'YlGnBu',
                 n_col=4)

3.Una ubicación específica

Extraer y visualizar series temporales.

Seleccionaremos las coordenadas del píxel número 144 y veremos su ubicación

xy<- xyFromCell(speed, 423)
xy <- SpatialPoints(xy)

#extraer por coordenada
#xy2 <- cbind(x = tu_valor_x, y = tu_valor_y)
#xy2 <- SpatialPoints(xy2)
ts.speed <- c(extract(speed, xy))

levelplot(speed[[1]], par.settings = YlOrRdTheme, margin= FALSE, main="Speed m/s") + latticeExtra::layer(sp.polygons(Ecuador, col='black'))+ latticeExtra::layer(sp.points(xy, col='red',pch = 19))

Ahora la serie temporal se extrae y se define como un tipo de dato temporal.

ts.direction <- c(extract(direction, xy))

plot(ts.speed, type="l", main="Speed")

start.date <- 2015 # Fecha de descarga final
finish.date <- 2019 # Fecha de descarga final

wind_s.ts <- ts(ts.speed,start=c(start.date, 1), end=c(finish.date, 12), frequency=12)

wind_d.ts <- ts(ts.direction,start=c(start.date, 1), end=c(finish.date, 12), frequency=12)

Usando la biblioteca highcharter, podemos obtener una visualización muy agradable de la serie temporal

highchart(type = "stock")%>%hc_add_series(wind_s.ts, name="Speed")

highchart(type = "stock")%>%hc_add_series(wind_d.ts, name="Direction")

Agregación Temporal

#Option 1 ts library. Seasonality Q1, Q2, Q3, Q4
wind_s.agg<- aggregate(wind_s.ts, nfrequency= 4, mean) # estacionalidad

wind_s.agg

##           Qtr1      Qtr2      Qtr3      Qtr4
## 2015 1.6425545 2.7999760 3.1471632 2.2009832
## 2016 1.7087085 2.4343747 2.9945822 0.9347907
## 2017 0.6539155 1.9412067 2.9283664 1.2596254
## 2018 1.3845367 2.0477719 3.0175552 1.1715427
## 2019 1.8665203 2.2756248 3.1290299 1.4944323

#Option 2 xts library
apply.quarterly(as.xts(wind_s.ts), mean)

##               [,1]
## Mar 2015 1.6425545
## Jun 2015 2.7999760
## Sep 2015 3.1471632
## Dec 2015 2.2009832
## Mar 2016 1.7087085
## Jun 2016 2.4343747
## Sep 2016 2.9945822
## Dec 2016 0.9347907
## Mar 2017 0.6539155
## Jun 2017 1.9412067
## Sep 2017 2.9283664
## Dec 2017 1.2596254
## Mar 2018 1.3845367
## Jun 2018 2.0477719
## Sep 2018 3.0175552
## Dec 2018 1.1715427
## Mar 2019 1.8665203
## Jun 2019 2.2756248
## Sep 2019 3.1290299
## Dec 2019 1.4944323

apply.yearly(as.xts(wind_s.ts), mean)

##              [,1]
## Dec 2015 2.447669
## Dec 2016 2.018114
## Dec 2017 1.695778
## Dec 2018 1.905352
## Dec 2019 2.191402

Inter-annual variability

barplot(apply.yearly(as.xts(wind_s.ts), mean))

Intra-annual variability

boxplot(wind_s.ts ~ cycle(wind_s.ts))

boxplot(wind_d.ts ~ cycle(wind_d.ts))

library(clifro) # windrose plot

df <- as.data.frame(speed_climatoloty, xy = TRUE) %>%
    melt(id.vars = c('x','y')) # avoid NULL values by cropping the images and not masking

df.d <- as.data.frame(direction_climatoloty, xy = TRUE) %>%
    melt(id.vars = c('x','y'))

df$direction<-df.d$value
df

windrose(speed = df$value,
                 direction = df$direction,
                 n_directions = 12,
                 speed_cuts = seq(0,4,1),
                  facet= df$variable,
                 ggtheme='minimal',
                 col_pal = 'YlGnBu',
                 n_col=4)

Weibull distribution

#http://www.di.fc.ul.pt/~jpn/r/distributions/fitting.html

plotdist(df$value, histo=TRUE, demp=TRUE)

# descdist(df$value)

weibull_fit<-  fitdist(df$value, "weibull")
summary(weibull_fit)#av_ws is the wind speed data

## Fitting of the distribution ' weibull ' by maximum likelihood 
## Parameters : 
##        estimate  Std. Error
## shape 0.8997051 0.007717847
## scale 1.0144725 0.014146746
## Loglikelihood:  -7634.787   AIC:  15273.57   BIC:  15287.33 
## Correlation matrix:
##           shape     scale
## shape 1.0000000 0.3389739
## scale 0.3389739 1.0000000

# plot(weibull_fit)
denscomp(weibull_fit)

Weibull bondad de ajuste Goodness-of-fit

fg<-fitdist(df$value, "gamma")
fln<-fitdist(df$value, "lnorm")

par(mfrow=c(2,2))
plot.legend<-c("Weibull","lognormal","gamma")
denscomp(list(weibull_fit,fln,fg), legendtext=plot.legend)

gofstat(list(weibull_fit, fg, fln), fitnames = c("weibull", "gamma", "lnorm"))

## Goodness-of-fit statistics
##                                   weibull       gamma       lnorm
## Kolmogorov-Smirnov statistic   0.09827509   0.1172141  0.06725339
## Cramer-von Mises statistic    29.35332850  41.1417252  9.38214652
## Anderson-Darling statistic   190.09697797 237.1824699 62.25942820
## 
## Goodness-of-fit criteria
##                                 weibull    gamma    lnorm
## Akaike's Information Criterion 15273.57 15414.12 13655.31
## Bayesian Information Criterion 15287.33 15427.87 13669.07

4.Comparación de sitios

xy2<- xyFromCell(speed, 144)
xy2 <- SpatialPoints(xy2)

levelplot(speed[[1]], par.settings = YlOrRdTheme, margin= FALSE, main="Speed m/s") + latticeExtra::layer(sp.polygons(Ecuador, col='black'))+ latticeExtra::layer(sp.points(xy, col='red',pch = 19))+ latticeExtra::layer(sp.points(xy2, col='red',pch = 19))

ts.speed2 <- c(extract(speed, xy2))
ts.direction2 <- c(extract(direction, xy2))

wind_s.ts2 <- ts(ts.speed2,start=c(start.date, 1), end=c(finish.date, 12), frequency=12)
wind_d.ts2 <- ts(ts.direction2,start=c(start.date, 1), end=c(finish.date, 12), frequency=12)

highchart(type = "chart")%>% hc_add_series(tapply(wind_s.ts2,cycle(wind_s.ts2), mean), name="Point2")%>% 
hc_add_series(tapply(wind_s.ts,cycle(wind_s.ts), mean), name="Point1")

Recurso eólico. Cubos espacio-temporales

Daniela Ballari

Librerias y Datos

Área de Estudio

ERA5 Land.c componentes u y v

V Meridional (N a S)

U Zonal (W a E, O a E)

Spatio-temporal datacube

1. Trabajando con un mes (1 mapa)

Leer los datos

Velocidad y dirección del viento

Densidad de energía eólica

Rosa de los Vientos

2.Series temporales de imágenes

Leer datos

Velocidad y dirección del viento

Climatologia

3.Una ubicación específica

Extraer y visualizar series temporales.

Agregación Temporal

Inter-annual variability

Intra-annual variability

Weibull distribution

Weibull bondad de ajuste Goodness-of-fit

4.Comparación de sitios