{ "cells": [ { "cell_type": "markdown", "id": "8c99517a", "metadata": {}, "source": [ "### What this script does\n", "- Loads corrected LWP parquet and reads all radiometer .dat files (SWR↓=SR15D1Dn_Irr, LWR↓=IR20Dn, LWR↑=IR20Up) for Mar–Jun.\n", "\n", "- Computes solar position/zenith with pvlib for each timestamp.\n", "\n", "- Aggregates and saves combined SWR and LWR parquet files.\n", "\n", "- Merges radiation with solar geometry and LWP_Corrected, filters low SWR, and makes scatter plots:\n", "\n", "- SWR↓ vs solar zenith (colored by LWP_Corrected; several styles).\n", "\n", "- Computes clear-sky GHI (Ineichen model), plots measured SWR vs clear-sky GHI vs zenith, and derives Clear-Sky Index (CSI).\n", "\n", "#### Edit before running\n", "1) Path to corrected LWP parquet (MWR):\n", "- base_dir_mwr = r\"C:\\path\\to\\your\\Microwave_radiometer\"\n", "- filename_lwp_corrected = \"Corrected_LWP_Data.parquet\"\n", "- parquet_path_lwp_corrected = os.path.join(base_dir_mwr, filename_lwp_corrected)\n", "\n", "2) Base folder with daily radiometer .dat files\n", "\n", " base_folder = r\"C:\\path\\to\\your\\radiometer\"\n", "\n", "3) Months you want to scan: \n", " - for month_folder in ['2024-03', '2024-04', '2024-05', '2024-06']:\n", "\n", " ...\n", "\n", "4) Site coordinates (used by pvlib for solar geometry / clear-sky GHI)\n", "- latitude = 52.6324\n", "- longitude = 4.7534\n", "- altitude = 0\n", "\n", "5) Column names in your .dat files (change if your logger uses different labels)\n", "- swr_column = 'SR15D1Dn_Irr' # Shortwave down (W/m^2)\n", "- lwr_column = 'IR20Dn' # Longwave down (W/m^2)\n", "- lwr_up_column = 'IR20Up' # Longwave up (W/m^2)\n", "\n", "6) Output parquet filenames (saved under base_folder)\n", "- parquet_file_swr = os.path.join(base_folder, 'Combined_SWR_Data.parquet')\n", "- parquet_file_lwr = os.path.join(base_folder, 'Combined_LWR_Data.parquet')\n", "\n", "7) Plot/filter settings\n" ] }, { "cell_type": "code", "execution_count": null, "id": "4a8f7bdb", "metadata": {}, "outputs": [], "source": [ "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import os\n", "import datetime\n", "\n", "!pip install pvlib\n", "\n", "from matplotlib.dates import DateFormatter\n", "from datetime import time\n", "import matplotlib.colors as mcolors\n" ] }, { "cell_type": "code", "execution_count": null, "id": "522299ed", "metadata": {}, "outputs": [], "source": [ "import pvlib\n" ] }, { "cell_type": "code", "execution_count": null, "id": "84205a2f", "metadata": {}, "outputs": [], "source": [ "#Edit before running!!\n", "# 1) Path to corrected LWP parquet (MWR)\n", "# Example: r\"D:\\Thesis\\data\\Microwave_radiometer\\Corrected_LWP_Data.parquet\"\n", "base_dir_mwr = r\"C:\\path\\to\\your\\Microwave_radiometer\"\n", "\n", "# Define the filenames for the Parquet files\n", "filename_lwp_corrected = 'Corrected_LWP_Data.parquet'\n", "\n", "\n", "# Create the full paths for loading the Parquet files\n", "parquet_path_lwp_corrected = os.path.join(base_dir_mwr, filename_lwp_corrected)\n", "\n", "# Load the Parquet files into DataFrames\n", "df_lwp_corrected = pd.read_parquet(parquet_path_lwp_corrected)\n", "\n", "# Display the first few rows of each DataFrame to verify\n", "print(\"Loaded Corrected LWP Data:\")\n", "print(df_lwp_corrected.head())\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "d4552709", "metadata": {}, "outputs": [], "source": [ "print(df_lwp_corrected)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "eb94a7fe", "metadata": { "scrolled": true }, "outputs": [], "source": [ "#Edit before running!!\n", "# Base folder path containing monthly folders\n", "## 2) Base folder with daily radiometer .dat files\n", "# Example: r\"D:\\Thesis\\data\\radiometer\"\n", "base_folder = r\"C:\\path\\to\\your\\radiometer\"\n", "\n", "# Initialize lists to store data\n", "all_data_swr_dn = []\n", "all_data_lwr_dn = []\n", "all_solar_angles = []\n", "\n", "# Location information for Alkmaar, Netherlands\n", "latitude = 52.6324 # degrees North\n", "longitude = 4.7534 # degrees East\n", "altitude = 0 # altitude above sea level in meters (assuming ground level)\n", "\n", "# Loop through each month's folder\n", "for month_folder in ['2024-03', '2024-04', '2024-05', '2024-06']:# '2024-07']:\n", " month_folder_path = os.path.join(base_folder, month_folder)\n", " \n", " # Check if the folder exists\n", " if not os.path.exists(month_folder_path):\n", " continue\n", " \n", " # Loop through each day's folder in the month\n", " for day_folder in os.listdir(month_folder_path):\n", " day_folder_path = os.path.join(month_folder_path, day_folder)\n", " \n", " # Get a list of all .dat files in the folder\n", " data_files_solar = [file for file in os.listdir(day_folder_path) if file.endswith('.dat')]\n", " \n", " # Loop through each file in the folder\n", " for file in data_files_solar:\n", " # Construct the full file path\n", " file_path = os.path.join(day_folder_path, file)\n", " \n", " # Read the data from the file\n", " data = pd.read_csv(file_path, skiprows=1, delimiter=',', encoding='latin1')\n", " \n", " # Convert TIMESTAMP column to datetime format with error handling\n", " data['TIMESTAMP'] = pd.to_datetime(data['TIMESTAMP'], errors='coerce')\n", " \n", " # Drop rows with NaT values in the TIMESTAMP column\n", " data.dropna(subset=['TIMESTAMP'], inplace=True)\n", " \n", " # Select relevant columns containing short wave down radiation and long wave down radiation data\n", " swr_column = 'SR15D1Dn_Irr' # Adjust column name if needed\n", " lwr_column = 'IR20Dn' # Adjust column name if needed\n", " lwr_up_column = 'IR20Up' # Upward long-wave radiation\n", "\n", " # Calculate solar zenith angle\n", " times = data['TIMESTAMP']\n", " solar_position = pvlib.solarposition.get_solarposition(times, latitude, longitude, altitude)\n", " \n", " # Convert the columns to numeric type\n", " data[swr_column] = pd.to_numeric(data[swr_column], errors='coerce')\n", " data[lwr_column] = pd.to_numeric(data[lwr_column], errors='coerce')\n", " data[lwr_up_column] = pd.to_numeric(data[lwr_up_column], errors='coerce')\n", "\n", " # Append the data to the lists\n", " all_data_swr_dn.append(data[['TIMESTAMP', swr_column]])\n", " all_data_lwr_dn.append(data[['TIMESTAMP', lwr_column, lwr_up_column]])\n", "\n", " all_solar_angles.append(solar_position)\n", "\n", "# Combine data from all files\n", "all_data_combined_swr = pd.concat(all_data_swr_dn, ignore_index=True)\n", "all_data_combined_lwr = pd.concat(all_data_lwr_dn, ignore_index=True)\n", "\n", "# Concatenate solar angles into a single DataFrame\n", "all_solar_angles_df = pd.concat(all_solar_angles, ignore_index=True)\n", "\n", "# Merge the dataframes on the TIMESTAMP column\n", "#final_combined_df = all_data_combined_swr.merge(all_data_combined_lwr, on='TIMESTAMP').merge(all_solar_angles_df, left_on='TIMESTAMP', right_on='apparent_zenith')\n", "\n", "# Print or process final_combined_df as needed\n", "#print(final_combined_df.head())\n" ] }, { "cell_type": "code", "execution_count": null, "id": "76d2ec8b", "metadata": {}, "outputs": [], "source": [ "# Save the combined SWR data to Parquet format\n", "parquet_file_swr = os.path.join(base_folder, 'Combined_SWR_Data.parquet')\n", "all_data_combined_swr.to_parquet(parquet_file_swr, compression='gzip')\n", "print(f\"Combined SWR DataFrame saved to {parquet_file_swr}\")\n", "\n", "# Save the combined LWR data to Parquet format\n", "parquet_file_lwr = os.path.join(base_folder, 'Combined_LWR_Data.parquet')\n", "all_data_combined_lwr.to_parquet(parquet_file_lwr, compression='gzip')\n", "print(f\"Combined LWR DataFrame saved to {parquet_file_lwr}\")" ] }, { "cell_type": "code", "execution_count": null, "id": "5e412126", "metadata": { "scrolled": true }, "outputs": [], "source": [ "#print(all_solar_angles_df)\n", "print(all_solar_angles_df)" ] }, { "cell_type": "code", "execution_count": null, "id": "15a49e94", "metadata": { "scrolled": false }, "outputs": [], "source": [ "# Filter out the zero values in shortwave down radiation\n", "filtered_data = all_data_combined_swr[all_data_combined_swr['SR15D1Dn_Irr'] > 10]\n", "filtered_angles = all_solar_angles_df.loc[filtered_data.index]\n", "\n", "# Plot the scatter plot excluding zero values in radiation\n", "plt.figure(figsize=(10, 6))\n", "\n", "# Adjust marker size (s), color (c), and transparency (alpha)\n", "plt.plot(filtered_angles['apparent_zenith'], filtered_data['SR15D1Dn_Irr'], 'ko', markersize=0.1, alpha=0.1)\n", "\n", "plt.xlabel('Apparent Solar Zenith Angle (degrees)')\n", "plt.ylabel('Shortwave Down Radiation (W/m^2)')\n", "plt.title('Shortwave Down Radiation vs Solar Zenith Angle in Alkmaar, Netherlands (March-June 2024)')\n", "# Save the plot as a PDF file\n", "#output_file = 'shortwave_radiation_vs_zenith_angle.png'\n", "#plt.savefig(output_file, format='png', dpi=300) # Save as PDF with 300 dpi for high quality\n", "\n", "# Show the plot (optional)\n", "plt.show()\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "6c7b45d6", "metadata": { "scrolled": true }, "outputs": [], "source": [ "#print(filtered_data)\n", "#print(filtered_angles)\n", "\n", "# Concatenate DataFrames along the columns\n", "merged_df = pd.concat([filtered_data, filtered_angles['apparent_zenith']], axis=1)\n", "\n", "# Print the merged DataFrame\n", "print(\"\\nMerged DataFrame:\")\n", "print(merged_df)" ] }, { "cell_type": "code", "execution_count": null, "id": "bc43aba8", "metadata": {}, "outputs": [], "source": [ "# Ensure TIMESTAMP columns are in datetime format for proper merging\n", "df_lwp_corrected['TIMESTAMP'] = pd.to_datetime(df_lwp_corrected['TIMESTAMP'])\n", "merged_df['TIMESTAMP'] = pd.to_datetime(merged_df['TIMESTAMP'])\n", "\n", "# Merge DataFrames on the 'TIMESTAMP' column\n", "merged_final_df = pd.merge(merged_df, df_lwp_corrected, on='TIMESTAMP', how='inner')\n", "\n", "# Print the final merged DataFrame\n", "print(\"\\nFinal Merged DataFrame:\")\n", "print(merged_final_df.head())" ] }, { "cell_type": "code", "execution_count": null, "id": "bd69f122", "metadata": {}, "outputs": [], "source": [ "filtered_df = merged_final_df[(merged_final_df['LWP_Corrected'] >= 0) & (merged_final_df['LWP_Corrected'] <= 50)]\n", "\n", "# Define LWP ranges and colors\n", "\n", "bin_size = 5\n", "lwp_bins = np.arange(0, 51, bin_size) \n", "#lwp_bins = np.arange(0, merged_final_df['LWP_Corrected'].max() + bin_size, bin_size) # Define ranges with step size of 200\n", "lwp_labels = [f'{b}-{b + bin_size}' for b in lwp_bins[:-1]]\n", "\n", "# Create a column for LWP range labels\n", "#merged_final_df['LWP_Range'] = pd.cut(merged_final_df['LWP_Corrected'], bins=lwp_bins, labels=lwp_labels, include_lowest=True)\n", "# Create a column for LWP range labels within the filtered DataFrame using .loc\n", "filtered_df.loc[:, 'LWP_Range'] = pd.cut(filtered_df['LWP_Corrected'], bins=lwp_bins, labels=lwp_labels, include_lowest=True)\n", "\n", "# Set up the color map\n", "cmap = plt.get_cmap('viridis') # Use a colormap from matplotlib\n", "norm = mcolors.BoundaryNorm(boundaries=lwp_bins, ncolors=cmap.N)\n", "\n", "# Plotting\n", "plt.figure(figsize=(12, 8))\n", "\n", "scatter = plt.scatter(\n", " filtered_df['apparent_zenith'],\n", " filtered_df['SR15D1Dn_Irr'],\n", " c=filtered_df['LWP_Corrected'],\n", " cmap=cmap,\n", " norm=norm,\n", " s=0.5,\n", " alpha=0.7\n", ")\n", "\n", "#scatter = plt.scatter(merged_final_df['apparent_zenith'],merged_final_df['SR15D1Dn_Irr'],c=merged_final_df['LWP_Corrected'],\n", " # cmap=cmap, norm=norm,s=0.5, # Size of the points alpha=0.7)\n", "\n", "# Add colorbar with labels\n", "cbar = plt.colorbar(scatter, boundaries=lwp_bins, ticks=lwp_bins, format='%1.1f')\n", "cbar.set_label('LWP Corrected (g/m²)')\n", "\n", "# Set colorbar ticks and labels\n", "cbar.set_ticks(lwp_bins[:-1] + bin_size / 2) # Set ticks to the middle of each bin\n", "cbar.set_ticklabels([f'{lwp_bins[i]}-{lwp_bins[i + 1]}' for i in range(len(lwp_bins) - 1)]) # Set labels\n", "\n", "# Labels and title\n", "plt.xlabel('Apparent Zenith (degrees)')\n", "plt.ylabel('Radiation (SR15D1Dn_Irr)')\n", "plt.title('Radiation vs Apparent Zenith with LWP Color Mapping')\n", "\n", "# Show grid\n", "plt.grid(True)\n", "\n", "# Show plot\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "id": "c7b9a406", "metadata": {}, "outputs": [], "source": [ "print(filtered_df)" ] }, { "cell_type": "code", "execution_count": null, "id": "f8ebb066", "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(12, 8))\n", "scatter = plt.scatter(\n", " filtered_df['apparent_zenith'],\n", " filtered_df['SR15D1Dn_Irr'],\n", " c=filtered_df['LWP_Corrected'],\n", " cmap='viridis',\n", " alpha=filtered_df['LWP_Corrected']/ filtered_df['LWP_Corrected'].max(),\n", " s=0.5\n", ")\n", "\n", "# Add color bar\n", "cbar = plt.colorbar(scatter)\n", "cbar.set_label('LWP Corrected (g/m²)')\n", "\n", "# Labels and title\n", "plt.xlabel('Apparent Zenith (degrees)')\n", "plt.ylabel('Radiation (SR15D1Dn_Irr)')\n", "plt.title('Radiation vs Apparent Zenith with Opacity-Based LWP Mapping')\n", "\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "id": "d3dc6dfd", "metadata": {}, "outputs": [], "source": [ "# Plot with continuous color mapping\n", "plt.figure(figsize=(12, 8))\n", "scatter = plt.scatter(\n", " filtered_df['apparent_zenith'],\n", " filtered_df['SR15D1Dn_Irr'],\n", " c=filtered_df['LWP_Corrected'],\n", " cmap='viridis', # Or any other color map\n", " s=0.5, # Size of the points\n", " alpha=0.7\n", ")\n", "\n", "# Add color bar\n", "cbar = plt.colorbar(scatter)\n", "cbar.set_label('LWP Corrected (g/m²)')\n", "\n", "# Labels and title\n", "plt.xlabel('Apparent Zenith (degrees)')\n", "plt.ylabel('Radiation (SR15D1Dn_Irr)')\n", "plt.title('Radiation vs Apparent Zenith with Continuous LWP Color Mapping')\n", "\n", "# Show plot\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "id": "ca619aba", "metadata": { "scrolled": false }, "outputs": [], "source": [ "'''\n", "# Plot shortwave down radiation vs solar zenith angle\n", "\n", "plt.figure(figsize=(10, 6))\n", "\n", "# Adjust marker size (s), color (c), and transparency (alpha)\n", "plt.scatter(all_solar_angles_df['apparent_zenith'], all_data_combined_swr['SR15D1Dn_Irr'], marker='o', s=0.1, c='black', alpha=0.1)\n", "\n", "plt.xlabel('Apparent Solar Zenith Angle (degrees)')\n", "plt.ylabel('Shortwave Down Radiation (W/m^2)')\n", "plt.title('Shortwave Down Radiation vs Solar Zenith Angle in Alkmaar, Netherlands (April-May 2024)')\n", "#plt.grid(True) # Enable grid for better visualization of data density\n", "plt.show()\n", "'''" ] }, { "cell_type": "code", "execution_count": null, "id": "c0337aee", "metadata": {}, "outputs": [], "source": [ "def calculate_clear_sky_sw_down_again(latitude, longitude, timestamps):\n", " # Create location object\n", " location = pvlib.location.Location(latitude, longitude)\n", " \n", " # Calculate clear-sky GHI using Ineichen model\n", " clearsky = location.get_clearsky(timestamps, model='ineichen')\n", " \n", " # Extract timestamps and clear-sky GHI values\n", " clear_sky_sw_down = clearsky['ghi']\n", " timestamps = clear_sky_sw_down.index\n", " \n", " # Create DataFrame with columns 'TIMESTAMP' and 'ghi'\n", " clear_sky_df = pd.DataFrame({\n", " 'TIMESTAMP': timestamps,\n", " 'ghi': clear_sky_sw_down.values\n", " })\n", " \n", " return clear_sky_df\n" ] }, { "cell_type": "code", "execution_count": null, "id": "d6325025", "metadata": {}, "outputs": [], "source": [ "latitude = 52.6324 # Latitude of Alkmaar, Netherlands\n", "longitude = 4.7534 # Longitude of Alkmaar, Netherlands\n", "\n", "# Assuming `all_data_combined_swr` is a DataFrame containing your data\n", "timestamps = pd.to_datetime(all_data_combined_swr['TIMESTAMP'], errors='coerce')\n", "timestamps = pd.DatetimeIndex(timestamps)\n", "\n", "clear_sky_df = calculate_clear_sky_sw_down_again(latitude, longitude, timestamps)" ] }, { "cell_type": "code", "execution_count": null, "id": "a2c4ab74", "metadata": {}, "outputs": [], "source": [ "location = pvlib.location.Location(latitude, longitude)\n", "solar_position = location.get_solarposition(timestamps)\n", "#solar_position=pvlib.solarposition.get_solarposition(timestamps, latitude, longitude)" ] }, { "cell_type": "code", "execution_count": null, "id": "0a37ec16", "metadata": { "scrolled": true }, "outputs": [], "source": [ "plt.figure(figsize=(10, 6))\n", "plt.plot(solar_position['zenith'], clear_sky_df['ghi'], 'o', markersize=1)\n", "plt.xlabel('Solar Zenith Angle (degrees)')\n", "plt.ylabel('Clear-Sky GHI (W/m^2)')\n", "plt.title('Clear-Sky GHI vs. Solar Zenith Angle')\n", "plt.grid(True)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "2ea9f651", "metadata": {}, "outputs": [], "source": [ "# Filter out zero values in measured shortwave down radiation\n", "\n", "# Calculate solar position for filtered data\n", "solar_position = pvlib.solarposition.get_solarposition(filtered_data['TIMESTAMP'], latitude, longitude)\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "3cd909aa", "metadata": {}, "outputs": [], "source": [ "print(solar_position)" ] }, { "cell_type": "code", "execution_count": null, "id": "5b6502a6", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Plot Clear-Sky GHI vs. Solar Zenith Angle\n", "plt.figure(figsize=(10, 6))\n", "\n", "# Plot Measured SW_dn with black small dots and reduced opacity\n", "#plt.plot(filtered_angles['apparent_zenith'], filtered_data['SR15D1Dn_Irr'], 'ko', markersize=0.1, alpha=0.1, label='Measured SW_dn')\n", "plt.plot(merged_df['apparent_zenith'], merged_df['SR15D1Dn_Irr'], 'ko', markersize=0.1, alpha=0.1, label='Measured SW_dn')\n", "\n", "# Plot Clear-Sky GHI with a red line on top\n", "plt.plot(merged_df['apparent_zenith'], clear_sky_df.loc[merged_df.index, 'ghi'], 'ro',markersize=0.1, alpha=0.5,label='Clear-Sky GHI', zorder=5)\n", "plt.xlabel('Solar Zenith Angle (degrees)')\n", "plt.ylabel('Shortwave Down Radiation (W/m^2)')\n", "plt.title('Shortwave Down Radiation vs Solar Zenith Angle in Alkmaar, Netherlands (March-July 2024)')\n", "#plt.legend()\n", "plt.grid(True)\n", "# Save the plot as a PDF file\n", "#output_file = 'shortwave_radiation_vs_zenith_angle_vs_csmodel.png'\n", "# Save the plot as a PNG file\n", "#plt.savefig(output_file, format='png', dpi=300) # Save as PNG with 300 dpi for high quality\n", "plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "id": "8f7db910", "metadata": {}, "outputs": [], "source": [ "\n", "# Ensure `ghi` from clear_sky_df is aligned with the timestamps in `merged_df`\n", "merged_df['ghi'] = clear_sky_df.loc[merged_df.index, 'ghi'] # Add 'ghi' values from clear_sky_df to merged_df\n", "merged_df['CSI'] = merged_df['SR15D1Dn_Irr'] / merged_df['ghi'] # Calculate CSI\n", "\n", "# Check the first few rows to confirm the new CSI column\n", "print(merged_df[['TIMESTAMP', 'SR15D1Dn_Irr', 'ghi', 'CSI']].head())" ] }, { "cell_type": "code", "execution_count": null, "id": "a939eb35", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Merge on TIMESTAMP with inner join to keep only overlapping times\n", "final_merged_df = pd.merge(filtered_df, merged_df[['TIMESTAMP', 'SR15D1Dn_Irr', 'ghi', 'CSI']], \n", " on='TIMESTAMP', how='inner')\n", "\n", "# Check the resulting DataFrame\n", "print(final_merged_df.head())" ] } ], "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.11.7" } }, "nbformat": 4, "nbformat_minor": 5 }