{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "da0a22b9-5a5d-4bd1-8e03-166029a08a69",
   "metadata": {},
   "source": [
    "# Instalar bibliotecas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a3aa15a4-36ea-42da-8744-3dfb2134abbc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# ! pip install windpowerlib"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0218dcc0-bbe1-4d59-b19b-64b71f4c8d7d",
   "metadata": {},
   "source": [
    "# Importar bibliotecas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58e01fca-5fbf-4570-9b96-66659307d56a",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np \n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.dates as mdates\n",
    "import requests\n",
    "\n",
    "from scipy.stats import weibull_min"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84a637e9-7833-418f-aae3-da8bfdb9e567",
   "metadata": {},
   "source": [
    "# Exercícios"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03344c46-d361-40ce-9bd9-372e823bae17",
   "metadata": {},
   "source": [
    "A distribuição de Weibull é usada na caracterização do recurso eólico. Esta distribuição descreve a frequência de ocorrência de cada velocidade, e é ajustada ao recurso de cada local através dos seus parâmetros $c$ e $k$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d163b205-6ecd-4e64-a8f1-d346d16b0eaa",
   "metadata": {},
   "source": [
    "### **1.** Para uma distribuição de Weibull com fator de escala $c=8$ m/s e fator de forma $k=1.6$ m/s, determina a distribuição de probabilidade da velocidade do vento, $f(u)$, e a distribuição de probabilidade acumulada, $F(u)$, para $u$ entre 0 e 25 m/s, em intervalos de 1 m/s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b9b575f8-605a-470b-be9a-7eff9ee25d80",
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 1.6\n",
    "c = 8  \n",
    "\n",
    "# valores discretos da velocidade do vento\n",
    "u = np.arange(0, 26, 1)\n",
    "\n",
    "df = pd.DataFrame(index=u)\n",
    "df.index.name = \"Velocidade do vento (m/s)\"\n",
    "\n",
    "# Distribuição de probabilidade\n",
    "df['f(u)'] = # TODO: determina f(u) e F(u)\n",
    "df['F(u)'] = \n",
    "\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a259e5b5-064a-4e3d-8acb-7c2c35ec29e8",
   "metadata": {},
   "source": [
    "#### a) Visualiza f(u) e F(u) graficamente, em separado."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e39a12e-bd33-4ab8-be09-4b92aefcfe5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot\n",
    "plt.figure(figsize=(10,5))\n",
    "plt.subplot(1,2,1)\n",
    "plt.plot(df['f(u)'], label=\"f(u)\")\n",
    "plt.xlim([0,25])\n",
    "plt.ylim(bottom=0)\n",
    "plt.xlabel(\"Velocidade do vento (m/s)\")\n",
    "plt.ylabel(\"Densidade de probabilidade\")\n",
    "plt.legend()\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.plot(df['F(u)'], label=\"F(u)\", color='g')\n",
    "plt.xlim([0,25])\n",
    "plt.ylim(bottom=0)\n",
    "plt.xlabel(\"Velocidade do vento (m/s)\")\n",
    "plt.ylabel(\"Probabilidade cumulativa\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf7d5966-5577-4af9-8573-6a101f521d54",
   "metadata": {},
   "source": [
    "#### b) Qual a moda da velocidade? E a média?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48e49364-4a31-4e73-9e8e-0f703fdfce06",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: determina a moda e a média\n",
    "moda = \n",
    "media = \n",
    "\n",
    "print(f\"Moda: {moda:.1f} m/s\")\n",
    "print(f\"Média: {media:.1f} m/s\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "65b09033-7699-48d4-805b-eee66a80bb19",
   "metadata": {},
   "source": [
    "## Parque Eólico de Marvila"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41d1b587-9970-4d08-b446-9f82edce1319",
   "metadata": {},
   "source": [
    "O [Parque Eólico de Marvila](https://www.openstreetmap.org/relation/14047195#map=13/39.57169/-8.65354) localiza-se perto de Mira de Aire, em Leiria. Tem uma potência instalada de 12 MW, com 6 turbinas Senvion MM92. Tem as coordenadas aproximadas:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a62aa903-22bd-447f-a227-43420dd22836",
   "metadata": {},
   "outputs": [],
   "source": [
    "latitude = 39.559912\n",
    "longitude = -8.697052"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed5987dc-1e15-438a-b020-4b2e07596e5c",
   "metadata": {},
   "source": [
    "Através do Open-Meteo, conseguimos dados horários de velocidade do vento a 10 e 100 m de altura para a localização deste parque eólico: tanto dados históricos como previsões. A imagem seguinte apresenta as previsões para a próxima semana."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a458fd7-c4cb-4ced-8ea0-2caaa28320fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "def km_h_to_m_s(data):\n",
    "    \"\"\"Converter km/h em m/s\"\"\"\n",
    "    return data*1000/3600"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2dbf76ef-1209-46fb-80ab-1f212a6edc7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open-Meteo API URL (Get hourly wind speeds at 10m & 100m)\n",
    "url = f\"https://api.open-meteo.com/v1/forecast?latitude={latitude}&longitude={longitude}&hourly=wind_speed_10m,wind_speed_100m&timezone=auto\"\n",
    "\n",
    "# Fetch data\n",
    "response = requests.get(url)\n",
    "data = response.json()\n",
    "\n",
    "# Extract wind speed data\n",
    "wind_speed_10m = data[\"hourly\"][\"wind_speed_10m\"]\n",
    "wind_speed_100m = data[\"hourly\"][\"wind_speed_100m\"]\n",
    "timestamps = pd.to_datetime(data[\"hourly\"][\"time\"])\n",
    "\n",
    "# Convert to Pandas DataFrame\n",
    "df = pd.DataFrame({\"Time\": timestamps, \"Wind Speed 10m\": wind_speed_10m, \"Wind Speed 100m\": wind_speed_100m})\n",
    "df.set_index(\"Time\", inplace=True)\n",
    "df = km_h_to_m_s(df) # converter km/h em m/s\n",
    "\n",
    "# Plot wind speed\n",
    "plt.figure(figsize=(10,5))\n",
    "plt.plot(df.index, df[\"Wind Speed 10m\"], label=\"10m de altura\", color=\"blue\")\n",
    "plt.plot(df.index, df[\"Wind Speed 100m\"], label=\"100m de altura\", color=\"red\", linestyle=\"dashed\")\n",
    "plt.xlabel(\"Data\")\n",
    "plt.ylabel(\"Velocidade do vento (m/s)\")\n",
    "plt.xlim([df.index.min(), df.index.max()])\n",
    "plt.ylim(bottom=0)\n",
    "plt.title(\"Previsão de velocidades do vento a 10m e 100m from Open-Meteo API\")\n",
    "plt.legend()\n",
    "plt.xticks(rotation=45)\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6c34bd6-3fc7-484c-b950-fda2c945e176",
   "metadata": {},
   "source": [
    "### **2.** Usa o Open-Meteo para obter um ano completo de dados históricos de velocidade do vento a 100 m de altura para a localização do Parque Eólico de Marvila."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cccff6ad-b7e5-4fc6-b65f-f24d01904df5",
   "metadata": {},
   "source": [
    "Vamos agora buscar dados históricos para esta localização para 10 anos (2014-2023) e para 100 m de altura."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7bebf514-85b0-4694-b87d-23d29304d920",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define date range for historical data (1 year)\n",
    "start_date = \"2014-01-01\"\n",
    "end_date = \"2023-12-31\"\n",
    "\n",
    "# Open-Meteo Archive API for past wind speeds\n",
    "url = f\"https://archive-api.open-meteo.com/v1/archive?latitude={latitude}&longitude={longitude}&start_date={start_date}&end_date={end_date}&hourly=wind_speed_100m&timezone=auto\"\n",
    "\n",
    "# Fetch data\n",
    "response = requests.get(url)\n",
    "data = response.json()\n",
    "\n",
    "# Extract wind speed data\n",
    "wind_speed_100m = data[\"hourly\"][\"wind_speed_100m\"]\n",
    "timestamps = pd.to_datetime(data[\"hourly\"][\"time\"])\n",
    "\n",
    "# Convert to Pandas DataFrame\n",
    "df_historico = pd.DataFrame({\"Time\": timestamps, \"Velocidade do vento 100m\": wind_speed_100m})\n",
    "df_historico.set_index(\"Time\", inplace=True)\n",
    "df_historico = km_h_to_m_s(df_historico)\n",
    "\n",
    "df_historico"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8dee825-61dc-4297-8194-1358e62d4793",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot historical wind speed\n",
    "plt.figure(figsize=(12,5))\n",
    "plt.plot(df_historico.index, df_historico[\"Velocidade do vento 100m\"])\n",
    "plt.xlabel(\"Data\")\n",
    "plt.ylabel(\"Velocidade do vento (m/s)\")\n",
    "plt.xticks(rotation=45)\n",
    "\n",
    "plt.xlim(df_historico.index.min(), df_historico.index.max())\n",
    "plt.ylim(bottom=0)\n",
    "\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dba6cf8f-0a08-43ed-9dce-8e532cfbebf0",
   "metadata": {},
   "source": [
    "A partir daqui podemos construir a distribuição de probabilidade da velocidade do vento."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b553658b-bf9b-4ff0-aae1-3a40397271fb",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Step 1: Extract wind speed data at 100m\n",
    "wind_speed_100m = df_historico[\"Velocidade do vento 100m\"]\n",
    "\n",
    "# Step 2: Define bins (1 m/s intervals from 0 to max wind speed)\n",
    "bins = np.arange(0, np.max(wind_speed_100m) + 1, 1)  # 1 m/s step\n",
    "\n",
    "# Step 3: Compute histogram (frequency count in each bin)\n",
    "hist, bin_edges = np.histogram(wind_speed_100m, bins=bins, density=True)\n",
    "\n",
    "df = pd.DataFrame({\n",
    "    \"Velocidade do vento (m/s)\": bin_edges[:-1],  # Exclude the last bin edge (it's the upper limit of the last bin)\n",
    "    \"f(u)\": hist\n",
    "}).set_index(\"Velocidade do vento (m/s)\")\n",
    "\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ea68361-7292-4677-8082-5bf45bbe2a3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,5))\n",
    "plt.bar(df.index, df['f(u)'], width=1.0, edgecolor=\"black\", alpha=0.7, color=\"skyblue\")\n",
    "\n",
    "plt.xlabel(\"Velocidade do vento a 100m (m/s)\")\n",
    "plt.ylabel(\"Frequência\")\n",
    "plt.grid(axis=\"y\", linestyle=\"--\", alpha=0.7)\n",
    "\n",
    "plt.xlim([-1, bin_edges[-1]])\n",
    "\n",
    "# Show plot\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c5d07826-8257-4b4c-93be-0c457e347966",
   "metadata": {},
   "source": [
    "A forma é semelhante à de uma distribuição de Weibull, certo? Vamos então descobrir quais os parâmetros da distribuição de Weibull mais adequados."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9f97923-eed0-4977-941d-a9325a6a9533",
   "metadata": {},
   "outputs": [],
   "source": [
    "velocidade_vento_100m = df_historico['Velocidade do vento 100m']\n",
    "\n",
    "# Step 2: Fit Weibull distribution using MLE\n",
    "shape, loc, scale = weibull_min.fit(velocidade_vento_100m, floc=0)\n",
    "k, c = shape, scale  # Weibull shape (k) and scale (c) parameters\n",
    "\n",
    "# Step 4: Compute Weibull PDF for overlay\n",
    "# x = np.linspace(0, np.max(wind_speed_100m), 100)  # Continuous range for smooth curve\n",
    "x = np.arange(0, np.max(wind_speed_100m), 1)\n",
    "pdf = weibull_min.pdf(x, k, scale=c)  # Weibull probability density function\n",
    "\n",
    "print(f\"Parâmetros estimados para a distribuição de Weibull: k = {k:.1f}, c = {c:.1f}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec5fcbef-894a-4545-9c61-8477f6384ce9",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,5))\n",
    "plt.bar(df.index, df['f(u)'], width=1.0, edgecolor=\"black\", alpha=0.7, color=\"skyblue\", label=\"Dados históricos\")\n",
    "plt.plot(x, pdf, \"r-\", lw=2, label=f\"Distribuição de Weibull ajustada (k={k:.1f}, c={c:.1f})\")\n",
    "\n",
    "plt.xlabel(\"Velocidade do vento a 100m (m/s)\")\n",
    "plt.ylabel(\"Frequência\")\n",
    "plt.grid(axis=\"y\", linestyle=\"--\", alpha=0.7)\n",
    "\n",
    "plt.xlim([-1, bin_edges[-1]])\n",
    "\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "85728c58-dba1-42ea-979c-bc3784a312e0",
   "metadata": {},
   "source": [
    "### **3.** Calcula a densidade de potência do vento, em kW/m2, para cada velocidade, e apresenta os resultados num gráfico."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5b3f28ec-eaeb-46c3-b116-5ae8a962c32b",
   "metadata": {},
   "outputs": [],
   "source": [
    "df['Densidade de potência do vento (kW/m2)'] = # TODO: calcular densidade de potência"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "919443d0-c7aa-4cb1-898b-539addfc4774",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: gráfico densidade de potência"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b305c03-21ff-4fb6-894a-9d02a536b3b6",
   "metadata": {},
   "source": [
    "### **4.** Compara a densidade de potência média com a densidade de potência à velocidade média. O que podes concluir?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6fd96cd0-7d75-4e57-be55-330521dae567",
   "metadata": {},
   "outputs": [],
   "source": [
    "dens_potencia_media = # TODO: calcula a densidade de potência média\n",
    "print(f\"Densidade de potência média: {dens_potencia_media:.3f} kW/m2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "edb0a548-7b89-489e-bc68-323d63d1954b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: calcula a velocidade média e a respetiva densidade de potência\n",
    "velocidade_media = \n",
    "dens_potencia = \n",
    "\n",
    "print(f\"Velocidade média: {velocidade_media:.1f} m/s\")\n",
    "print(f\"Densidade de potência à velocidade média: {dens_potencia:.3f} kW/m2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bfbcf0c-82d8-4add-ac69-6c70052f7ad1",
   "metadata": {},
   "source": [
    "### **5.** Análise de turbina"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a79c58e-f549-4520-bdb0-e348a25f0f9b",
   "metadata": {},
   "source": [
    "Vamos agora analisar uma turbina Senvion MM92, como as que estão presentes no parque eólico em estudo. Para isso, vamos usar a biblioteca `windpowerlib`, e considerar que o rotor se encontra a 100 m de altura."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6be3e04c-fe9f-46c5-b7cc-5efe1880f0f5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from windpowerlib import WindTurbine, ModelChain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "938532b7-a49e-4ce0-8591-04160154de71",
   "metadata": {},
   "outputs": [],
   "source": [
    "turbina = WindTurbine(turbine_type = 'MM92/2050', hub_height = 100)\n",
    "turbina"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "638aa715-68b2-47f9-b6d3-4e52e8010216",
   "metadata": {},
   "source": [
    "Repara nas características da turbina acima. Vais precisar do diâmetro do rotor mais à frente."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4dd1e3f5-3fef-4772-9541-b1826914f1a5",
   "metadata": {},
   "source": [
    "#### a) Indica a potência nominal da turbina e a sua área de varrimento."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "120d03bb-6443-457d-a7e1-003af10c5673",
   "metadata": {},
   "outputs": [],
   "source": [
    "P_nominal = turbina.nominal_power / 1e6\n",
    "print(f\"Potência nominal: {P_nominal:.2f} MW\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f104cd2c-bcf4-4e77-989d-9e10f92e4bb1",
   "metadata": {},
   "outputs": [],
   "source": [
    "area_turbina = # TODO: calcula a área de varrimento da turbina\n",
    "print(f\"Área de varrimento: {area_turbina:.2f} m2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d99377f-68d8-47ca-b2e9-45766ba68095",
   "metadata": {},
   "source": [
    "#### b) Obtém a curva característica da turbina. Apresenta-a num gráfico."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49ed42a4-32fb-4bcc-aff8-3d81640e7114",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_turbina = turbina.power_curve\n",
    "\n",
    "df_turbina.value = turbina.power_curve.value / 1e6\n",
    "df_turbina = df_turbina.rename({\n",
    "    'wind_speed': 'Velocidade do vento (m/s)',\n",
    "    'value': 'Potência da turbina (MW)'\n",
    "}, axis=1).set_index('Velocidade do vento (m/s)')\n",
    "\n",
    "df_turbina"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f81aaea-c29e-4922-b3f4-146703d617a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: completa o gráfico da curva de potência da turbina\n",
    "\n",
    "plt.xlim([df_turbina.index.min(), df_turbina.index.max()+5])\n",
    "plt.ylim(bottom=0)\n",
    "\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2838da93-ec8a-46a9-8ce2-8ba293e77e25",
   "metadata": {},
   "source": [
    "Repara que não há dados acima de uma determinada velocidade. Isto é porque para ventos muito elevados, as turbinas têm de ser desligadas por motivos de segurança."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9892c669-0c2b-4764-ae13-170ebce7ff0c",
   "metadata": {},
   "source": [
    "#### c) Compara graficamente a curva característica da turbina com a distribuição de probabilidade da velocidade do vento. O que podes concluir?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c5a507e-84b4-40ad-9bb3-86bfe2bc1a2d",
   "metadata": {},
   "source": [
    "Vamos primeiro juntar as duas dataframes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9db53fb4-ebb9-42e9-9572-d804ae7ee78e",
   "metadata": {},
   "outputs": [],
   "source": [
    "df['Potência da turbina (MW)'] = df_turbina['Potência da turbina (MW)']\n",
    "df['Potência da turbina (MW)'] = df['Potência da turbina (MW)'].fillna(0)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4fa40e20-ff2f-4b86-9def-d6f059dbfc34",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the figure and axis\n",
    "fig, ax1 = plt.subplots(figsize=(10, 6))\n",
    "\n",
    "# Plot Fração solar on the first y-axis\n",
    "ax1.plot(df['f(u)'], label=\"f(u)\", color='tab:blue')\n",
    "ax1.set_xlabel('Velocidade do vento (m/s)')\n",
    "ax1.set_ylabel('Densidade de probabilidade', color='tab:blue')\n",
    "ax1.tick_params(axis='y', labelcolor='tab:blue')\n",
    "ax1.set_ylim(bottom=0)\n",
    "\n",
    "# Create a second y-axis for Custo (€/kWh)\n",
    "ax2 = ax1.twinx()\n",
    "ax2.plot(df['Potência da turbina (MW)'], color='tab:red', label='Curva de Potência')\n",
    "ax2.set_ylabel('Potência da turbina (MW)', color='tab:red')\n",
    "ax2.tick_params(axis='y', labelcolor='tab:red')\n",
    "ax2.set_ylim(bottom=0)\n",
    "\n",
    "plt.xlim([df.index.min(), df.index.max()])\n",
    "ax1.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1ea1c705-ad02-4a00-afe4-9d3b27cf0493",
   "metadata": {},
   "source": [
    "#### d) Determina a potência média de funcionamento da turbina nesta localização."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56f57e94-fe54-440e-b995-2e4325a7ff85",
   "metadata": {},
   "outputs": [],
   "source": [
    "P_media = # TODO: determina a potência média da turbina\n",
    "print(f\"Potência média: {P_media:.2f} MW\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3825aa96-ec7c-4324-80bd-756b82eedb31",
   "metadata": {},
   "source": [
    "#### e) Determina a potência do vento em MW para a área de varrimento da turbina para cada velocidade."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "95a7f1b0-1e09-42ca-ad99-0e5218725f1f",
   "metadata": {},
   "outputs": [],
   "source": [
    "df['Potência do vento (MW)'] = # TODO: calcula a potência em MW\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62565591-44dd-42ea-97d4-bbaffd463200",
   "metadata": {},
   "source": [
    "#### f) Determina o rendimento da turbina para cada velocidade e apresenta o resultado graficamente. Compara os valores obtidos com o limite de Betz, e comenta a forma da curva obtida."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de5daf52-b1a5-4677-a74e-3d4048bf8134",
   "metadata": {},
   "outputs": [],
   "source": [
    "df['Rendimento'] = # TODO: determinar potência da turbina\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2795f83d-7282-43c3-b11f-ee51633836f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(df['Rendimento'], label='Rendimento da turbina')\n",
    "plt.axhline(16/27, color='red', label='Limite de Betz', linestyle='--')\n",
    "\n",
    "plt.xlabel('Velocidade do vento (m/s)')\n",
    "plt.ylabel('Rendimento')\n",
    "\n",
    "plt.xlim([df.index.min(), 30])\n",
    "plt.ylim(bottom=0)\n",
    "\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2d5240e-25ed-4d9c-8f24-7c1c6cc0c27f",
   "metadata": {},
   "source": [
    "#### g) Determina a energia gerada pela turbina ao longo de um ano, por velocidade. Compara graficamente com a curva de potência da turbina e comenta."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59436ff1-8800-4235-93d0-8b9d1fb6f1d1",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df['Energia gerada (MWh)'] = # TODO: calcula a energia gerada\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "04f5dd2b-82ec-4f80-820c-425f247bece5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the figure and axis\n",
    "fig, ax1 = plt.subplots(figsize=(10, 6))\n",
    "\n",
    "ax1.bar(df.index, df['Energia gerada (MWh)'], color='tab:red', alpha=0.7, label='Curva de Potência')\n",
    "ax1.set_xlabel('Velocidade do vento (m/s)')\n",
    "ax1.set_ylabel('Energia gerada (MWh)', color='tab:red')\n",
    "ax1.tick_params(axis='y', labelcolor='tab:red')\n",
    "ax1.set_ylim(bottom=0)\n",
    "\n",
    "# Plot Fração solar on the first y-axis\n",
    "ax2 = ax1.twinx()\n",
    "ax2.plot(df['Potência da turbina (MW)'], label=\"f(u)\", color='tab:blue')\n",
    "ax2.set_ylabel('Potência da turbina (MW)', color='tab:blue')\n",
    "ax2.tick_params(axis='y', labelcolor='tab:blue')\n",
    "ax2.set_ylim(bottom=0)\n",
    "\n",
    "plt.xlim([df.index.min(), 30])\n",
    "ax1.grid()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d15b5e6c-f2a5-421e-9787-a1f324feda69",
   "metadata": {},
   "source": [
    "#### h) Qual a energia total gerada ao longo de um ano?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e611b833-797d-4506-9978-824d726b2ad3",
   "metadata": {},
   "outputs": [],
   "source": [
    "E_gerada = # TODO: calcula a energia total anual\n",
    "print(f\"Energia total gerada: {round(E_gerada)} MWh\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "670b82c4-1f5e-408d-a97a-619f7208386d",
   "metadata": {},
   "source": [
    "#### i) Determina o rendimento médio da turbina."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7342e9d7-a048-40d3-9de1-150de79247ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "rendimento_medio = # TODO: calcula o rendimento\n",
    "print(f\"Rendimento médio: {rendimento_medio:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ea20c5c8-e34f-4a06-a7a5-be4be240498f",
   "metadata": {},
   "source": [
    "#### j) Determina o fator de capacidade da turbina. Qual dos dois valores, rendimento médio ou fator de capacidade, te parece mais útil para avaliar a performance da turbina num determinado local? Porquê?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84173f06-e787-4356-a715-2e791063c5e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "Fc = # TODO: calcula o fator de capacidade\n",
    "print(f\"Fator de capacidade: {Fc:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7939925-6f0b-4d07-a291-23f72f023457",
   "metadata": {},
   "source": [
    "## Comparação de cenários"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50153861-1839-44aa-9611-11276e6b70df",
   "metadata": {},
   "source": [
    "As funções definidas abaixo vão servir para compararmos o recurso e geração eólicos em diferentes cenários."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fdf70a5e-1e09-4b5d-868d-305a523ebb2f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_wind_frequency(latitude, longitude, start_date, end_date, height, bin_width=1):\n",
    "    # Step 1: Fetch wind data from Open-Meteo API\n",
    "    url = f\"https://archive-api.open-meteo.com/v1/archive?latitude={latitude}&longitude={longitude}&start_date={start_date}&end_date={end_date}&hourly=wind_speed_10m,wind_speed_100m&timezone=auto\"\n",
    "    \n",
    "    response = requests.get(url)\n",
    "    if response.status_code != 200:\n",
    "        raise Exception(f\"API request failed with status code {response.status_code}\")\n",
    "\n",
    "    # Convert JSON response to DataFrame\n",
    "    data = response.json()\n",
    "    wind_speed_10m = km_h_to_m_s(np.array(data[\"hourly\"][\"wind_speed_10m\"]))\n",
    "    wind_speed_100m = km_h_to_m_s(np.array(data[\"hourly\"][\"wind_speed_100m\"]))\n",
    "\n",
    "    # Step 2: Compute wind shear exponent (α) for each timestamp\n",
    "    h1, h2 = 10, 100\n",
    "    alpha = np.log(wind_speed_100m / wind_speed_10m) / np.log(h2 / h1)\n",
    "\n",
    "    # Step 3: Estimate wind speed at the desired height\n",
    "    wind_speed_at_height = wind_speed_10m * (height / h1) ** alpha\n",
    "    wind_speed_at_height = wind_speed_at_height[~np.isnan(wind_speed_at_height)]\n",
    "\n",
    "    # Step 4: Create a frequency distribution (histogram)\n",
    "    max_speed = np.nanmax(wind_speed_at_height)\n",
    "    bins = np.arange(0, max_speed + bin_width, bin_width)\n",
    "    hist, bin_edges = np.histogram(wind_speed_at_height, bins=bins, density=True)\n",
    "\n",
    "    # Step 5: Create DataFrame for observed frequency distribution\n",
    "    bin_centers = bins[:-1]  # Use bin edges directly instead of midpoints\n",
    "    df_observed = pd.DataFrame({\n",
    "        \"Velocidade do vento (m/s)\": bin_centers,\n",
    "        f\"f(u) (h={height}m)\": hist\n",
    "    }).set_index(\"Velocidade do vento (m/s)\")\n",
    "\n",
    "    return df_observed\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e78adf81-94d4-415a-a70b-ce5ea06effe8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_wind_turbine_energy(turbine, wind_distribution):\n",
    "    # Step 1: Get the turbine power curve\n",
    "    df = turbine.power_curve.copy()\n",
    "    df = df.rename(columns={\"wind_speed\": \"Velocidade do vento (m/s)\",\n",
    "                                            \"value\": \"Potência da turbina (MW)\"}).set_index(\"Velocidade do vento (m/s)\")\n",
    "\n",
    "    # Step 2: Ensure power output is only taken for exact wind speeds in the power curve\n",
    "    power_output = wind_distribution.index.map(lambda ws: df[\"Potência da turbina (MW)\"].get(ws, 0))\n",
    "\n",
    "    # Step 3: Compute total annual energy generated (MWh)\n",
    "    hours_per_year = 8760  # Total hours in a year\n",
    "    total_energy = np.nansum(power_output * wind_distribution.values) * hours_per_year  # MWh\n",
    "\n",
    "    # Step 4: Compute capacity factor (%)\n",
    "    nominal_power = turbine.nominal_power / 1e6  # Convert W to MW\n",
    "    capacity_factor = (total_energy / (nominal_power * hours_per_year))\n",
    "\n",
    "    return {\n",
    "        \"Energia gerada (MWh)\": total_energy,\n",
    "        \"Fator de capacidade\": capacity_factor\n",
    "    }\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f64ba138-4716-467a-a1fa-ddd875566144",
   "metadata": {},
   "source": [
    "### **6.** Considera turbinas no Parque Eólico de Marvila, com o rotor a diferentes alturas."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "195389c2-8bd0-4a50-88aa-7e7e0ccd7dc7",
   "metadata": {},
   "source": [
    "#### a) Determina a distribuição do vento para as diferentes alturas e apresenta os resultados num gráfico. Comenta."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6c1c03d-9ab2-4c06-b2e8-f69b7a06fd14",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "heights = [50, 80, 100, 120]  # List of heights to iterate over\n",
    "\n",
    "# Initialize an empty DataFrame to store all results\n",
    "df_all_heights = pd.DataFrame()\n",
    "\n",
    "# Loop over each height\n",
    "for height in heights:\n",
    "    print(f\"Processing wind speed at {height}m...\")\n",
    "    \n",
    "    # Get Weibull frequency distribution for the current height\n",
    "    df_wind = get_wind_frequency(latitude, longitude, start_date, end_date, height)\n",
    "\n",
    "    # Rename column for clarity\n",
    "    df_wind = df_wind.rename(columns={df_wind.columns[0]: f\"f(u) (h={height}m)\"})\n",
    "\n",
    "    # Merge with the main DataFrame\n",
    "    if df_all_heights.empty:\n",
    "        df_all_heights = df_wind  # First iteration, initialize DataFrame\n",
    "    else:\n",
    "        df_all_heights = df_all_heights.join(df_wind, how=\"outer\")  # Merge results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f28976fb-5739-4414-8220-8b238b78a0dd",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "df_all_heights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e88a4554-2896-4103-bd2a-b2ccb29979e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,5))\n",
    "\n",
    "for height in heights:\n",
    "    plt.plot(df_all_heights.index, df_all_heights[f\"f(u) (h={height}m)\"], label=f\"{height}m\")\n",
    "\n",
    "# Labels and title\n",
    "plt.xlabel(\"Velocidade do vento (m/s)\")\n",
    "plt.ylabel(\"Frequência\")\n",
    "plt.legend()\n",
    "plt.grid(True, linestyle=\"--\", alpha=0.7)\n",
    "\n",
    "# Show the plot\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3a9aa83-aaa7-4675-9a15-6d55ea591e6e",
   "metadata": {},
   "source": [
    "#### b) Considerando o mesmo tipo de turbinas que se encontram presentes no Parque Eólico de Marvila, calcula para as diversas alturas a energia total gerada ao longo de um ano e o fator de capacidade."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fd79d05-b3b8-40a3-acc5-09a21ad415ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an empty list to store results\n",
    "results_list = []\n",
    "heights = []  # List to store the extracted height values\n",
    "\n",
    "for col in df_all_heights.columns:\n",
    "    height = int(col.split(\"=\")[-1].replace(\"m)\", \"\"))  # Extract numeric height\n",
    "    results = compute_wind_turbine_energy(turbina, df_all_heights[col])\n",
    "    results[\"Altura (m)\"] = height  # Store the height as a number\n",
    "    results_list.append(results)  # Append results to the list\n",
    "\n",
    "# Convert results list into a DataFrame\n",
    "df_results = pd.DataFrame(results_list)\n",
    "\n",
    "# Set \"Height\" as the index and sort by height\n",
    "df_results.set_index(\"Altura (m)\", inplace=True)\n",
    "df_results.sort_index(inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e11f669-df31-405e-ad4a-3a10010c93da",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "356b41fb-ee0e-4f65-9bf3-b78b393c5b7e",
   "metadata": {},
   "source": [
    "### **7.** Repete a questão anterior, mas desta vez para a localização do Windfloat, o primeiro projeto eólico offshore flutuante da Europa Continental."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2f89c814-bc82-4bc5-8898-b9cd3b2947c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "latitude_offshore = 41.670843\n",
    "longitude_offshore = -9.283087"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c760ccd-4921-4587-807a-6b91e79c622e",
   "metadata": {},
   "outputs": [],
   "source": [
    "heights = [50, 80, 100, 120]  # List of heights to iterate over\n",
    "\n",
    "# Initialize an empty DataFrame to store all results\n",
    "df_all_heights_offshore = pd.DataFrame()\n",
    "\n",
    "# Loop over each height\n",
    "for height in heights:\n",
    "    print(f\"Processing wind speed at {height}m...\")\n",
    "    \n",
    "    # Get Weibull frequency distribution for the current height\n",
    "    df_weibull = get_wind_frequency(latitude_offshore, longitude_offshore, start_date, end_date, height)\n",
    "\n",
    "    # Rename column for clarity\n",
    "    df_weibull = df_weibull.rename(columns={df_weibull.columns[0]: f\"f(u) (h={height}m)\"})\n",
    "\n",
    "    # Merge with the main DataFrame\n",
    "    if df_all_heights_offshore.empty:\n",
    "        df_all_heights_offshore = df_weibull  # First iteration, initialize DataFrame\n",
    "    else:\n",
    "        df_all_heights_offshore = df_all_heights_offshore.join(df_weibull, how=\"outer\")  # Merge results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ecd05f4d-0de4-4ffc-8f27-52fa610dec25",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_all_heights_offshore"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e373bb19-9943-4812-8d14-c54a9cd66af7",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,5))\n",
    "\n",
    "for height in heights:\n",
    "    plt.plot(df_all_heights_offshore.index, df_all_heights_offshore[f\"f(u) (h={height}m)\"], label=f\"{height}m\")\n",
    "\n",
    "# Labels and title\n",
    "plt.xlabel(\"Velocidade do vento (m/s)\")\n",
    "plt.ylabel(\"Frequência\")\n",
    "plt.legend()\n",
    "plt.grid(True, linestyle=\"--\", alpha=0.7)\n",
    "\n",
    "# Show the plot\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "530e0621-9d73-4024-9454-2ed52e8a1d0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an empty list to store results\n",
    "results_list = []\n",
    "heights = []  # List to store the extracted height values\n",
    "\n",
    "for col in df_all_heights_offshore.columns:\n",
    "    height = int(col.split(\"=\")[-1].replace(\"m)\", \"\"))  # Extract numeric height\n",
    "    results = compute_wind_turbine_energy(turbina, df_all_heights_offshore[col])\n",
    "    results[\"Altura (m)\"] = height  # Store the height as a number\n",
    "    results_list.append(results)  # Append results to the list\n",
    "\n",
    "# Convert results list into a DataFrame\n",
    "df_results_offshore = pd.DataFrame(results_list)\n",
    "\n",
    "# Set \"Height\" as the index and sort by height\n",
    "df_results_offshore.set_index(\"Altura (m)\", inplace=True)\n",
    "df_results_offshore.sort_index(inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1803def8-9aa4-44ef-967c-92dfd0117dea",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_results_offshore"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14fd0c78-1500-4607-b2f3-c09b1694b22d",
   "metadata": {},
   "source": [
    "#### a) Apresenta um gráfico a comparar as distribuições de velocidade do vento para o caso onshore e offshore. Comenta."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ed4632f-b4c0-49d0-8bcd-a13c753f039d",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 6))\n",
    "\n",
    "# Plot each height for df_all_heights\n",
    "for column in df_all_heights.columns:\n",
    "    plt.plot(df_all_heights.index, df_all_heights[column], label=f\"Onshore - {column}\", linestyle='-', alpha=0.7)\n",
    "\n",
    "# Plot each height for df_all_heights_offshore\n",
    "for column in df_all_heights_offshore.columns:\n",
    "    plt.plot(df_all_heights_offshore.index, df_all_heights_offshore[column], label=f\"Offshore - {column}\", linestyle='--', alpha=0.7)\n",
    "\n",
    "# Labels and title\n",
    "plt.xlabel(\"Velocidade do vento (m/s)\")\n",
    "plt.ylabel(\"Frequência\")\n",
    "plt.legend()\n",
    "plt.grid(True, linestyle=\"--\", alpha=0.6)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63c0d512-a11f-436e-af2b-baedc6174574",
   "metadata": {},
   "source": [
    "#### b) Apresenta um gráfico a comparar os fatores de capacidade do caso onshore e offshore ao longo das várias alturas. Comenta."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d6212c1b-6d19-448f-b2df-672f3adb95ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(df_results[\"Fator de capacidade\"], label=\"Parque Eólico de Marvila\")\n",
    "plt.plot(df_results_offshore[\"Fator de capacidade\"], label=\"Windfloat\")\n",
    "\n",
    "plt.ylim(bottom=0)\n",
    "\n",
    "plt.xlabel('Altura do rotor (m)')\n",
    "plt.ylabel('Fator de capacidade')\n",
    "\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "53b209ed-8794-4d99-ac91-d861ac1d6b48",
   "metadata": {},
   "source": [
    "# Bónus"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce05de98-0db0-4fae-8b50-7fffd5395c24",
   "metadata": {},
   "source": [
    "Para um valor extra, constrói um gráfico diferente dos anteriores que mostre algo interessante relacionado com energia eólica, e explica o que ele demonstra :)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12a3cd85-3048-4c46-b791-dee85a5cbec4",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}