Mostrando entradas con la etiqueta seq. Mostrar todas las entradas
Mostrando entradas con la etiqueta seq. Mostrar todas las entradas

2021-01-01

Cómo crear una serie temporal con intervalos de 30 minutos

Title

Problema

Queremos crear una serie temporal de intervalos de 30 minutos.

Ejemplo

[1] "2017-01-01 00:00:00 UTC"
[2] "2017-01-01 00:30:00 UTC"
[3] "2017-01-01 01:00:00 UTC"
[4] "2017-01-01 01:30:00 UTC"
[5] "2017-01-01 02:00:00 UTC"
[6] "2017-01-01 02:30:00 UTC"

Solución

Usamos la función seq y especificando los minutos en el argumento by, así como el uso horario "UTC". Tecleamos ?seq.POSIXt para obtener más detalles sobre el argumento que podemos especificar como cadena sde texto:

A character string, containing one of "sec", "min", "hour", "day", "DSTday", "week", "month", "quarter" or "year". This can optionally be preceded by a (positive or negative) integer and a space, or followed by "s".

seq(as.POSIXct("2017-01-01", tz = "UTC"),
    as.POSIXct("2017-01-02", tz = "UTC"),
    by = "30 min")

Referencias

2020-12-23

How to create a time series by 30 minute intervals

Title

Problem

We want to create a time series by 30 minute intervals.

Example

[1] "2017-01-01 00:00:00 UTC"
[2] "2017-01-01 00:30:00 UTC"
[3] "2017-01-01 01:00:00 UTC"
[4] "2017-01-01 01:30:00 UTC"
[5] "2017-01-01 02:00:00 UTC"
[6] "2017-01-01 02:30:00 UTC"

Solution

We use the function seq and specify minutes in the by argument, and pass the time zone "UTC". Type ?seq.POSIXt for more details about the by argument specified as a character string:

A character string, containing one of "sec", "min", "hour", "day", "DSTday", "week", "month", "quarter" or "year". This can optionally be preceded by a (positive or negative) integer and a space, or followed by "s".

seq(as.POSIXct("2017-01-01", tz = "UTC"),
    as.POSIXct("2017-01-02", tz = "UTC"),
    by = "30 min")

References

2019-02-08

Calcular y representar la duración del día con R

Problema

Queremos calcular y representar la duración del día con R en función de unas coordinadas geográficas.

Solución

  • 1. Calculamos la salida y puesta de sol
  • Primero calculamos la salida y la puesta de sol con la función getSunlightTimes del paquete suncalc. Indicamos el intervalo deseado, las coordinadas (latitud y longitud), y el huso horario (tz, time zone) correspondiente.

    library(suncalc) 
    library(tidyverse)
    library(scales)
    df <-
      getSunlightTimes(
        date = seq.Date(as.Date("2017-12-01"), as.Date("2018-12-31"), by = 1),
        keep = c("sunrise", "sunriseEnd", "sunset", "sunsetStart"),
        lat = 39.8628,
        lon = 4.0273,
        tz = "CET"
      )
    
  • 2. Gráfico de la salida y puesta de sol
  • Necesitamos manipular los datos originales para calcular la diferencia entre el inicio del día, y la salida y puesta de sol. Después representamos las dos nuevas variables usando geom_ribbon de ggplot2. Luego personalizamos los ejes y el título.

    # Amanecer/ocaso
    df %>%
      mutate(
        date = as.POSIXct(date) - 12 * 60 * 60 ,
        sunrise = sunrise - date,
        sunset =  sunset - date,
      ) %>%
      ggplot() +
      geom_ribbon(aes(x = date, ymin = sunrise, ymax = sunset),
                  fill = "#FDE725FF",
                  alpha = .8) + # "#ffeda0"
      scale_x_datetime(
        breaks = seq(as.POSIXct(min(df$date)), as.POSIXct(max(df$date)), "month"),
        expand = c(0, 0),
        labels = date_format("%b %y"),
        minor_breaks = NULL
      ) +
      scale_y_continuous(
        limits = c(0, 24),
        breaks = seq(0, 24, 2),
        expand = c(0, 0),
        minor_breaks = NULL
      ) +
      labs(
        x = "Date",
        y = "Hours",
        title = sprintf(
          "Sunrise and Sunset for %s\n%s ",
          "Toledo (Spain)",
          paste0(as.Date(range(df$date)), sep = " ", collapse = "to ")
        )
      ) +
      theme(
        panel.background = element_rect(fill = "#180F3EFF"),
        panel.grid = element_line(colour = "grey", linetype = "dashed")
      )
    
  • 3. Duración del día
  • Muy similar al gráfico anterior. Ahora solamente necesitamos calcular la duración del día day_length y representar los resultados con geom_area and geom_line.

    df %>%
      mutate(
        date = as.POSIXct(date),
        day_length = as.numeric(sunset - sunrise)
      ) %>%
      ggplot(aes(x = date, y = day_length)) +
      geom_area(fill = "#FDE725FF", alpha = .4) +
      geom_line(color = "#525252") +
      scale_x_datetime(
        expand = c(0, 0),
        labels = date_format("%b '%y"),
        breaks =  seq(as.POSIXct(min(df$date)), as.POSIXct(max(df$date)), "month"),
        minor_breaks = NULL
      ) +
      scale_y_continuous(
        limits = c(0, 24),
        breaks = seq(0, 24, 2),
        expand = c(0, 0),
        minor_breaks = NULL
      ) +
      labs(x = "Date", y = "Hours", title = "Toledo (Spain) - Daytime duration") +
      theme_bw()
    

    Entradas relacionadas

    Referencias

    2019-02-04

    Calculate and plot sunrise and sunset times with R

    Problem

    We would like to calculate and plot the sunrise and sunset times based on any location's latitude and longitude coordinates with R.

    Solution

  • 1. Compute sunrise and sunset times
  • First we calculate the sunrise and sunset times using the function getSunlightTimes from the package suncalc. We pass the desired date interval, the appropiate latitude and longitude coordinates, and time zone (tz).

    library(suncalc) 
    library(tidyverse)
    library(scales)
    df <-
      getSunlightTimes(
        date = seq.Date(as.Date("2017-12-01"), as.Date("2018-12-31"), by = 1),
        keep = c("sunrise", "sunriseEnd", "sunset", "sunsetStart"),
        lat = 39.8628,
        lon = 4.0273,
        tz = "CET"
      )
    
  • 2. Sunrise and sunset times plot
  • We need to manipulate the original data frame to calculate the difference between midnight start of day, and the sunrise and sunset times. Then we plot those two new variables using geom_ribbon from ggplot2. We further customize the axes, and title.

    # Sunrise/set
    df %>%
      mutate(
        date = as.POSIXct(date) - 12 * 60 * 60 ,
        sunrise = sunrise - date,
        sunset =  sunset - date,
      ) %>%
      ggplot() +
      geom_ribbon(aes(x = date, ymin = sunrise, ymax = sunset),
                  fill = "#FDE725FF",
                  alpha = .8) + # "#ffeda0"
      scale_x_datetime(
        breaks = seq(as.POSIXct(min(df$date)), as.POSIXct(max(df$date)), "month"),
        expand = c(0, 0),
        labels = date_format("%b %y"),
        minor_breaks = NULL
      ) +
      scale_y_continuous(
        limits = c(0, 24),
        breaks = seq(0, 24, 2),
        expand = c(0, 0),
        minor_breaks = NULL
      ) +
      labs(
        x = "Date",
        y = "Hours",
        title = sprintf(
          "Sunrise and Sunset for %s\n%s ",
          "Toledo (Spain)",
          paste0(as.Date(range(df$date)), sep = " ", collapse = "to ")
        )
      ) +
      theme(
        panel.background = element_rect(fill = "#180F3EFF"),
        panel.grid = element_line(colour = "grey", linetype = "dashed")
      )
    
  • 3. Daytime duration
  • Very similar to the preceding plot. This time we only need to calculate the day_length and plot the results using geom_area and geom_line.

    df %>%
      mutate(
        date = as.POSIXct(date),
        day_length = as.numeric(sunset - sunrise)
      ) %>%
      ggplot(aes(x = date, y = day_length)) +
      geom_area(fill = "#FDE725FF", alpha = .4) +
      geom_line(color = "#525252") +
      scale_x_datetime(
        expand = c(0, 0),
        labels = date_format("%b '%y"),
        breaks =  seq(as.POSIXct(min(df$date)), as.POSIXct(max(df$date)), "month"),
        minor_breaks = NULL
      ) +
      scale_y_continuous(
        limits = c(0, 24),
        breaks = seq(0, 24, 2),
        expand = c(0, 0),
        minor_breaks = NULL
      ) +
      labs(x = "Date", y = "Hours", title = "Toledo (Spain) - Daytime duration") +
      theme_bw()
    

    Related posts

    References

    2018-09-23

    Gráficos de An Introduction to Statistical Learning - Figura 2.6.

    Gráfico a replicar

    Continuamos con la serie iniciada sobre la creación de gráficos del libro An Introduction to Statistical Learning. En esta ocasión replicaremos los gráficos de la figura 2.6. Utiliza el conjunto de datos Income2 (renta). El gráfico representa los valores observados de renta (en miles de dólares) en función de los años de educación y la antigüedad de 30 individuos. Los puntos rojos representan los valores observados para esas tres variables. La superficie amarilla representa un modelo de suavizado con thin-plate spline. En este caso no hay líneas verticales negras que representaban el error asociado con cada observación, pues la superficie thin-plate spline se adapta perfectamente a los datos.

    Solución

    Como utilizamos el paquete mgcv que se utiliza para regresión no para interpolación por lo que tenemos que ajustarlo para que la superficie thin plate spline pase por todos los puntos.

    1. Desactivamos la aproximación de bajo rango que utiliza bs = 'tp' estableciendo el parámetro k con el número exacto de datos de la muestra.
    2. xt <- unique(income_2[c("Education", "Seniority")]) 
      nrow(xt)
      
      [1] 30
      
    3. Empleamos sp = 0 para desactivar la penalización de la spline.

  • Primera aproximación
  • model <- gam(Income ~ s(Education, Seniority, k = 30, sp = 0),
                               data = income_2)
    x <- range(income_2$Education)
    x <- seq(x[1], x[2], length.out=30)
    y <- range(income_2$Seniority)
    y <- seq(y[1], y[2], length.out=30)
    z <- outer(x, y, 
               function(Education,Seniority)
                         predict(model, data.frame(Education,Seniority)))
    p <- persp(x, y, z, theta = 30, phi = 30, 
               col = "yellow", expand = 0.5, shade = 0.2, 
               xlab = "Education", ylab = "Seniority", zlab = "Income")
    obs <-  trans3d(income_2$Education, income_2$Seniority, income_2$Income, p)
    pred <-  trans3d(income_2$Education, income_2$Seniority, fitted(model), p)
    points(obs, col = "red", pch = 16)
    segments(obs$x, obs$y, pred$x, pred$y)
    
    Como se puede apreciar en el gráfico, hemos eliminado los errores y se adapta perfectamente a los puntos. Sin embargo, se puede apreciar que la superficie es muy ondulada pues spline es isotrópica o radial. Como dos variables difieren considerablemente en escala, necesitaremos estandarizarlas.

    with(income_2, plot(Education, Seniority, asp = 1))
    
    Las estandarizamos y representamos nuevamente:

    xt_scaled <- scale(xt)
    dat <- data.frame(xt_scaled, Income = income_2$Income)
    with(dat, plot(Education, Seniority, asp = 1))
    
    # Ajustamos el modelo a los datos escalados
    interpolation_model <- gam(Income ~ s(Education, Seniority, k = 30, sp = 0),
                               data = dat)
    # Creamos las coordenadas para el gráfico
    x <- range(dat$Education)
    x <- seq(x[1], x[2], length.out=30)
    y <- range(dat$Seniority)
    y <- seq(y[1], y[2], length.out=30)
    z <- outer(x, y, 
               function(Education,Seniority)
                         predict(interpolation_model, data.frame(Education,Seniority)))
    
    # Volvemos a transformar las coordenadas x e y a su escala original.
    # No es necesario transformar los valores esperados (predicted values, pred)
    scaled_center <- attr(xt_scaled, "scaled:center")
    scaled_scale <- attr(xt_scaled, "scaled:scale")
    xx <- x * scaled_scale[1] + scaled_center[1]
    yy <- y * scaled_scale[2] + scaled_center[2]
    
    # Usamos `xx`, `yy` y `z`
    p <- persp(xx, yy, z, theta = 30, phi = 30,
               col = "yellow",expand = 0.5, shade = 0.2,
               xlab = "Education", ylab = "Seniority", zlab = "Income")
    obs <-  trans3d(income_2$Education, income_2$Seniority, income_2$Income, p)
    pred <-  trans3d(income_2$Education, income_2$Seniority, fitted(interpolation_model), p)
    points(obs, col = "red", pch = 16)
    segments(obs$x, obs$y, pred$x, pred$y)
    

    Entradas relacionadas

    Referencias

    2018-09-15

    Gráficos de An Introduction to Statistical Learning - Figura 2.5.

    Gráfico a replicar

    Continuamos con la serie iniciada sobre la creación de gráficos del libro An Introduction to Statistical Learning. En esta ocasión replicaremos los gráficos de la figura 2.5. Utiliza el conjunto de datos Income2 (renta). El gráfico representa los valores observados de renta (en miles de dólares) en función de los años de educación y la antigüedad de 30 individuos. Los puntos rojos representan los valores observados para esas tres variables. La superficie amarilla representa un modelo de suavizado con thin-plate spline. En este caso . Las líneas verticales negras representan el error asociado con cada observación.

    Solución

    Prácticamente idéntica a las dos entradas anteriores, con la diferencia de que en lugar de la función loess o lm para ajustar un modelo lineal, empleamos otros modelos de suavizado con splines, con diferentes parámetros (tensor product, o bivariate).

  • Tensor product
  • Bivariate
  • income_2 <- read.csv("http://www-bcf.usc.edu/~gareth/ISL/Income2.csv")
    library(mgcv)
    model1 <- gam(Income ~ te(Education, Seniority, bs = "tp"), data = income_2) # tensor product
    model1 <- gam(Income ~ s(Education, Seniority, bs = "tp"), data = income_2) # bivariate
    x <- range(income_2$Education)
    x <- seq(x[1], x[2], length.out=30)
    y <- range(income_2$Seniority)
    y <- seq(y[1], y[2], length.out=30)
    z <- outer(x,y,
               function(Education,Seniority)
                         predict(model1, data.frame(Education,Seniority)))
    p <- persp(x,y,z, theta=30, phi=30,
               col="yellow",expand = 0.5,shade = 0.2,
               xlab="Education", ylab="Seniority", zlab="Income")
    obs <- trans3d(income_2$Education, income_2$Seniority,income_2$Income,p)
    pred <- trans3d(income_2$Education, income_2$Seniority,fitted(model1),p)
    points(obs, col="red",pch=16)
    segments(obs$x, obs$y, pred$x, pred$y)
    

    Entradas relacionadas

    Referencias

    Nube de datos