Modeling Offshore Wind Potential: OPG’s Toronto Islands Proposal

January 20XX

wind energy
renewable energy
data visualization
data analysis
simulation
energy modeling
Lake Ontario
Toronto Islands
turbine feasibility
air density
weather data
power generation
R (Programming Language)
Author

A. Srikanth

Published

August 29, 2025

Project Spotlight

Important

The work was originally implemented in Python and adapted for R.

Context

Here’s the setup: imagine 49 wind turbines arrayed in a 7-by-7 grid in Lake Ontario, just south of the Toronto Islands. I worked with a synthetic January weather record data set to drive the model and kept the turbine parameters grounded in utility-scale specs for our client, Ontario Power Generation. The weather is made up; the physics and unit handling are not. The aim is simple – estimate plausible output for a full array in winter and see how wind variability shows up in the numbers.

Code
suppressPackageStartupMessages({
  library(readr); library(dplyr); library(scales)
  library(plotly); library(htmltools); library(stringr)
})

bg_col      <- "#FAF8F1"
point_color <- "#751F2C"
font_family <- "Ramabhadra"

`%||%` <- function(x, y) if (is.null(x) || length(x) == 0) y else x

on_lat_range <- c(41.6, 56.9)
on_lon_range <- c(-95.2, -74.3)

csv_path <- "data/OntarioPowerGeneration_Ontario_Mapping.csv"
turb <- read_csv(csv_path, show_col_types = FALSE)

names(turb) <- tolower(trimws(names(turb)))
pick_col <- function(nms, candidates){ hit <- intersect(candidates, nms); if(length(hit)) hit[1] else NA_character_ }

col_lat  <- pick_col(names(turb), c("lat","latitude"))
col_lon  <- pick_col(names(turb), c("lon","long","longitude"))
col_name <- pick_col(names(turb), c("name","id","turbine"))
col_cap  <- pick_col(names(turb), c("capacity_mw","rated_mw","nameplate_mw","capacity","power_mw"))
col_stat <- pick_col(names(turb), c("status","type"))
if (is.na(col_lat) || is.na(col_lon)) stop("CSV must include latitude/longitude columns.")

turb_clean <- turb %>%
  transmute(
    lat    = as.numeric(.data[[col_lat]]),
    lon    = as.numeric(.data[[col_lon]]),
    name   = if (!is.na(col_name)) as.character(.data[[col_name]]) else NA_character_,
    cap    = if (!is.na(col_cap))  suppressWarnings(as.numeric(.data[[col_cap]])) else NA_real_,
    status = if (!is.na(col_stat)) as.character(.data[[col_stat]]) else NA_character_
  ) %>%
  filter(is.finite(lat), is.finite(lon)) %>%
  mutate(
    turb_no   = dplyr::row_number(),
    size_plot = if (any(is.finite(cap))) rescale(sqrt(pmax(cap, 0)), to = c(6, 22)) else 14,
    label_cap = ifelse(is.finite(cap), paste0(number(cap, accuracy = 0.1), " MW"), "N/A"),
    label_stat= ifelse(is.na(status) | status == "", "", paste0("<br>STATUS: ", htmltools::htmlEscape(status))),
    hover = paste0(
      "<span style='color:#CC6F77;'><b>TURBINE #", turb_no, "</b></span><br>",
      "<span style='color:slategray;'>LAT: ", round(lat, 4), " | LON: ", round(lon, 4), "</span>",
      label_stat
    )
  )

fig_on <- plot_ly(
  turb_clean,
  type = "scattergeo",
  mode = "markers",
  lat  = ~lat, lon = ~lon,
  text = ~hover, hoverinfo = "text",
  marker = list(size = ~size_plot, color = point_color, line = list(width = 0), opacity = 0.9,
                sizemode = "diameter", sizemin = 4)
) %>%
  layout(
    width = 790,
    font = list(family = font_family, size = 14),
    paper_bgcolor = bg_col, plot_bgcolor = bg_col,
    margin = list(l = 6, r = 6, t = 6, b = 6),
    hoverlabel = list(font = list(family = font_family, size = 14, color = "#313131"),
                      bgcolor = "#FFF", namelength = -1),
    hovermode = "closest",
    geo = list(
      scope = "north america",
      projection = list(type = "mercator"),
      # Explicitly show ALL of Ontario on load:
      lataxis = list(range = on_lat_range),
      lonaxis = list(range = on_lon_range),
      # Styling
      showland = TRUE, landcolor = "#FDEBED",
      subunitcolor = "#cccccc", countrycolor = "#cccccc",
      showsubunits = TRUE, showlakes = TRUE, lakecolor = "#e6f2ff",
      showcoastlines = TRUE, coastlinecolor = "slategray", coastlinewidth = 0.8,
      resolution = 50
    ),
    dragmode = "pan"
  ) %>%
  config(
    responsive = FALSE, scrollZoom = FALSE, doubleClick = TRUE,
    displayModeBar = TRUE, displaylogo = FALSE,
    modeBarButtonsToRemove = list("zoom2d","pan2d","select2d","lasso2d",
                                  "zoomIn2d","zoomOut2d","toggleSpikelines","toImage")
  )

browsable(
  tagList(
    tags$style(HTML("
      #on-map {
        border-radius: 6px;
        overflow-x: auto;    /* horizontal scroll when viewport is narrow */
        overflow-y: hidden;
        -webkit-overflow-scrolling: touch;
        background: #FAF8F1;
        padding: 0.75rem;
        width: 100%;
        max-width: 790px;
        margin: 0 auto;
      }
      #on-map .js-plotly-plot, #on-map .plot-container { min-width: 790px; }
      @media (max-width: 790px) { #on-map { padding: 0.5rem; } }
    ")),
    tags$div(id = "on-map", fig_on)
  )
)

Objectives

I wanted to answer three main questions:

  1. How much power could these turbines generate under typical January conditions?
  2. How does changing wind speed affect output?
  3. How do assumptions about the turbine itself (e.g., its limits and maximum capacity) affect the total?

Data Sources

The dataset included hourly values for wind speed, temperature, and air pressure. I also added turbine characteristics like cut-in speed (i.e., the minimum wind required to start turning), cut-out speed (i.e., the maximum safe wind before shutting down), and rated capacity (i.e., the maximum power the turbine can produce). Together, these formed the backbone of the calculations.

Data Preparation

The raw numbers weren’t in the right form to plug straight into equations. First, I calculated air density for each hour using temperature and pressure. Denser air carries more energy. Then I converted wind speeds from kilometres per hour (km/h) to metres per second (m/s). That step is critical because the turbine power formula uses SI units – the International System of Units common in science and engineering (metres, kilograms, seconds, etc.) – to keep calculations consistent.

\[ \rho = \frac{p_{\text{stn}} \times 1000}{287.05\,(T + 273.15)} \]

  • \(\rho\) — air density (\(\mathrm{kg\,m^{-3}}\))
  • \(p_{\text{stn}}\) — station pressure in \(\mathrm{kPa}\) (multiply by \(1000\) to convert to \(\mathrm{Pa}\))
  • \(T\) — air temperature in \(^{\circ}\mathrm{C}\) (add \(273.15\) to convert to \(\mathrm{K}\))
  • \(287.05\) — specific gas constant for dry air (\(\mathrm{J\,kg^{-1}\,K^{-1}}\))
Code
suppressPackageStartupMessages({
  library(readr)
  library(dplyr)
  library(scales)
  library(plotly)
  library(htmltools)
  library(DT)
})

bg_col      <- "#FAF8F1"
point_color <- "#751F2C"
font_family <- "Ramabhadra"

.fixed_fig_width  <- 1100L
.fixed_fig_height <- 390L
.wrapper_max_w    <- 790L

csv_path <- "data/OntarioPowerGeneration_IslandAirportWeather.csv"

column_names <- c(
  "month","day","time","temp","temp_flag","dew","dew_pt_temp_flag",
  "rel_hum","rel_hum_flag","wind_dir","wind_dir_flag","wind_spd_kmh",
  "wind_spd_kmh_flag","visibility","visibility_flag","stn_press",
  "stn_press_flag","hmdx","hmdx_flag","wind_chill","wind_chill_flag","weather"
)

wx_raw <- readr::read_csv(
  file = csv_path,
  skip = 1,
  col_names = column_names,
  na = c("", "NA"),
  show_col_types = FALSE
)

safe_num <- function(x) suppressWarnings(as.numeric(x))

wx <- wx_raw %>%
  dplyr::mutate(
    temp          = safe_num(temp),
    stn_press     = safe_num(stn_press),
    wind_spd_kmh  = safe_num(wind_spd_kmh),
    rel_hum       = safe_num(rel_hum),
    visibility    = safe_num(visibility),

    air_density_kg_m3 = dplyr::if_else(
      !is.na(temp) & !is.na(stn_press),
      (stn_press * 1000) / (287.05 * (temp + 273.15)),
      NA_real_
    ),

    hover = paste0(
      "<span style='color:#CC6F77;'><b>AIR DENSITY POINT</b></span><br>",
      "TEMP: ", ifelse(is.na(temp), "NA", sprintf("%.1f", temp)), " °C<br>",
      "PRESSURE: ",
      ifelse(is.na(stn_press), "NA", scales::number(stn_press, accuracy = 0.01)),
      " kPa<br>",
      "DENSITY: ",
      ifelse(is.na(air_density_kg_m3), "NA", sprintf("%.2f", air_density_kg_m3)),
      " kg/m³<br>",
      "<span style='color:slategray;'>",
      "WIND: ", ifelse(is.na(wind_spd_kmh), "NA", paste0(round(wind_spd_kmh), " km/h")),
      " | HUMIDITY: ", ifelse(is.na(rel_hum), "NA", paste0(round(rel_hum), "%")),
      " | VIS: ", ifelse(is.na(visibility), "NA", paste0(round(visibility, 1), " km")),
      ifelse(is.na(weather) | weather == "", "", paste0(" | WX: ", weather)),
      "</span>"
    )
  )

fig <- plotly::plot_ly(
  wx, type = "scatter", mode = "markers",
  x = ~temp, y = ~air_density_kg_m3, text = ~hover, hoverinfo = "text",
  marker = list(size = 9, color = point_color, line = list(width = 0), opacity = 0.9,
                sizemode = "diameter", sizemin = 4)
) %>%
  layout(
    width  = .fixed_fig_width,
    height = .fixed_fig_height,
    font = list(family = font_family, size = 12),
    paper_bgcolor = bg_col, plot_bgcolor = bg_col,
    margin = list(l = 18, r = 18, t = 36, b = 36),
    legend = list(orientation = "v", x = 1.05, y = 1, xanchor = "left", yanchor = "top",
                  font = list(size = 16), title = list(text = "")),
    xaxis = list(title = list(text = "TEMPERATURE (°C)", standoff = 20),
                 showspikes = TRUE, tickfont = list(size = 16), gridcolor = "#E8E8E8",
                 zeroline = FALSE, fixedrange = TRUE),
    yaxis = list(title = list(text = "AIR DENSITY (kg/m<sup>3</sup>)", standoff = 20),
                 showspikes = TRUE, tickfont = list(size = 16), automargin = TRUE, fixedrange = TRUE),
    uniformtext = list(minsize = 16),
    hovermode = "y unified",
    hoverlabel = list(
      font = list(family = font_family, size = 14, color = "#313131"),
      bgcolor = "#FFF",
      namelength = -1
    ),
    dragmode = FALSE
  ) %>%
  config(
    responsive = FALSE,
    scrollZoom = TRUE, doubleClick = FALSE,
    modeBarButtonsToRemove = list(
      "zoom2d","pan2d","select2d","lasso2d","zoomIn2d","zoomOut2d",
      "autoScale2d","resetScale2d","toggleSpikelines","toImage"
    ),
    displaylogo = FALSE, displayModeBar = TRUE, showTips = FALSE
  )

htmltools::browsable(
  htmltools::tagList(
    htmltools::tags$style(htmltools::HTML(sprintf("
      /* Scrollable wrapper for the plot */
      #temp-plot-wrap {
        border-radius: 6px;
        overflow: hidden;                /* clip vertical overflow */
        overflow-x: auto;                /* horizontal scroll when narrow */
        -webkit-overflow-scrolling: touch;
        background: %s;
        padding: 0.75rem;
        width: 100%%;
        max-width: %dpx;                 /* cap container width to create overflow */
        margin: 0 auto 1rem;
      }
      /* Ensure inner plot area is at least the fixed figure width */
      #temp-plot-wrap > div { min-width: %dpx; }

      /* Scrollable wrapper for the table (optional but consistent) */
      #temp-table-wrap {
        border-radius: 6px;
        overflow: hidden;
        overflow-x: auto;
        -webkit-overflow-scrolling: touch;
        background: %s;
        padding: 0.75rem;
        width: 100%%;
        max-width: %dpx;
        margin: 0 auto 1rem;
      }
      #temp-table-wrap > div { min-width: %dpx; }

      @media (max-width: %dpx) {
        #temp-plot-wrap, #temp-table-wrap { padding: 0.5rem; }
      }

      /* DataTable alignment */
      #temp-table table.dataTable thead th { text-align: left !important; }
      #temp-table table.dataTable tbody td { text-align: right !important; }
      #temp-table table.dataTable thead th:first-child,
      #temp-table table.dataTable tbody td:first-child { text-align: left !important; }
    ",
      bg_col, .wrapper_max_w, .fixed_fig_width,
      bg_col, .wrapper_max_w, .fixed_fig_width,
      .wrapper_max_w
    ))),
    htmltools::tags$div(id = "temp-plot-wrap", fig),

    htmltools::tags$div(
      id = "temp-table-wrap",
      htmltools::tags$div(
        id = "temp-table",
        DT::datatable(
          wx %>%
            dplyr::transmute(
              `TEMPERATURE (°C)`        = temp,
              `STATION PRESSURE (kPa)`  = stn_press,
              `AIR DENSITY (kg/m³)`     = ifelse(is.na(air_density_kg_m3), NA, sprintf("%.2f", air_density_kg_m3))
            ) %>%
            dplyr::slice(1:9),
          rownames = TRUE,
          options = list(
            autoWidth = TRUE,
            dom = "tip",
            scrollY = "135px",
            scrollCollapse = TRUE,
            paging = FALSE
          )
        )
      )
    )
  )
)
744 rows × 24 columns


The scatterplot of air density versus temperature (above) shows the expected trend: colder air is denser, and as temperatures rise, density drops.

Methodology

For every hour of data, I applied the turbine power equation. If the wind was too weak (below cut-in) or too strong (above cut-out), the power was set to zero. When the wind was in range, power increased with the cube of wind speed but was capped once it hit the turbine’s rated capacity. The outputs were calculated in watts, then scaled up to megawatts for clarity.

\[ P(V)= \begin{cases} 0, & V < V_{ci} \\ \min\!\left(\tfrac{1}{2}\,\rho\,A\,C_p\,V^3,\; P_{\text{rated}}\right), & V_{ci} \le V \le V_{co} \\ 0, & V > V_{co} \end{cases} \]

  • \(V\) — wind speed (\(\mathrm{m\,s^{-1}}\))
  • \(\rho\) — air density (\(\mathrm{kg\,m^{-3}}\))
  • \(A\) — rotor-swept area (\(\mathrm{m^{2}}\))
  • \(C_p\) — power coefficient
  • \(V_{ci}, V_{co}\) — cut-in and cut-out speeds
  • \(P_{\mathrm{rated}}\) — turbine nameplate power

Analysis

Looking at the hourly results, you can see how much variability comes from wind alone. A calm day contributes almost nothing, while a windy spell produces spikes of several thousand megawatt-hours across the farm. By multiplying one turbine’s output by 49, the project captured the effect of running a full array instead of a single unit.

The bar-and-line chart shows this relationship in action. Daily energy production (bars) rises and falls alongside the wind speed trend (line). It’s not a groundbreaking observation—higher wind means more power—but it’s still satisfying to see the data confirm the intuition. In total, the farm would have produced about 43,078 MWh in January, enough to power thousands of homes for a month!


Code
library(readr)
library(dplyr)
library(DT)
library(htmltools)

weather <- read_csv("data/OntarioPowerGeneration_IslandAirportWeather.csv")

weather <- weather %>%
  mutate(`wind_spd_ms` = round(`wind_spd_kmh` * 0.277778, 2))

browsable(
  tagList(
    tags$style(HTML("
      #temp-plot, #wind-table {
        border-radius: 6px;
        overflow: hidden;
        background: #FAF8F1;
        padding: 0.75rem;
        width: 100%;
        max-width: 900px;
        margin: 0 auto 1rem;
      }
      @media (max-width: 790px) {
        #temp-plot, #wind-table { padding: 0.5rem; }
      }
      /* Header left aligned */
      #wind-table table.dataTable thead th {
        text-align: left !important;
      }
      /* Body cells right aligned */
      #wind-table table.dataTable tbody td {
        text-align: right !important;
      }
      /* Row numbers (first col) left aligned */
      #wind-table table.dataTable thead th:first-child,
      #wind-table table.dataTable tbody td:first-child {
        text-align: left !important;
      }
    ")),
    tags$div(
      id = "wind-table",
      DT::datatable(
        weather %>%
          select(
            "WIND SPEED (km/h)" = wind_spd_kmh,
            "WIND SPEED (m/s)"  = wind_spd_ms
          ) %>%
          dplyr::slice(1:9),
        rownames = TRUE,
        options = list(
          autoWidth = TRUE,
          dom = "tip",
          scrollY = "135px",
          scrollCollapse = TRUE,
          paging = FALSE
        )
      )
    )
  )
)
744 rows × 23 columns


Code
suppressPackageStartupMessages({
  library(readr)
  library(dplyr)
  library(htmltools)
  library(DT)
  library(stringr)
})

wx_path    <- "data/OntarioPowerGeneration_IslandAirportWeather.csv"
specs_path <- "data/OntarioPowerGeneration_WindTurbineSpecifications.csv"

safe_num <- function(x) suppressWarnings(as.numeric(x))

parse_num <- function(x) {
  x <- as.character(x %||% NA_character_)
  x <- gsub("[^0-9eE+\\-\\.]", "", x)
  out <- suppressWarnings(as.numeric(x))
  if (length(out) == 0) NA_real_ else out
}

find_spec <- function(specs, pattern) {
  hit <- specs %>%
    filter(str_detect(tolower(field), tolower(pattern))) %>%
    slice(1) %>%
    pull(values)
  if (length(hit) == 0) return(NA_real_)
  parse_num(hit[1])
}

`%||%` <- function(x,y) if (is.null(x) || length(x)==0) y else x

wx_cols <- c(
  "month","day","time","temp","temp_flag","dew","dew_pt_temp_flag",
  "rel_hum","rel_hum_flag","wind_dir","wind_dir_flag","wind_spd_kmh",
  "wind_spd_kmh_flag","visibility","visibility_flag","stn_press",
  "stn_press_flag","hmdx","hmdx_flag","wind_chill","wind_chill_flag","weather"
)

wx <- readr::read_csv(
  file = wx_path,
  col_names = wx_cols,
  skip = 1,
  na = c("", "NA"),
  show_col_types = FALSE
) %>%
  mutate(
    temp         = safe_num(temp),
    stn_press    = safe_num(stn_press),
    wind_spd_kmh = safe_num(wind_spd_kmh),
    `Wind Spd (m/s)` = ifelse(is.finite(wind_spd_kmh), wind_spd_kmh * 0.277778, NA_real_),
    `air_density (kg/m^3)` = ifelse(
      is.finite(temp) & is.finite(stn_press),
      (stn_press * 1000) / (287.05 * (temp + 273.15)),
      NA_real_
    )
  )

specs <- readr::read_csv(specs_path, show_col_types = FALSE)
names(specs) <- tolower(names(specs))

cut_in       <- find_spec(specs, "cut[- ]?in windspeed|cut[- ]?in")
cut_out      <- find_spec(specs, "cut[- ]?out windspeed|cut[- ]?out")
cp           <- find_spec(specs, "maximum power coefficient|\\bcp\\b")
area_m2      <- find_spec(specs, "turbine swept area|swept area|rotor area")
p_nom_mw     <- find_spec(specs, "turbine nominal power|rated power|nameplate")
p_nom_watt   <- ifelse(is.finite(p_nom_mw), p_nom_mw * 1e6, NA_real_)

if (!is.finite(cut_in))     cut_in     <- -Inf
if (!is.finite(cut_out))    cut_out    <-  Inf
if (!is.finite(cp))         cp         <- NA_real_
if (!is.finite(area_m2))    area_m2    <- NA_real_
if (!is.finite(p_nom_watt)) p_nom_watt <- NA_real_

wx_power <- wx %>%
  mutate(
    .valid_ws  = is.finite(`Wind Spd (m/s)`) &
                 `Wind Spd (m/s)` >= cut_in &
                 `Wind Spd (m/s)` <= cut_out,
    .raw_power = ifelse(
      .valid_ws & is.finite(`air_density (kg/m^3)`) & is.finite(cp) & is.finite(area_m2),
      `air_density (kg/m^3)` * 0.5 * cp * area_m2 * (`Wind Spd (m/s)`^3),
      0
    ),
    `Power (W)` = ifelse(
      .valid_ws & is.finite(p_nom_watt),
      pmin(.raw_power, p_nom_watt),
      ifelse(.valid_ws, .raw_power, 0)
    )
  ) %>%
  select(`Wind Spd (m/s)`, `air_density (kg/m^3)`, `Power (W)`)

wx_power_fmt <- wx_power %>%
  mutate(
    `WIND SPEED (m/s)`              = ifelse(is.na(`Wind Spd (m/s)`), NA, sprintf("%.2f", `Wind Spd (m/s)`)),
    `AIR DENSITY (kg/m<sup>3</sup>)` = ifelse(is.na(`air_density (kg/m^3)`), NA, sprintf("%.2f", `air_density (kg/m^3)`)),
    `POWER (W)`                     = ifelse(is.na(`Power (W)`), NA, format(`Power (W)`, scientific = TRUE, digits = 3))
  ) %>%
  select(`WIND SPEED (m/s)`, `AIR DENSITY (kg/m<sup>3</sup>)`, `POWER (W)`)

browsable(
  tagList(
    tags$style(HTML("
      #wind-power-wrap {
        border-radius: 6px;
        overflow: hidden;
        background: #FAF8F1;
        padding: 0.75rem;
        width: 100%;
        max-width: 900px;
        margin: 0 auto 1rem;
      }
      #wind-power-table table.dataTable thead th { text-align: left !important; }
      #wind-power-table table.dataTable tbody td { text-align: right !important; }
      #wind-power-table table.dataTable thead th:first-child,
      #wind-power-table table.dataTable tbody td:first-child { text-align: left !important; }
    ")),
    tags$div(
      id = "wind-power-wrap",
      tags$div(
        id = "wind-power-table",
        # ---- Display table ----
        DT::datatable(
          wx_power_fmt %>% dplyr::slice(1:9),
          rownames = TRUE,
          escape = FALSE,
          options = list(
            autoWidth = TRUE,
            dom = "tip",
            scrollY = "135px",
            scrollCollapse = TRUE,
            paging = FALSE
          )
        )
      )
    )
  )
)
744 rows × 3 columns


Code
wx_power_all <- wx_power %>%
  mutate(
    `Power_All_Turbines_MW` = (`Power (W)` * 49) / 1e6
  )

total_mwh <- sum(wx_power_all$Power_All_Turbines_MW, na.rm = TRUE)
cat(sprintf("Total energy produced by the windfarm in January: %.2f MWh\n", total_mwh))
Total energy produced by the windfarm in January: 43077.82 MWh


Code
suppressPackageStartupMessages({
  library(readr); library(dplyr); library(stringr)
  library(lubridate); library(plotly); library(htmltools)
})

wx_path    <- "data/OntarioPowerGeneration_IslandAirportWeather.csv"
specs_path <- "data/OntarioPowerGeneration_WindTurbineSpecifications.csv"

safe_num  <- function(x) suppressWarnings(as.numeric(x))
parse_num <- function(x){
  x <- gsub("[^0-9eE+\\-\\.]", "", as.character(x))
  suppressWarnings(as.numeric(x))
}
find_spec <- function(specs, pattern){
  hit <- specs %>% dplyr::filter(stringr::str_detect(tolower(field), tolower(pattern))) %>%
    dplyr::slice(1) %>% dplyr::pull(values)
  if (!length(hit)) return(NA_real_)
  parse_num(hit[1])
}

wx_cols <- c(
  "month","day","time","temp","temp_flag","dew","dew_pt_temp_flag",
  "rel_hum","rel_hum_flag","wind_dir","wind_dir_flag","wind_spd_kmh",
  "wind_spd_kmh_flag","visibility","visibility_flag","stn_press",
  "stn_press_flag","hmdx","hmdx_flag","wind_chill","wind_chill_flag","weather"
)

wx <- readr::read_csv(
  file = wx_path, col_names = wx_cols, skip = 1, na = c("", "NA"), show_col_types = FALSE
) %>%
  dplyr::mutate(
    temp          = safe_num(temp),
    stn_press     = safe_num(stn_press),
    wind_spd_kmh  = safe_num(wind_spd_kmh),
    `Wind Spd (m/s)` = dplyr::if_else(is.finite(wind_spd_kmh), wind_spd_kmh * 0.277778, NA_real_),
    `air_density (kg/m^3)` = dplyr::if_else(
      is.finite(temp) & is.finite(stn_press),
      (stn_press * 1000) / (287.05 * (temp + 273.15)),
      NA_real_
    )
  )

specs <- readr::read_csv(specs_path, show_col_types = FALSE)
names(specs) <- tolower(names(specs))

cut_in     <- find_spec(specs, "cut[- ]?in windspeed|cut[- ]?in")
cut_out    <- find_spec(specs, "cut[- ]?out windspeed|cut[- ]?out")
cp         <- find_spec(specs, "maximum power coefficient|\\bcp\\b")
area_m2    <- find_spec(specs, "turbine swept area|rotor area|swept area")
p_nom_mw   <- find_spec(specs, "turbine nominal power|rated power|nameplate")
p_nom_watt <- ifelse(is.finite(p_nom_mw), p_nom_mw * 1e6, NA_real_)
if (!is.finite(cut_in))  cut_in  <- -Inf
if (!is.finite(cut_out)) cut_out <-  Inf

wx <- wx %>%
  dplyr::mutate(
    .valid_ws  = is.finite(`Wind Spd (m/s)`) &
                 `Wind Spd (m/s)` >= cut_in &
                 `Wind Spd (m/s)` <= cut_out,
    .raw_power = dplyr::if_else(
      .valid_ws & is.finite(`air_density (kg/m^3)`) & is.finite(cp) & is.finite(area_m2),
      `air_density (kg/m^3)` * 0.5 * cp * area_m2 * (`Wind Spd (m/s)`^3),
      0
    ),
    `Power (W)` = dplyr::if_else(
      .valid_ws & is.finite(p_nom_watt),
      pmin(.raw_power, p_nom_watt),
      dplyr::if_else(.valid_ws, .raw_power, 0)
    ),
    Power_All_Turbines_MW = (`Power (W)` * 49) / 1e6
  )

wx_daily <- wx %>%
  dplyr::mutate(date = lubridate::make_date(2025, month, day)) %>%
  dplyr::filter(month == 1) %>%
  dplyr::group_by(date) %>%
  dplyr::summarise(
    Daily_MWh = sum(Power_All_Turbines_MW, na.rm = TRUE),
    Daily_Wind_kmh = mean(wind_spd_kmh, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  dplyr::arrange(date)

bg_col      <- "#FAF8F1"
point_color <- "#751F2C"
font_family <- "Ramabhadra"

.fixed_fig_width <- 1100L
.fixed_fig_height <- 390L

fig_jan <- plotly::plot_ly(
  wx_daily,
  x = ~date, y = ~Daily_MWh,
  type = "bar",
  name = "Power (MWh)",
  marker = list(color = point_color),
  hovertemplate = "POWER: %{y:.1f}<extra></extra>"
) %>%
  add_trace(
    x = ~date, y = ~Daily_Wind_kmh,
    type = "scatter", mode = "lines+markers",
    name = "Wind Speed (km/h)",
    yaxis = "y2",
    line = list(width = 3),
    marker = list(size = 6),
    hovertemplate = "WIND SPEED: %{y:.1f} km/h<extra></extra>"
  ) %>%
  layout(
    width  = .fixed_fig_width,
    height = .fixed_fig_height,
    font = list(family = font_family, size = 12),
    paper_bgcolor = bg_col,
    plot_bgcolor  = bg_col,
    margin = list(l = 18, r = 18, t = 36, b = 36),
    showlegend = TRUE,
    legend = list(orientation = "v", x = 1.05, y = 1, xanchor = "left", yanchor = "bottom", font = list(size = 16), title = list(text = "")),
    xaxis = list(
      title = list(text = "DATE", standoff = 20),
      tickformat = "%b %d",
      showspikes = TRUE,
      tickfont = list(size = 16),
      gridcolor = "#E8E8E8",
      zeroline = FALSE,
      fixedrange = TRUE
    ),
    yaxis = list(
      title = list(text = "POWER (MWh)", standoff = 20),
      showspikes = TRUE,
      tickfont = list(size = 16),
      automargin = TRUE,
      fixedrange = TRUE
    ),
    yaxis2 = list(
      title = list(text = "WIND (km/h)", standoff = 20),
      overlaying = "y",
      side = "right",
      showgrid = FALSE,
      tickfont = list(size = 16),
      fixedrange = TRUE
    ),
    uniformtext = list(minsize = 16),
    hovermode = "x unified",
    hoverlabel = list(
      font = list(family = font_family, size = 14, color = "#313131"),
      bgcolor = "#FFF",
      namelength = -1
    ),
    dragmode = FALSE
  ) %>%
  config(
    responsive = FALSE,
    scrollZoom = FALSE, doubleClick = FALSE,
    modeBarButtonsToRemove = list(
      "zoom2d","pan2d","select2d","lasso2d",
      "zoomIn2d","zoomOut2d","autoScale2d","resetScale2d",
      "toggleSpikelines","toImage"
    ),
    displaylogo = FALSE, displayModeBar = TRUE, showTips = FALSE
  )

htmltools::browsable(
  htmltools::tagList(
    htmltools::tags$style(htmltools::HTML(sprintf("
      #jan-plot-wrap {
        border-radius: 6px;
        overflow: hidden;
        overflow-x: auto;
        -webkit-overflow-scrolling: touch;
        background: %s;
        padding: 0.75rem;
        width: 100%%;
        max-width: 790px;
        margin: 0 auto 1rem;
      }
      #jan-plot-wrap > div {
        min-width: %dpx;
      }
      @media (max-width: 790px) {
        #jan-plot-wrap { padding: 0.5rem; }
      }
    ", bg_col, .fixed_fig_width))),
    htmltools::tags$div(id = "jan-plot-wrap", fig_jan)
  )
)

Results & Next Steps

The results underline two points. First, wind is highly variable, even within a single month, and that variability drives power output. Second, scaling up to dozens of turbines smooths out some of that variability, producing a meaningful and steady supply over time.

If this were a real feasibility study, the next steps would involve testing across multiple months or years, adding turbine placement effects, and layering in grid integration considerations. Even with a single month of simulated data, the exercise shows how basic physics, weather data, and engineering limits combine to give a picture of what a wind project might deliver.