{ "cells": [ { "cell_type": "markdown", "id": "5b7d33c8", "metadata": {}, "source": [ "### What this notebook does\n", "\n", "This notebook (PMT-profiles-fast.ipynb) loads and processes meteorological measurement data collected from a mast, a radiometer, and a sonic instrument for a specific date.\n", "\n", "It automatically:\n", "\n", "- Reads all .dat files from your local folders for the selected date\n", "- Converts and cleans timestamped measurements\n", "\n", "- Calculates: Air temperature at multiple heights (2–10 m), Wind speed, pressure, humidity, and air density and other thermodynamic and other quantities from the instruments. \n", "\n", "- Creates time-series plots for temperature, wind speed, air density, and related variables\n", "\n", "- Merges all data into a single 10-minute averaged CSV file for further analysis\n", "\n", "You only need to edit:\n", "\n", "- date_str = 'YYYY-MM-DD' # the date to process\n", "- data_folder = f\"C:\\\\path\\\\to\\\\your\\\\data\\\\mast\\\\{month_str}\\\\{date_str}\"\n", "- data_folder_sonic = f\"C:\\\\path\\\\to\\\\your\\\\radiometer\\\\{month_str}\\\\{date_str}\"\n", "- data_dir = f\"C:\\\\path\\\\to\\\\your\\\\Sonic\\\\{month_str}\\\\{date_str}\"" ] }, { "cell_type": "code", "execution_count": null, "id": "0e3dcdce", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:52:45.616125Z", "iopub.status.busy": "2025-06-13T18:52:45.616125Z", "iopub.status.idle": "2025-06-13T18:52:48.245912Z", "shell.execute_reply": "2025-06-13T18:52:48.245912Z" }, "papermill": { "duration": 2.663584, "end_time": "2025-06-13T18:52:48.247419", "exception": false, "start_time": "2025-06-13T18:52:45.583835", "status": "completed" }, "tags": [] }, "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", "from matplotlib.dates import DateFormatter\n", "\n", "from datetime import time\n", "import pvlib\n", "import seaborn as sns\n", "from sklearn.linear_model import LinearRegression\n", "import matplotlib.dates as mdates\n", "from sklearn.metrics import mean_absolute_error, mean_squared_error" ] }, { "cell_type": "code", "execution_count": null, "id": "2aa7dd7c", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:52:48.279468Z", "iopub.status.busy": "2025-06-13T18:52:48.279468Z", "iopub.status.idle": "2025-06-13T18:52:48.286718Z", "shell.execute_reply": "2025-06-13T18:52:48.286718Z" }, "papermill": { "duration": 0.023276, "end_time": "2025-06-13T18:52:48.286718", "exception": false, "start_time": "2025-06-13T18:52:48.263442", "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "# Constants\n", "Cp = 1005 # Specific heat capacity of dry air at constant pressure (J/kg/K)\n", "g = 9.81 # Acceleration due to gravity (m/s^2)\n", "\n", "A=6.11*100 #Pa\n", "beta=0.067 #K^-1\n", "Ttrip=273.16 #K\n", "epsilon=0.622\n", "sigma=5.67e-8 #W*m^-2*K^-4\n", "Lv = 2.5e6 # Latent heat of vaporization in J/kg\n", "\n", "Rd=287.04\n", "Rv=461.5\n", "\n", "rho_atm=1.225 #kg/m^3\n", "m_co2=0.044 #kg/mole molecular mass CO2\n", "m_atm=0.028 #molecular mass atmosphere" ] }, { "cell_type": "code", "execution_count": null, "id": "0251418c", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:52:48.310960Z", "iopub.status.busy": "2025-06-13T18:52:48.310960Z", "iopub.status.idle": "2025-06-13T18:52:48.327045Z", "shell.execute_reply": "2025-06-13T18:52:48.326038Z" }, "papermill": { "duration": 0.031841, "end_time": "2025-06-13T18:52:48.327045", "exception": false, "start_time": "2025-06-13T18:52:48.295204", "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "# Function to calculate vapor pressure (e) from relative humidity (RH) and temperature (T)\n", "def calculate_vapor_pressure(RH, T):\n", " #es = A*np.exp(beta*(T+273.15-Ttrip)) # Saturation vapor pressure in Pa\n", " es=610.78*np.exp(17.2694*(T+273.15-Ttrip)/(T+273.15-35.86))\n", " e = RH * es / 100 # Vapor pressure in Pa\n", " return e\n", "\n", "# Function to calculate specific humidity (qv) from vapor pressure (e) and atmospheric pressure (p)\n", "def calculate_specific_humidity(e, p):\n", " #qv = epsilon * e *1000 / (p*100) # Specific humidity in g/kg\n", " qv=e*1000/(p*100+(((Rv/Rd)-1)*(p*100-e)))\n", " return qv\n", "\n", "#def calculate_saturation_specific_humidity(e, p):\n", " # qs = epsilon * (e*100/RH) *1000 / (p*100) # Specific humidity in g/kg\n", " # return qs\n", "\n", "def calculate_saturation_specific_humidity(T, p):\n", " es=610.78*np.exp(17.2694*(T+273.15-Ttrip)/(T+273.15-35.86))\n", " qs=es*1000/(((Rv/Rd)*(p*100-es))+es)\n", " return qs\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "id": "7d45267c", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:52:48.358678Z", "iopub.status.busy": "2025-06-13T18:52:48.358678Z", "iopub.status.idle": "2025-06-13T18:53:11.521674Z", "shell.execute_reply": "2025-06-13T18:53:11.521674Z" }, "papermill": { "duration": 23.178731, "end_time": "2025-06-13T18:53:11.521674", "exception": false, "start_time": "2025-06-13T18:52:48.342943", "status": "completed" }, "scrolled": true, "tags": [] }, "outputs": [], "source": [ "\n", "# ⚠️ Edit this path to match your local setup before running.\n", "\n", "date_str = '2024-05-23' # ← only line you normally edit\n", "month_str = date_str[:7] # '2024-05'\n", "\n", "# Example placeholder path (safe to show publicly)\n", "data_folder = f\"C:\\\\path\\\\to\\\\your\\\\data\\\\mast\\\\{month_str}\\\\{date_str}\"\n", "\n", "# Get a list of all .dat files in the folder\n", "data_files = [file for file in os.listdir(data_folder) if file.endswith('.dat')]\n", "\n", "# Heights corresponding to each temperature measurement (in meters)\n", "heights = [2,2.99,4.47,6.69,10]\n", "\n", "\n", "# Initialize an empty DataFrame to store the combined data\n", "all_data = []\n", "\n", "# Initialize an empty list to store the dry static energy data for each height\n", "all_dse_data = []\n", "all_virtual_dse_data = [] # Initialize a list for virtual dry static energy\n", "\n", "\n", "# Initialize an empty list to store the wind speed data for each height\n", "all_ws_data = []\n", "\n", "# Initialize an empty list to store the pressure data for height 2.99\n", "all_pressure_data = [] #mbar=hPa\n", "\n", "# Initialize an empty list to store the specific humidity data for each height\n", "all_qv_data = []\n", "\n", "all_qs_data=[]\n", "\n", "# Initialize an empty list to store the relative humidity data for each height\n", "all_rh_data = []\n", "\n", "all_winddr_data=[]\n", "\n", "all_windsdv_data=[]\n", "\n", "na_counts = 0 # Initialize a counter for NaT values\n", "\n", "# Loop through each file in the folder\n", "for file in data_files:\n", " # Construct the full file path\n", " file_path = os.path.join(data_folder, file)\n", " \n", " # Read the data from the file\n", " data = pd.read_csv(file_path, skiprows=1)\n", " \n", " # Convert TIMESTAMP column to datetime format with error handling\n", " data['TIMESTAMP'] = pd.to_datetime(data['TIMESTAMP'], errors='coerce')\n", " \n", " # Count NaT values\n", " na_counts += data['TIMESTAMP'].isna().sum()\n", " # Drop rows with NaT values in the TIMESTAMP column\n", " data.dropna(subset=['TIMESTAMP'], inplace=True)\n", " \n", " # Select relevant columns containing temperature data\n", " temperature_columns = ['AirTC_E5567_Avg', 'AirTC_E5568_Avg', 'AirTC_E5569_Avg', 'AirTC_E5570_Avg','AirTC_E5571_Avg'] # Adjust column names if needed\n", " \n", " # Loop through each column and convert it to numeric type\n", " for column in temperature_columns:\n", " data[column] = pd.to_numeric(data[column], errors='coerce')\n", " # Append the data to the list\n", " all_data.append(data)\n", " # Print the number of NaT values\n", " print(f\"Number of NaT values in TIMESTAMP column: {na_counts}\")\n", " \n", " \n", " \n", " # Calculate dry static energy for each height\n", " for column, height in zip(temperature_columns, heights):\n", " # Calculate temperature in Kelvin\n", " data[f'Temperature_K_{height}'] = data[column] + 273.15\n", " \n", " # Calculate dry static energy\n", " data[f'Dry_Static_Energy_{height}'] = data[f'Temperature_K_{height}'] + g * height/Cp\n", " \n", " # Append the dry static energy data to the list\n", " all_dse_data.append(data)\n", " \n", " \n", " \n", " # Plot time series of temperature measurements for each file\n", " plt.figure(figsize=(10, 6))\n", " for column in temperature_columns:\n", " plt.plot(data['TIMESTAMP'], data[column], label=column)\n", " \n", " # Format the x-axis to show time in HH:MM format\n", " date_format = DateFormatter('%H:%M')\n", " plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", " plt.xlabel('Time (UTC)')\n", " plt.ylabel('Temperature (°C)')\n", " plt.title(f'Time Series of Temperature Measurements - {file}')\n", " plt.legend()\n", " plt.grid(True)\n", " plt.xticks(rotation=45)\n", " \n", " plt.tight_layout()\n", " plt.savefig(os.path.join(data_folder, f'{file}_temperature_plot.jpg'))\n", " plt.close()\n", " \n", " # Select relevant columns containing wind speed data\n", " ws_columns = ['WS_ms_D15008_Avg', 'WS_ms_D15014_Avg', 'WS_ms_D15463_Avg']\n", " \n", " # Loop through each column and convert it to numeric type\n", " for column1 in ws_columns:\n", " data[column1] = pd.to_numeric(data[column1], errors='coerce')\n", " # Append the data to the list\n", " all_ws_data.append(data)\n", " \n", " \n", " # Plot time series of wind speed measurements for each file\n", " plt.figure(figsize=(10, 6))\n", " for column1 in ws_columns:\n", " plt.plot(data['TIMESTAMP'], data[column1], label=column1)\n", " \n", " # Format the x-axis to show time in HH:MM format\n", " date_format = DateFormatter('%H:%M')\n", " plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", " plt.xlabel('Time (UTC)')\n", " plt.ylabel('Wind speed (m/s)')\n", " plt.title(f'Time Series of Wind Speed Measurements - {file}')\n", " plt.legend()\n", " plt.grid(True)\n", " plt.xticks(rotation=45)\n", " \n", " plt.tight_layout()\n", " plt.savefig(os.path.join(data_folder, f'{file}_windspeed_plot.jpg'))\n", " plt.close()\n", " \n", " \n", " # Select relevant columns containing pressyre data\n", " p_column = 'BP_mbar_Avg'\n", " \n", " # Loop through each column and convert it to numeric type\n", " data[p_column] = pd.to_numeric(data[p_column], errors='coerce')\n", " # Append the data to the list\n", " all_pressure_data.append(data)\n", " \n", " #select RH columns\n", " rh_columns = ['RH_E5567_Avg', 'RH_E5568_Avg', 'RH_E5569_Avg', 'RH_E5570_Avg', 'RH_E5571_Avg']\n", " \n", " # Loop through each column and convert it to numeric type\n", " for column in rh_columns:\n", " data[column] = pd.to_numeric(data[column], errors='coerce')\n", " \n", " \n", " # Create empty columns for specific humidity corresponding to each height\n", " for height in heights:\n", " data[f'qv_{height}m'] = np.nan\n", " \n", " # Calculate and plot specific humidity for each height\n", " for temp_column, rh_column, height in zip(temperature_columns, rh_columns, heights):\n", " # Calculate vapor pressure (e) using the temperature for each hour\n", " T = data[temp_column] # Temperature in degrees Celsius\n", " RH = data[rh_column] # Relative humidity as percentage\n", " \n", " e = calculate_vapor_pressure(RH, T)\n", " \n", " # Atmospheric pressure (p) at height 2.99m\n", " p = data[p_column]\n", " \n", " # Calculate specific humidity (qv) in g/kg\n", " qv = calculate_specific_humidity(e, p)\n", " \n", " qs= calculate_saturation_specific_humidity(T,p)\n", " \n", " \n", " # Assign specific humidity data to the corresponding column\n", " data[f'qv_{height}m'] = qv\n", " \n", " data[f'qs_{height}m'] = qs\n", "\n", " \n", " # Calculate virtual dry static energy\n", " data[f'Virtual_Dry_Static_Energy_{height}'] = data[f'Dry_Static_Energy_{height}'] + data[f'Temperature_K_{height}']*(((Rv / Rd) - 1) * (data[f'qv_{height}m']/1000)) # Convert g/kg to kg/kg by dividing by 1000\n", " \n", " # Append the virtual DSE data for the current height\n", " all_virtual_dse_data.append(data[['TIMESTAMP'] + [f'Virtual_Dry_Static_Energy_{height}' for height in heights]])\n", " # Append the specific humidity data to the list\n", " all_qv_data.append(data[['TIMESTAMP'] + [f'qv_{height}m' for height in heights]])\n", " \n", " all_qs_data.append(data[['TIMESTAMP'] + [f'qs_{height}m' for height in heights]])\n", "\n", " # Plot specific humidity for each height\n", " plt.figure(figsize=(10, 6))\n", " for height in heights:\n", " plt.plot(data['TIMESTAMP'], data[f'qv_{height}m'], label=f'Specific Humidity - Height: {height}m')\n", " plt.xlabel('Time (UTC)')\n", " plt.ylabel('Specific Humidity (g/kg)')\n", " plt.title(f'Specific Humidity (qv) - {file}')\n", " plt.legend()\n", " plt.grid(True)\n", " plt.xticks(rotation=45)\n", " plt.tight_layout()\n", " \n", " # Save the plot\n", " plot_name = f'{file.split(\".\")[0]}_specific_humidity_plot.jpg'\n", " plt.savefig(os.path.join(data_folder, plot_name))\n", " plt.close()\n", " \n", " \n", " \n", " wind_dir_columns= ['WindDir_D15008_Avg','WindDir_D15014_Avg','WindDir_D15463_Avg']\n", " \n", " wind_dirsdv_columns =['WindDir_D15008_StDev','WindDir_D15014_StDev', 'WindDir_D15463_StDev']\n", " \n", " for column in wind_dir_columns:\n", " data[column] = pd.to_numeric(data[column], errors='coerce')\n", " # Append the data to the list\n", " all_winddr_data.append(data)\n", " \n", " for column in wind_dirsdv_columns:\n", " data[column] = pd.to_numeric(data[column], errors='coerce')\n", " # Append the data to the list\n", " \n", " all_windsdv_data.append(data)\n", " \n", " # Loop through each column and convert it to numeric type\n", " for column in rh_columns:\n", " data[column] = pd.to_numeric(data[column], errors='coerce')\n", " all_rh_data.append(data)\n", " \n", " \n", "# Concatenate all data into a single DataFrame\n", "all_dse_data = pd.concat(all_dse_data)\n", "\n", "all_ws_data= pd.concat(all_ws_data)\n", "\n", "all_pressure_data= pd.concat(all_pressure_data)\n", "\n", "all_winddr_data=pd.concat(all_winddr_data)\n", "\n", "all_windsdv_data=pd.concat(all_windsdv_data)\n", "\n", "all_qv_data = pd.concat(all_qv_data, ignore_index=True)\n", "\n", "# Concatenate all virtual DSE data into a single DataFrame\n", "all_virtual_dse_data_df = pd.concat(all_virtual_dse_data, ignore_index=True)\n", "\n", "all_rh_data= pd.concat(all_rh_data,ignore_index=True)\n", "all_qs_data= pd.concat(all_qs_data,ignore_index=True)\n", "\n", "combined_data = pd.concat(all_data)\n", "print(all_qv_data)" ] }, { "cell_type": "code", "execution_count": null, "id": "a8b7d0a7", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:50:44.499992Z", "iopub.status.busy": "2025-06-13T18:50:44.499992Z", "iopub.status.idle": "2025-06-13T18:50:44.526776Z", "shell.execute_reply": "2025-06-13T18:50:44.526776Z" }, "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "print(all_virtual_dse_data_df)" ] }, { "cell_type": "code", "execution_count": null, "id": "c4882a0e", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:50:44.563677Z", "iopub.status.busy": "2025-06-13T18:50:44.563677Z", "iopub.status.idle": "2025-06-13T18:50:51.164702Z", "shell.execute_reply": "2025-06-13T18:50:51.164702Z" }, "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "\n", "# Plot time series of temperature measurements for the entire day\n", "plt.figure(figsize=(10, 6))\n", "for column in temperature_columns:\n", " plt.plot(combined_data['TIMESTAMP'], combined_data[column], label=column)\n", " \n", " \n", "# Format the x-axis to show time in HH:MM format\n", "date_format = DateFormatter('%H:%M')\n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", "plt.xlabel('Time (UTC)')\n", "plt.ylabel('Temperature (°C)')\n", "plt.title('Time Series of Temperature Measurements for the Whole Day')\n", "plt.legend()\n", "plt.grid(True)\n", "plt.xticks(rotation=45)\n", "plt.tight_layout()\n", "\n", "# Save the plot as a JPEG file in the same folder\n", "#output_file = os.path.join(data_folder, 'temperature_time_series.jpg')\n", "#plt.savefig(output_file, format='jpg')\n", "# Save the plot as a PDF file in the same folder\n", "output_file = os.path.join(data_folder, 'temperature_time_series.pdf')\n", "plt.savefig(output_file, format='pdf', dpi=300) # Save as PDF for better quality\n", "# Show the plot (optional)\n", "plt.show()\n", "\n", "# Plot the dry static energy over the whole day for each height\n", "plt.figure(figsize=(10, 6))\n", "for height in heights:\n", " plt.plot(all_dse_data['TIMESTAMP'], all_dse_data[f'Dry_Static_Energy_{height}'], label=f'{height}m')\n", "\n", " \n", "# Format the x-axis to show time in HH:MM format\n", "date_format = DateFormatter('%H:%M')\n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", "plt.xlabel('Time (UTC)')\n", "plt.ylabel('Dry Static Energy (over heat capacity) (K)')\n", "plt.title('Dry Static Energy Variation Over the Whole Day')\n", "plt.legend(title='Height', loc='upper right')\n", "plt.grid(True)\n", "plt.xticks(rotation=45)\n", "plt.tight_layout()\n", "\n", "# Save the plot as a JPEG file in the same folder\n", "#output_file = os.path.join(data_folder, 'drystaticenergy_time_series.jpg')\n", "#plt.savefig(output_file, format='jpg')\n", "\n", "# Save the plot as a PDF file in the same folder\n", "output_file = os.path.join(data_folder, 'drystaticenergy_time_series.pdf')\n", "plt.savefig(output_file, format='pdf', dpi=300) # Save as PDF for better quality\n", "# Show the plot (optional)\n", "plt.show()\n", "\n", "# Plot the virtual dry static energy over the whole day for each height\n", "plt.figure(figsize=(10, 6))\n", "for height in heights:\n", " plt.plot(all_virtual_dse_data_df['TIMESTAMP'], all_virtual_dse_data_df[f'Virtual_Dry_Static_Energy_{height}'], label=f'{height}m')\n", "\n", "# Format the x-axis to show time in HH:MM format\n", "date_format = DateFormatter('%H:%M')\n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", "plt.xlabel('Time (UTC)')\n", "plt.ylabel('Virtual Dry Static Energy (over heat capacity) (K)')\n", "plt.title('Virtual Dry Static Energy Variation Over the Whole Day')\n", "plt.legend(title='Height', loc='upper right')\n", "plt.grid(True)\n", "plt.xticks(rotation=45)\n", "plt.tight_layout()\n", "\n", "# Save the plot as a PNG file in the same folder\n", "output_file = os.path.join(data_folder, 'virtual_dry_static_energy_time_series.png')\n", "plt.savefig(output_file, format='png', dpi=300) # Save as PNG for better quality\n", "\n", "# Show the plot (optional)\n", "plt.show()\n", "\n", "\n", "\n", "# Plot time series of temperature measurements for the entire day\n", "plt.figure(figsize=(10, 6))\n", "for column in ws_columns:\n", " plt.plot(all_ws_data['TIMESTAMP'], all_ws_data[column], label=column)\n", " \n", " \n", " \n", "# Format the x-axis to show time in HH:MM format\n", "date_format = DateFormatter('%H:%M')\n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", "#plt.ylim(0, 15) \n", "\n", "plt.xlabel('Time')\n", "plt.ylabel('Wind Speed (m/s)')\n", "plt.title('Time Series of Wind Speed Measurements for the Whole Day')\n", "plt.legend()\n", "plt.grid(True)\n", "plt.xticks(rotation=45)\n", "plt.tight_layout()\n", "\n", "# Save the plot as a JPEG file in the same folder\n", "#output_file = os.path.join(data_folder, 'windspeed_time_series.jpg')\n", "#plt.savefig(output_file, format='jpg')\n", "\n", "output_file = os.path.join(data_folder, 'windspeed_time_series.pdf')\n", "plt.savefig(output_file, format='pdf', dpi=300) # Save as PDF for better quality\n", "\n", "# Show the plot (optional)\n", "plt.show()\n", "\n", "# Plot time series of temperature measurements for the entire day\n", "plt.figure(figsize=(10, 6))\n", "\n", "plt.plot(all_pressure_data['TIMESTAMP'], all_pressure_data[p_column], label=p_column)\n", " \n", " \n", "# Format the x-axis to show time in HH:MM format\n", "date_format = DateFormatter('%H:%M')\n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", "plt.xlabel('Time')\n", "plt.ylabel('Pressure at 2.99 m(mbar)')\n", "plt.title('Time Series of Pressure Measurements for the Whole Day')\n", "plt.legend()\n", "plt.grid(True)\n", "plt.xticks(rotation=45)\n", "plt.tight_layout()\n", "\n", "# Save the plot as a JPEG file in the same folder\n", "#output_file = os.path.join(data_folder, 'pressure_time_series.jpg')\n", "#plt.savefig(output_file, format='jpg')\n", "\n", "output_file = os.path.join(data_folder, 'pressure_time_series.pdf')\n", "plt.savefig(output_file, format='pdf', dpi=300) # Save as PDF for better quality\n", "# Show the plot (optional)\n", "plt.show()\n", "\n", "#print(all_pressure_data[p_column].mean())\n" ] }, { "cell_type": "code", "execution_count": null, "id": "f9cbf4c1", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:50:51.247998Z", "iopub.status.busy": "2025-06-13T18:50:51.247998Z", "iopub.status.idle": "2025-06-13T18:50:56.022560Z", "shell.execute_reply": "2025-06-13T18:50:56.022560Z" }, "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "\n", "# 1. Plot time series of temperature measurements for the entire day\n", "plt.figure(figsize=(10, 6))\n", "for column in temperature_columns:\n", " plt.plot(combined_data['TIMESTAMP'], combined_data[column], label=column)\n", "\n", "# Format the x-axis to show time in HH:MM format\n", "date_format = DateFormatter('%H:%M')\n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", "# Set labels, title, and legend with larger font sizes\n", "plt.xlabel('Time (UTC)', fontsize=16)\n", "plt.ylabel('Temperature (°C)', fontsize=16)\n", "plt.title('Time Series of Temperature Measurements for the Whole Day', fontsize=18)\n", "plt.legend(fontsize=14)\n", "plt.grid(True)\n", "plt.xticks(rotation=45, fontsize=12)\n", "plt.yticks(fontsize=12)\n", "\n", "plt.tight_layout()\n", "\n", "# Save the plot as a JPEG file\n", "output_file = os.path.join(data_folder, 'temperature_time_series.jpg')\n", "plt.savefig(output_file, format='jpg', dpi=300)\n", "plt.show()\n", "\n", "\n", "# 2. Plot the dry static energy over the whole day for each height\n", "plt.figure(figsize=(10, 6))\n", "for height in heights:\n", " plt.plot(all_dse_data['TIMESTAMP'], all_dse_data[f'Dry_Static_Energy_{height}'], label=f'{height}m')\n", "\n", "# Format the x-axis to show time in HH:MM format\n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", "# Set labels, title, and legend with larger font sizes\n", "plt.xlabel('Time (UTC)', fontsize=16)\n", "plt.ylabel('Dry Static Energy (over heat capacity) (K)', fontsize=16)\n", "plt.title('Dry Static Energy Variation Over the Whole Day', fontsize=18)\n", "plt.legend(title='Height', loc='upper right', fontsize=14, title_fontsize=14)\n", "plt.grid(True)\n", "plt.xticks(rotation=45, fontsize=12)\n", "plt.yticks(fontsize=12)\n", "\n", "plt.tight_layout()\n", "\n", "# Save the plot as a JPEG file\n", "output_file = os.path.join(data_folder, 'drystaticenergy_time_series.jpg')\n", "plt.savefig(output_file, format='jpg', dpi=300)\n", "plt.show()\n", "\n", "\n", "# 3. Plot time series of wind speed measurements for the entire day\n", "plt.figure(figsize=(10, 6))\n", "for column in ws_columns:\n", " plt.plot(all_ws_data['TIMESTAMP'], all_ws_data[column], label=column)\n", "\n", "# Format the x-axis to show time in HH:MM format\n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", "# Set labels, title, and legend with larger font sizes\n", "plt.xlabel('Time', fontsize=16)\n", "plt.ylabel('Wind Speed (m/s)', fontsize=16)\n", "plt.title('Time Series of Wind Speed Measurements for the Whole Day', fontsize=18)\n", "plt.legend(fontsize=14)\n", "plt.grid(True)\n", "plt.xticks(rotation=45, fontsize=12)\n", "plt.yticks(fontsize=12)\n", "\n", "plt.tight_layout()\n", "\n", "# Save the plot as a JPEG file\n", "output_file = os.path.join(data_folder, 'windspeed_time_series.jpg')\n", "plt.savefig(output_file, format='jpg', dpi=300)\n", "plt.show()\n", "\n", "\n", "# 4. Plot time series of pressure measurements for the entire day\n", "plt.figure(figsize=(10, 6))\n", "plt.plot(all_pressure_data['TIMESTAMP'], all_pressure_data[p_column], label=p_column)\n", "\n", "# Format the x-axis to show time in HH:MM format\n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", "# Set labels, title, and legend with larger font sizes\n", "plt.xlabel('Time', fontsize=16)\n", "plt.ylabel('Pressure at 2.99 m (mbar)', fontsize=16)\n", "plt.title('Time Series of Pressure Measurements for the Whole Day', fontsize=18)\n", "plt.legend(fontsize=14)\n", "plt.grid(True)\n", "plt.xticks(rotation=45, fontsize=12)\n", "plt.yticks(fontsize=12)\n", "\n", "plt.tight_layout()\n", "\n", "# Save the plot as a JPEG file\n", "output_file = os.path.join(data_folder, 'pressure_time_series.jpg')\n", "plt.savefig(output_file, format='jpg', dpi=300)\n", "plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "id": "f6a56a89", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:50:56.141971Z", "iopub.status.busy": "2025-06-13T18:50:56.141971Z", "iopub.status.idle": "2025-06-13T18:50:58.733512Z", "shell.execute_reply": "2025-06-13T18:50:58.733512Z" }, "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "\n", "all_pressure_data['rho_air']=all_pressure_data[p_column]*100/(Rd*(all_pressure_data['AirTC_E5568_Avg']+273.15))\n", "all_pressure_data['rho_air_Tv']=all_pressure_data[p_column]*100/(Rd*(all_pressure_data['AirTC_E5568_Avg']+273.15)*(1+(((Rv/Rd)-1)*(all_pressure_data['qv_2.99m']/1000))))\n", "#print(all_pressure_data)\n", "# Plot time series of temperature measurements for the entire day\n", "plt.figure(figsize=(10, 6))\n", "\n", "plt.plot(all_pressure_data['TIMESTAMP'], all_pressure_data['rho_air'],label='rho_air from T')\n", "plt.plot(all_pressure_data['TIMESTAMP'], all_pressure_data['rho_air_Tv'],label='rho_air from T_v')\n", " \n", " \n", "# Format the x-axis to show time in HH:MM format\n", "date_format = DateFormatter('%H:%M')\n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", "plt.xlabel('Time')\n", "plt.ylabel('rho_air (kg/m^3)')\n", "plt.title('Time Series of Density of Air at 2.99 m')\n", "#plt.legend()\n", "plt.grid(True)\n", "plt.xticks(rotation=45)\n", "plt.tight_layout()\n", "plt.legend()\n", "# Save the plot as a JPEG file in the same folder\n", "output_file = os.path.join(data_folder, 'density_air_time_series.pdf')\n", "#plt.savefig(output_file, format='jpg')\n", "plt.savefig(output_file, format='pdf', dpi=300) # Save as PDF for better quality\n", "\n", "# Show the plot (optional)\n", "plt.show()\n", "\n", "# Plot specific humidity for each height\n", "plt.figure(figsize=(10, 6))\n", "for height in heights:\n", " plt.plot(all_qv_data['TIMESTAMP'], all_qv_data[f'qv_{height}m'], label=f'{height}m')\n", " \n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", "\n", "plt.xlabel('Time (UTC)')\n", "plt.ylabel('Specific Humidity, q_v (g/kg)')\n", "plt.title('Specific Humidity Variation Over the Whole Day')\n", "plt.legend(title='Height', loc='upper right')\n", "\n", "#plt.legend()\n", "plt.grid(True)\n", "plt.xticks(rotation=45)\n", "plt.tight_layout()\n", "\n", "# Save the plot\n", "plot_name = 'specific_humidity_time_series.pdf'\n", "plt.savefig(os.path.join(data_folder, plot_name),dpi=300)\n", "#plt.savefig(output_file, format='pdf', dpi=300) # Save as PDF for better quality\n", "\n", "plt.show()\n", "\n", "# Plot relative humidity for each height\n", "plt.figure(figsize=(10, 6))\n", "for column in rh_columns:\n", " plt.plot(all_rh_data['TIMESTAMP'], all_rh_data[column], label=column)#f'Relative Humidity - Height: {height}m')\n", " \n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", "\n", "plt.xlabel('Time (UTC)')\n", "plt.ylabel('Relative Humidity (%)')\n", "plt.title('Relative Humidity Variation Over the Whole Day')\n", "plt.legend()\n", "plt.grid(True)\n", "plt.xticks(rotation=45)\n", "plt.tight_layout()\n", "\n", "# Save the plot\n", "plot_name = 'relative_humidity_time_series.pdf'\n", "plt.savefig(os.path.join(data_folder, plot_name),dpi=300)\n", "\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "376b2b88", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:50:58.879951Z", "iopub.status.busy": "2025-06-13T18:50:58.879951Z", "iopub.status.idle": "2025-06-13T18:51:02.437403Z", "shell.execute_reply": "2025-06-13T18:51:02.437403Z" }, "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "\n", "# Define date format for x-axis\n", "date_format = DateFormatter('%H:%M')\n", "\n", "# 1. Plot time series of air density measurements\n", "all_pressure_data['rho_air'] = all_pressure_data[p_column]*100/(Rd*(all_pressure_data['AirTC_E5568_Avg']+273.15))\n", "all_pressure_data['rho_air_Tv'] = all_pressure_data[p_column]*100/(Rd*(all_pressure_data['AirTC_E5568_Avg']+273.15)*(1+(((Rv/Rd)-1)*(all_pressure_data['qv_2.99m']/1000))))\n", "\n", "plt.figure(figsize=(10, 6))\n", "plt.plot(all_pressure_data['TIMESTAMP'], all_pressure_data['rho_air'], label='rho_air from T')\n", "plt.plot(all_pressure_data['TIMESTAMP'], all_pressure_data['rho_air_Tv'], label='rho_air from T_v')\n", "\n", "# Format the x-axis to show time in HH:MM format\n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", "# Set labels, title, and legend with larger font sizes\n", "plt.xlabel('Time (UTC)', fontsize=16)\n", "plt.ylabel('rho_air (kg/m³)', fontsize=16)\n", "plt.title('Time Series of Density of Air at 2.99 m', fontsize=18)\n", "plt.legend(fontsize=14)\n", "plt.grid(True)\n", "plt.xticks(rotation=45, fontsize=12)\n", "plt.yticks(fontsize=12)\n", "\n", "plt.tight_layout()\n", "\n", "# Save the plot as a JPEG file\n", "output_file = os.path.join(data_folder, 'density_air_time_series.jpg')\n", "plt.savefig(output_file, format='jpg', dpi=300)\n", "plt.show()\n", "\n", "\n", "# 2. Plot specific humidity for each height\n", "plt.figure(figsize=(10, 6))\n", "for height in heights:\n", " plt.plot(all_qv_data['TIMESTAMP'], all_qv_data[f'qv_{height}m'], label=f'{height}m')\n", "\n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", "# Set labels, title, and legend with larger font sizes\n", "plt.xlabel('Time (UTC)', fontsize=16)\n", "plt.ylabel('Specific Humidity, q_v (g/kg)', fontsize=16)\n", "plt.title('Specific Humidity Variation Over the Whole Day', fontsize=18)\n", "plt.legend(title='Height', loc='upper right', fontsize=14, title_fontsize=14)\n", "plt.grid(True)\n", "plt.xticks(rotation=45, fontsize=12)\n", "plt.yticks(fontsize=12)\n", "\n", "plt.tight_layout()\n", "\n", "# Save the plot as a JPEG file\n", "plot_name = 'specific_humidity_time_series.jpg'\n", "plt.savefig(os.path.join(data_folder, plot_name), format='jpg', dpi=300)\n", "plt.show()\n", "\n", "\n", "# 3. Plot relative humidity for each height\n", "plt.figure(figsize=(10, 6))\n", "for column in rh_columns:\n", " plt.plot(all_rh_data['TIMESTAMP'], all_rh_data[column], label=column)\n", "\n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", "# Set labels, title, and legend with larger font sizes\n", "plt.xlabel('Time (UTC)', fontsize=16)\n", "plt.ylabel('Relative Humidity (%)', fontsize=16)\n", "plt.title('Relative Humidity Variation Over the Whole Day', fontsize=18)\n", "plt.legend(fontsize=14)\n", "plt.grid(True)\n", "plt.xticks(rotation=45, fontsize=12)\n", "plt.yticks(fontsize=12)\n", "\n", "plt.tight_layout()\n", "\n", "# Save the plot as a JPEG file\n", "plot_name = 'relative_humidity_time_series.jpg'\n", "plt.savefig(os.path.join(data_folder, plot_name), format='jpg', dpi=300)\n", "plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "id": "a1f053c9", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:51:02.593525Z", "iopub.status.busy": "2025-06-13T18:51:02.592527Z", "iopub.status.idle": "2025-06-13T18:51:02.647518Z", "shell.execute_reply": "2025-06-13T18:51:02.645498Z" }, "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "\n", "# Convert TIMESTAMP to datetime if it's not already\n", "all_pressure_data['TIMESTAMP'] = pd.to_datetime(all_pressure_data['TIMESTAMP'], errors='coerce')\n", "\n", "# Set the TIMESTAMP column as the index\n", "all_pressure_data.set_index('TIMESTAMP', inplace=True)\n", "\n", "# Select only numeric columns for resampling\n", "numeric_columns = all_pressure_data.select_dtypes(include='number').columns\n", "all_pressure_data_numeric = all_pressure_data[numeric_columns]\n", "\n", "# Resample the data to 10-minute intervals and calculate the mean for each interval\n", "all_pressure_data_10min = all_pressure_data_numeric.resample('10T').mean()\n", "\n", "# Reset the index to make TIMESTAMP a column again\n", "all_pressure_data_10min.reset_index(inplace=True)\n", "\n", "# Display the resulting DataFrame\n", "print(all_pressure_data_10min)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "fbd81267", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:51:02.812661Z", "iopub.status.busy": "2025-06-13T18:51:02.812661Z", "iopub.status.idle": "2025-06-13T18:51:02.818385Z", "shell.execute_reply": "2025-06-13T18:51:02.818385Z" }, "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "print(all_pressure_data_10min.columns)" ] }, { "cell_type": "markdown", "id": "a6d9ed0c", "metadata": { "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "tags": [] }, "source": [ "### Solar/IR" ] }, { "cell_type": "code", "execution_count": null, "id": "661c2e1d", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:51:03.139654Z", "iopub.status.busy": "2025-06-13T18:51:03.139654Z", "iopub.status.idle": "2025-06-13T18:51:22.679631Z", "shell.execute_reply": "2025-06-13T18:51:22.675517Z" }, "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "# Directory containing the sonic data files\n", "# ⚠️ Edit this path to match your local setup before running.\n", "data_folder_sonic = f\"C:\\\\path\\\\to\\\\your\\\\radiometer\\\\{month_str}\\\\{date_str}\"\n", "\n", "\n", "# Get a list of all .dat files in the folder\n", "data_files_sonic = [file for file in os.listdir(data_folder_sonic) if file.endswith('.dat')]\n", "\n", "# Initialize an empty list to store the short wave down radiation data for each file\n", "all_data_swr = []\n", "\n", "all_data_swr_net=[]\n", "\n", "all_data_swr_up = []\n", "\n", "\n", "all_data_lwr_net=[]\n", "\n", "all_data_lwr_down = []\n", "\n", "all_data_lwr_up = []\n", "\n", "all_data_albedo=[] #albedo at surface, Swup/Swdown\n", "filtered_albedo = [] # albedo values when Swdown > 10\n", "\n", "\n", "# Loop through each file in the folder\n", "for file in data_files_sonic:\n", " # Construct the full file path\n", " file_path = os.path.join(data_folder_sonic, 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 data\n", " swr_column = 'SR15D1Dn_Irr' # Adjust column name if needed\n", " swr_up_column = 'SR15D1Up_Irr'\n", " swr_net_column = 'NetRs'\n", " # Select relevant columns containing long wave down radiation data\n", "\n", " lwr_net_column= 'NetRl'\n", " lwr_down_column= 'IR20Dn'\n", " lwr_up_column= 'IR20Up'\n", " \n", " \n", " albedo_column = 'Albedo'\n", " \n", " # Convert the column to numeric type\n", " data[swr_column] = pd.to_numeric(data[swr_column], errors='coerce')\n", " \n", " \n", " # Append the data to the list\n", " all_data_swr.append(data)\n", " \n", " # Plot and save the time series of short wave down radiation measurements for each file\n", " plt.figure(figsize=(10, 6))\n", " plt.plot(data['TIMESTAMP'], data[swr_column], label='Short Wave Down Radiation')\n", " plt.xlabel('Time')\n", " plt.ylabel('Short Wave Down Radiation')\n", " plt.title(f'Time Series of Short Wave Down Radiation - {file}')\n", " plt.legend()\n", " plt.grid(True)\n", " plt.xticks(rotation=45)\n", " plt.tight_layout()\n", " plt.savefig(os.path.join(data_folder_sonic, f'{os.path.splitext(file)[0]}_swr_plot.jpg'))\n", " plt.close()\n", " \n", " # Convert the column to numeric type \n", " data[lwr_net_column] = pd.to_numeric(data[lwr_net_column], errors='coerce')\n", "\n", " all_data_lwr_net.append(data)\n", " # Plot and save the time series of short wave down radiation measurements for each file\n", " plt.figure(figsize=(10, 6))\n", " plt.plot(data['TIMESTAMP'], data[lwr_net_column], label='Long Wave Net Radiation')\n", " plt.xlabel('Time')\n", " plt.ylabel('Long Wave Net Radiation, LW_Net[w/m^2]')\n", " plt.title(f'Time Series of Net Long Wave Radiation - {file}')\n", " plt.legend()\n", " plt.grid(True)\n", " plt.xticks(rotation=45)\n", " plt.tight_layout()\n", " plt.savefig(os.path.join(data_folder_sonic, f'{os.path.splitext(file)[0]}_lwr_net_plot.jpg'))\n", " plt.close()\n", " \n", " # Convert the column to numeric type \n", " data[lwr_down_column] = pd.to_numeric(data[lwr_down_column], errors='coerce')\n", "\n", " all_data_lwr_down.append(data)\n", " \n", " # Convert the column to numeric type \n", " data[lwr_up_column] = pd.to_numeric(data[lwr_up_column], errors='coerce')\n", "\n", " all_data_lwr_up.append(data)\n", " \n", " \n", " # Convert the column to numeric type \n", " data[swr_up_column] = pd.to_numeric(data[swr_up_column], errors='coerce')\n", "\n", " all_data_swr_up.append(data)\n", " \n", " # Convert the column to numeric type \n", " data[swr_net_column] = pd.to_numeric(data[swr_net_column], errors='coerce')\n", "\n", " all_data_swr_net.append(data)\n", " \n", " ## albedo \n", " data[albedo_column]= pd.to_numeric(data[albedo_column], errors='coerce')\n", " \n", " all_data_albedo.append(data)\n", " # Filter albedo values where short wave down radiation is larger than 10\n", " # Filter albedo values where short wave down radiation is larger than 10\n", " filtered_albedo_data = data[data[swr_column] > 10][[albedo_column, 'TIMESTAMP']]\n", " filtered_albedo.append(filtered_albedo_data)\n", "\n", "# Combine data from all files\n", "all_data_combined_swr = pd.concat(all_data_swr, ignore_index=True)\n", "\n", "all_data_combined_lwr_net = pd.concat(all_data_lwr_net, ignore_index=True)\n", "\n", "all_data_combined_lwr_down = pd.concat(all_data_lwr_down, ignore_index=True)\n", "\n", "all_data_combined_lwr_up = pd.concat(all_data_lwr_up, ignore_index=True)\n", "\n", "\n", "all_data_combined_swr_up = pd.concat(all_data_swr_up, ignore_index=True)\n", "\n", "all_data_combined_swr_net = pd.concat(all_data_swr_net, ignore_index=True)\n", "\n", "\n", "all_data_combined_albedo= pd.concat(all_data_albedo, ignore_index=True)\n", "filtered_albedo_combined = pd.concat(filtered_albedo, ignore_index=True)\n" ] }, { "cell_type": "code", "execution_count": null, "id": "dd27dfe7", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:51:22.863690Z", "iopub.status.busy": "2025-06-13T18:51:22.862680Z", "iopub.status.idle": "2025-06-13T18:51:22.883459Z", "shell.execute_reply": "2025-06-13T18:51:22.881430Z" }, "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "scrolled": true, "tags": [] }, "outputs": [], "source": [ "'''\n", "# Plot and save the combined time series of albedo measurements\n", "plt.figure(figsize=(10, 6))\n", "plt.plot(filtered_albedo_combined['TIMESTAMP'], 1/filtered_albedo_combined[albedo_column], label='Albedo')\n", "\n", "#plt.plot(all_data_combined_albedo['TIMESTAMP'], all_data_combined_swr_up[swr_up_column]/all_data_combined_swr[swr_column], label='Albedo calc')\n", "\n", "# Format the x-axis to show time in HH:MM format\n", "date_format = DateFormatter('%H:%M')\n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", "# Set y-axis limits from -0.1 to 1.1\n", "#plt.ylim(-1.1, 0.5)\n", "\n", "plt.xlabel('Time')\n", "plt.ylabel('Surface Albedo')\n", "plt.title('Combined Time Series of Surface ALbedo')\n", "plt.legend()\n", "plt.grid(True)\n", "plt.xticks(rotation=45)\n", "plt.tight_layout()\n", "plt.savefig(os.path.join(data_folder_sonic, 'combined_albedo_plot.pdf'),dpi=300)\n", "\n", "plt.show()\n", "\n", "# Plot and save the combined time series of short wave down radiation measurements\n", "plt.figure(figsize=(10, 6))\n", "plt.plot(all_data_combined_swr['TIMESTAMP'], all_data_combined_swr[swr_column], label='Short Wave Down Radiation')\n", "\n", "# Format the x-axis to show time in HH:MM format\n", "date_format = DateFormatter('%H:%M')\n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "plt.xlabel('Time')\n", "plt.ylabel('Short Wave Down Radiation, SW_dn [W/m^2]')\n", "plt.title('Combined Time Series of Short Wave Down Radiation [W/m^2]')\n", "plt.legend()\n", "plt.grid(True)\n", "plt.xticks(rotation=45)\n", "plt.tight_layout()\n", "plt.savefig(os.path.join(data_folder_sonic, 'combined_swr_plot.jpg'))\n", "plt.show()\n", "\n", "# Create a figure with four subplots arranged in a 2x2 grid\n", "fig, axs = plt.subplots(2, 2, figsize=(14, 12), sharex=True)\n", "\n", "# Plot net long wave radiation\n", "axs[0, 0].plot(all_data_combined_swr_net['TIMESTAMP'], all_data_combined_swr_net[swr_net_column], label='Net Short Wave Radiation')\n", "axs[0, 0].set_ylabel('Net SW Radiation [W/m^2]')\n", "axs[0, 0].set_title('Net Short Wave Radiation')\n", "axs[0, 0].legend()\n", "axs[0, 0].grid(True)\n", "\n", "# Plot long wave radiation up\n", "axs[0, 1].plot(all_data_combined_swr_up['TIMESTAMP'], all_data_combined_swr_up[swr_up_column], label='Short Wave Radiation Up, SW_up', color='orange')\n", "axs[0, 1].set_ylabel('SW Radiation Up [W/m^2]')\n", "axs[0, 1].set_title('Short Wave Radiation Up')\n", "axs[0, 1].legend()\n", "axs[0, 1].grid(True)\n", "\n", "# Plot long wave radiation down\n", "axs[1, 0].plot(all_data_combined_swr['TIMESTAMP'], all_data_combined_swr[swr_column], label='Short Wave Radiation Down, SW_dn', color='green')\n", "axs[1, 0].set_ylabel('SW Radiation Down [W/m^2]')\n", "axs[1, 0].set_xlabel('Time')\n", "axs[1, 0].set_title('Short Wave Radiation Down')\n", "axs[1, 0].legend()\n", "axs[1, 0].grid(True)\n", "\n", "# Plot all values combined\n", "axs[1, 1].plot(all_data_combined_swr_net['TIMESTAMP'], all_data_combined_swr_net[swr_net_column], label='Net Short Wave Radiation')\n", "axs[1, 1].plot(all_data_combined_swr_up['TIMESTAMP'], all_data_combined_swr_up[swr_up_column], label='Short Wave Radiation Up, SW_up', color='orange')\n", "axs[1, 1].plot(all_data_combined_swr['TIMESTAMP'], all_data_combined_swr[swr_column], label='Short Wave Radiation Down, SW_dn', color='green')\n", "axs[1, 1].set_xlabel('Time')\n", "axs[1, 1].set_ylabel('SW Radiation [W/m^2]')\n", "axs[1, 1].set_title('Combined Short Wave Radiation')\n", "axs[1, 1].legend()\n", "axs[1, 1].grid(True)\n", "\n", "# Format the x-axis to show time in HH:MM format\n", "date_format = DateFormatter('%H:%M')\n", "for ax in axs.flat:\n", " ax.xaxis.set_major_formatter(date_format)\n", " ax.tick_params(rotation=45)\n", "\n", "plt.tight_layout()\n", "plt.savefig(os.path.join(data_folder_sonic, 'combined_swr_plots.pdf'),dpi=300)\n", "plt.show()\n", "\n", "# Plot and save the combined time series of net long wave down radiation measurements\n", "plt.figure(figsize=(10, 6))\n", "plt.plot(all_data_combined_lwr_net['TIMESTAMP'], all_data_combined_lwr_net[lwr_net_column], label='Net Long Wave Radiation')\n", "\n", "plt.plot(all_data_combined_lwr_up['TIMESTAMP'], all_data_combined_lwr_up[lwr_up_column], label='Long Wave Radiation Up, LW_up')\n", "\n", "plt.plot(all_data_combined_lwr_down['TIMESTAMP'], all_data_combined_lwr_down[lwr_down_column], label='Long Wave Radiation Down, LW_dn')\n", "\n", "# Format the x-axis to show time in HH:MM format\n", "date_format = DateFormatter('%H:%M')\n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "plt.xlabel('Time')\n", "plt.ylabel('Net Long Wave Radiation, LW_net [W/m^2]')\n", "plt.title('Combined Time Series of Long Wave Radiation [W/m^2]')\n", "plt.legend()\n", "plt.grid(True)\n", "plt.xticks(rotation=45)\n", "plt.tight_layout()\n", "plt.savefig(os.path.join(data_folder_sonic, 'combined_lwr_plot.jpg'))\n", "plt.show()\n", "\n", "# Create a figure with four subplots arranged in a 2x2 grid\n", "fig, axs = plt.subplots(2, 2, figsize=(14, 12), sharex=True)\n", "\n", "# Plot net long wave radiation\n", "axs[0, 0].plot(all_data_combined_lwr_net['TIMESTAMP'], all_data_combined_lwr_net[lwr_net_column], label='Net Long Wave Radiation')\n", "axs[0, 0].set_ylabel('Net LW Radiation [W/m^2]')\n", "axs[0, 0].set_title('Net Long Wave Radiation')\n", "axs[0, 0].legend()\n", "axs[0, 0].grid(True)\n", "\n", "# Plot long wave radiation up\n", "axs[0, 1].plot(all_data_combined_lwr_up['TIMESTAMP'], all_data_combined_lwr_up[lwr_up_column], label='Long Wave Radiation Up, LW_up', color='orange')\n", "axs[0, 1].set_ylabel('LW Radiation Up [W/m^2]')\n", "axs[0, 1].set_title('Long Wave Radiation Up')\n", "axs[0, 1].legend()\n", "axs[0, 1].grid(True)\n", "\n", "# Plot long wave radiation down\n", "axs[1, 0].plot(all_data_combined_lwr_down['TIMESTAMP'], all_data_combined_lwr_down[lwr_down_column], label='Long Wave Radiation Down, LW_dn', color='green')\n", "axs[1, 0].set_ylabel('LW Radiation Down [W/m^2]')\n", "axs[1, 0].set_xlabel('Time')\n", "axs[1, 0].set_title('Long Wave Radiation Down')\n", "axs[1, 0].legend()\n", "axs[1, 0].grid(True)\n", "\n", "# Plot all values combined\n", "axs[1, 1].plot(all_data_combined_lwr_net['TIMESTAMP'], all_data_combined_lwr_net[lwr_net_column], label='Net Long Wave Radiation')\n", "axs[1, 1].plot(all_data_combined_lwr_up['TIMESTAMP'], all_data_combined_lwr_up[lwr_up_column], label='Long Wave Radiation Up, LW_up', color='orange')\n", "axs[1, 1].plot(all_data_combined_lwr_down['TIMESTAMP'], all_data_combined_lwr_down[lwr_down_column], label='Long Wave Radiation Down, LW_dn', color='green')\n", "axs[1, 1].set_xlabel('Time')\n", "axs[1, 1].set_ylabel('LW Radiation [W/m^2]')\n", "axs[1, 1].set_title('Combined Long Wave Radiation')\n", "axs[1, 1].legend()\n", "axs[1, 1].grid(True)\n", "\n", "# Format the x-axis to show time in HH:MM format\n", "date_format = DateFormatter('%H:%M')\n", "for ax in axs.flat:\n", " ax.xaxis.set_major_formatter(date_format)\n", " ax.tick_params(rotation=45)\n", "\n", "plt.tight_layout()\n", "plt.savefig(os.path.join(data_folder_sonic, 'combined_lwr_plots.pdf'),dpi=300)\n", "plt.show()\n", "'''" ] }, { "cell_type": "code", "execution_count": null, "id": "967bf774", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:51:23.053783Z", "iopub.status.busy": "2025-06-13T18:51:23.053783Z", "iopub.status.idle": "2025-06-13T18:51:46.136219Z", "shell.execute_reply": "2025-06-13T18:51:46.135100Z" }, "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "scrolled": true, "tags": [] }, "outputs": [], "source": [ "# Plot and save the combined time series of albedo measurements\n", "plt.figure(figsize=(10, 6))\n", "plt.plot(filtered_albedo_combined['TIMESTAMP'], 1/filtered_albedo_combined[albedo_column], \n", " label='Albedo', color='blue')\n", "\n", "# Format the x-axis to show time in HH:MM format\n", "date_format = DateFormatter('%H:%M')\n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", "# Set y-axis limits if needed (optional)\n", "# plt.ylim(-1.1, 0.5)\n", "\n", "# Set title and labels with larger font sizes\n", "plt.xlabel('Time', fontsize=16)\n", "plt.ylabel('Surface Albedo', fontsize=16)\n", "plt.title('Combined Time Series of Surface Albedo', fontsize=18)\n", "plt.legend(fontsize=14)\n", "plt.grid(True)\n", "plt.xticks(rotation=45, fontsize=12)\n", "plt.yticks(fontsize=12)\n", "\n", "plt.tight_layout()\n", "plt.savefig(os.path.join(data_folder_sonic, 'combined_albedo_plot.png'), dpi=300)\n", "plt.show()\n", "\n", "# Plot and save the combined time series of short wave down radiation measurements\n", "plt.figure(figsize=(10, 6))\n", "plt.plot(all_data_combined_swr['TIMESTAMP'], all_data_combined_swr[swr_column], \n", " label='Short Wave Down Radiation', color='blue')\n", "\n", "# Format the x-axis to show time in HH:MM format\n", "date_format = DateFormatter('%H:%M')\n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", "# Set title and labels with larger font sizes\n", "plt.xlabel('Time', fontsize=16)\n", "plt.ylabel('Short Wave Down Radiation, SW_dn [W/m²]', fontsize=16)\n", "plt.title('Combined Time Series of Short Wave Down Radiation [W/m²]', fontsize=18)\n", "plt.legend(fontsize=14)\n", "plt.grid(True)\n", "plt.xticks(rotation=45, fontsize=12)\n", "plt.yticks(fontsize=12)\n", "\n", "plt.tight_layout()\n", "plt.savefig(os.path.join(data_folder_sonic, 'combined_swr_plot.png'), dpi=300)\n", "plt.show()\n", "\n", "# Create a figure with four subplots arranged in a 2x2 grid\n", "fig, axs = plt.subplots(2, 2, figsize=(14, 12), sharex=True)\n", "\n", "# Plot net long wave radiation\n", "axs[0, 0].plot(all_data_combined_swr_net['TIMESTAMP'], all_data_combined_swr_net[swr_net_column], \n", " label='Net Short Wave Radiation', color='blue')\n", "axs[0, 0].set_ylabel('Net SW Radiation [W/m²]', fontsize=14)\n", "axs[0, 0].set_title('Net Short Wave Radiation', fontsize=16)\n", "axs[0, 0].legend(fontsize=12)\n", "axs[0, 0].grid(True)\n", "\n", "# Plot long wave radiation up\n", "axs[0, 1].plot(all_data_combined_swr_up['TIMESTAMP'], all_data_combined_swr_up[swr_up_column], \n", " label='Short Wave Radiation Up, SW_up', color='orange')\n", "axs[0, 1].set_ylabel('SW Radiation Up [W/m²]', fontsize=14)\n", "axs[0, 1].set_title('Short Wave Radiation Up', fontsize=16)\n", "axs[0, 1].legend(fontsize=12)\n", "axs[0, 1].grid(True)\n", "\n", "# Plot long wave radiation down\n", "axs[1, 0].plot(all_data_combined_swr['TIMESTAMP'], all_data_combined_swr[swr_column], \n", " label='Short Wave Radiation Down, SW_dn', color='green')\n", "axs[1, 0].set_ylabel('SW Radiation Down [W/m²]', fontsize=14)\n", "axs[1, 0].set_xlabel('Time', fontsize=16)\n", "axs[1, 0].set_title('Short Wave Radiation Down', fontsize=16)\n", "axs[1, 0].legend(fontsize=12)\n", "axs[1, 0].grid(True)\n", "\n", "# Plot all values combined\n", "axs[1, 1].plot(all_data_combined_swr_net['TIMESTAMP'], all_data_combined_swr_net[swr_net_column], \n", " label='Net Short Wave Radiation', color='blue')\n", "axs[1, 1].plot(all_data_combined_swr_up['TIMESTAMP'], all_data_combined_swr_up[swr_up_column], \n", " label='Short Wave Radiation Up, SW_up', color='orange')\n", "axs[1, 1].plot(all_data_combined_swr['TIMESTAMP'], all_data_combined_swr[swr_column], \n", " label='Short Wave Radiation Down, SW_dn', color='green')\n", "axs[1, 1].set_xlabel('Time', fontsize=16)\n", "axs[1, 1].set_ylabel('SW Radiation [W/m²]', fontsize=14)\n", "axs[1, 1].set_title('Combined Short Wave Radiation', fontsize=16)\n", "axs[1, 1].legend(fontsize=12)\n", "axs[1, 1].grid(True)\n", "\n", "# Format the x-axis to show time in HH:MM format\n", "date_format = DateFormatter('%H:%M')\n", "for ax in axs.flat:\n", " ax.xaxis.set_major_formatter(date_format)\n", " ax.tick_params(rotation=45)\n", "\n", "plt.tight_layout()\n", "plt.savefig(os.path.join(data_folder_sonic, 'combined_swr_plots.png'), dpi=300)\n", "plt.show()\n", "\n", "# Plot and save the combined time series of net long wave down radiation measurements\n", "plt.figure(figsize=(10, 6))\n", "plt.plot(all_data_combined_lwr_net['TIMESTAMP'], all_data_combined_lwr_net[lwr_net_column], \n", " label='Net Long Wave Radiation', color='blue')\n", "plt.plot(all_data_combined_lwr_up['TIMESTAMP'], all_data_combined_lwr_up[lwr_up_column], \n", " label='Long Wave Radiation Up, LW_up', color='orange')\n", "plt.plot(all_data_combined_lwr_down['TIMESTAMP'], all_data_combined_lwr_down[lwr_down_column], \n", " label='Long Wave Radiation Down, LW_dn', color='green')\n", "\n", "# Format the x-axis to show time in HH:MM format\n", "date_format = DateFormatter('%H:%M')\n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", "# Set title and labels with larger font sizes\n", "plt.xlabel('Time', fontsize=16)\n", "plt.ylabel('Net Long Wave Radiation, LW_net [W/m²]', fontsize=16)\n", "plt.title('Combined Time Series of Long Wave Radiation [W/m²]', fontsize=18)\n", "plt.legend(fontsize=14)\n", "plt.grid(True)\n", "plt.xticks(rotation=45, fontsize=12)\n", "plt.yticks(fontsize=12)\n", "\n", "plt.tight_layout()\n", "plt.savefig(os.path.join(data_folder_sonic, 'combined_lwr_plot.png'), dpi=300)\n", "plt.show()\n", "\n", "# Create a figure with four subplots arranged in a 2x2 grid for long wave radiation\n", "fig, axs = plt.subplots(2, 2, figsize=(14, 12), sharex=True)\n", "\n", "# Plot net long wave radiation\n", "axs[0, 0].plot(all_data_combined_lwr_net['TIMESTAMP'], all_data_combined_lwr_net[lwr_net_column], \n", " label='Net Long Wave Radiation', color='blue')\n", "axs[0, 0].set_ylabel('Net LW Radiation [W/m²]', fontsize=14)\n", "axs[0, 0].set_title('Net Long Wave Radiation', fontsize=16)\n", "axs[0, 0].legend(fontsize=12)\n", "axs[0, 0].grid(True)\n", "\n", "# Plot long wave radiation up\n", "axs[0, 1].plot(all_data_combined_lwr_up['TIMESTAMP'], all_data_combined_lwr_up[lwr_up_column], \n", " label='Long Wave Radiation Up, LW_up', color='orange')\n", "axs[0, 1].set_ylabel('LW Radiation Up [W/m²]', fontsize=14)\n", "axs[0, 1].set_title('Long Wave Radiation Up', fontsize=16)\n", "axs[0, 1].legend(fontsize=12)\n", "axs[0, 1].grid(True)\n", "\n", "# Plot long wave radiation down\n", "axs[1, 0].plot(all_data_combined_lwr_down['TIMESTAMP'], all_data_combined_lwr_down[lwr_down_column], \n", " label='Long Wave Radiation Down, LW_dn', color='green')\n", "axs[1, 0].set_ylabel('LW Radiation Down [W/m²]', fontsize=14)\n", "axs[1, 0].set_xlabel('Time', fontsize=16)\n", "axs[1, 0].set_title('Long Wave Radiation Down', fontsize=16)\n", "axs[1, 0].legend(fontsize=12)\n", "axs[1, 0].grid(True)\n", "\n", "# Plot all values combined\n", "axs[1, 1].plot(all_data_combined_lwr_net['TIMESTAMP'], all_data_combined_lwr_net[lwr_net_column], \n", " label='Net Long Wave Radiation', color='blue')\n", "axs[1, 1].plot(all_data_combined_lwr_up['TIMESTAMP'], all_data_combined_lwr_up[lwr_up_column], \n", " label='Long Wave Radiation Up, LW_up', color='orange')\n", "axs[1, 1].plot(all_data_combined_lwr_down['TIMESTAMP'], all_data_combined_lwr_down[lwr_down_column], \n", " label='Long Wave Radiation Down, LW_dn', color='green')\n", "axs[1, 1].set_xlabel('Time', fontsize=16)\n", "axs[1, 1].set_ylabel('LW Radiation [W/m²]', fontsize=14)\n", "axs[1, 1].set_title('Combined Long Wave Radiation', fontsize=16)\n", "axs[1, 1].legend(fontsize=12)\n", "axs[1, 1].grid(True)\n", "\n", "# Format the x-axis to show time in HH:MM format\n", "date_format = DateFormatter('%H:%M')\n", "for ax in axs.flat:\n", " ax.xaxis.set_major_formatter(date_format)\n", " ax.tick_params(rotation=45)\n", "\n", "plt.tight_layout()\n", "plt.savefig(os.path.join(data_folder_sonic, 'combined_lwr_plots.png'), dpi=300)\n", "plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "id": "1ee8e9cb", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:51:46.358917Z", "iopub.status.busy": "2025-06-13T18:51:46.358917Z", "iopub.status.idle": "2025-06-13T18:51:47.396169Z", "shell.execute_reply": "2025-06-13T18:51:47.396169Z" }, "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "# Calculate surface temperature from LWR up data\n", "all_data_combined_lwr_up['Surface_Temperature'] = (all_data_combined_lwr_up['IR20Up'] / sigma) ** 0.25\n", "#all_data_combined_lwr_up = pd.concat(all_data_lwr_up, ignore_index=True)\n", "\n", "# Plot and save the time series of surface temperature\n", "plt.figure(figsize=(10, 6))\n", "plt.plot(all_data_combined_lwr_up['TIMESTAMP'], all_data_combined_lwr_up['Surface_Temperature'], label='Surface Temperature')\n", "plt.xlabel('Time')\n", "plt.ylabel('Surface Temperature (K)')\n", "plt.title('Time Series of Surface Temperature')\n", "plt.legend()\n", "plt.grid(True)\n", "plt.xticks(rotation=45)\n", "plt.tight_layout()\n", "plt.savefig(os.path.join(data_folder_sonic, 'surface_temperature_plot.jpg'))\n", "plt.show()\n", "#print(all_data_combined_lwr_up)\n", "#print(\"Surface temperature calculation completed and plot saved.\")" ] }, { "cell_type": "code", "execution_count": null, "id": "6486ada6", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:51:47.651536Z", "iopub.status.busy": "2025-06-13T18:51:47.651536Z", "iopub.status.idle": "2025-06-13T18:51:48.597085Z", "shell.execute_reply": "2025-06-13T18:51:48.597085Z" }, "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "# Resample to 1-minute intervals and calculate the mean\n", "#all_data_combined_lwr_up = pd.concat(all_data_lwr_up, ignore_index=True)\n", "\n", "all_data_combined_lwr_up['TIMESTAMP'] = pd.to_datetime(all_data_combined_lwr_up['TIMESTAMP'], errors='coerce')\n", "all_data_combined_lwr_up.set_index('TIMESTAMP', inplace=True)\n", "\n", "\n", "#print(all_data_combined_lwr_up)\n", "resampled_surface_temp = all_data_combined_lwr_up['Surface_Temperature'].resample('1T').mean()\n", "print(resampled_surface_temp)\n", "\n", "\n", "# Plot both temperatures\n", "plt.figure(figsize=(10, 6))\n", "plt.plot(combined_data['TIMESTAMP'],combined_data['AirTC_E5567_Avg']+273.15, label='T_2m (K) from mast')\n", "plt.plot(resampled_surface_temp, label='T_srf from LW_up (K)') # Convert to °C\n", "plt.xlabel('Time')\n", "plt.ylabel('Temperature (°C)')\n", "plt.title('Time Series of T_2m and Surface Temperature')\n", "plt.legend()\n", "plt.grid(True)\n", "plt.xticks(rotation=45)\n", "plt.tight_layout()\n", "plt.savefig(os.path.join(data_folder_sonic, 'T2m_and_surface_temperature_plot.jpg'))\n", "plt.show()\n", "plt.close()" ] }, { "cell_type": "code", "execution_count": null, "id": "a0ca30f1", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:51:48.830369Z", "iopub.status.busy": "2025-06-13T18:51:48.830369Z", "iopub.status.idle": "2025-06-13T18:52:03.724648Z", "shell.execute_reply": "2025-06-13T18:52:03.724648Z" }, "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "tags": [] }, "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') #Ineichen model\n", " \n", " # Extract clear-sky GHI values\n", " # clear_sky_sw_down = clearsky['ghi']\n", " \n", " #return clear_sky_sw_down\n", "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", "# Example usage\n", "latitude = 52.6324 # Latitude of Alkmaar, Netherlands\n", "longitude = 4.7534 # Longitude of Alkmaar, Netherlands\n", "#timestamps = pd.date_range(start='2024-05-02 00:00:00', end='2024-05-02 23:59:59', freq='1S') # Frequency set to 1 second\n", "timestamps = pd.to_datetime(all_data_combined_swr['TIMESTAMP'], errors='coerce')# Use actual data timestamps\n", "# Convert timestamps to datetime index\n", "timestamps = pd.DatetimeIndex(timestamps)\n", "\n", "# Call the function to calculate clear-sky SW down radiation\n", "clear_sky_sw_down = calculate_clear_sky_sw_down_again(latitude, longitude, timestamps)\n", "\n", "\n", "# Plot the clear-sky SW down radiation\n", "plt.figure(figsize=(10, 6))\n", "# Plot measured short wave down radiation\n", "plt.plot(all_data_combined_swr['TIMESTAMP'], all_data_combined_swr[swr_column], \n", " label='Measured SW Down Radiation', color='blue')\n", "\n", "plt.plot(clear_sky_sw_down['TIMESTAMP'], clear_sky_sw_down['ghi'], \n", " label='Clear-Sky SW Down Radiation (pvlib)', color='red')\n", "\n", "# Format the x-axis to show time in HH:MM format\n", "date_format = DateFormatter('%H:%M')\n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", "# Set title and labels with larger font sizes\n", "plt.xlabel('Time', fontsize=16)\n", "plt.ylabel('SW Down Radiation (W/m²)', fontsize=16)\n", "plt.title('SW Down Radiation - May 11, 2024', fontsize=18)\n", "plt.legend(fontsize=14)\n", "plt.grid(True)\n", "plt.xticks(rotation=45, fontsize=12)\n", "plt.yticks(fontsize=12)\n", "\n", "plt.tight_layout()\n", "plt.savefig(os.path.join(data_folder_sonic, 'combined_swr_plot_with_clear_sky.jpg'))\n", "plt.show()\n", "\n", "# Plot the clear-sky GHI\n", "plt.figure(figsize=(10, 6))\n", "# Plot measured short wave down radiation\n", "plt.plot(all_data_combined_swr['TIMESTAMP'], all_data_combined_swr[swr_column], \n", " label='Measured SW Down Radiation', color='blue')\n", "\n", "plt.plot(clear_sky_sw_down['TIMESTAMP'], clear_sky_sw_down['ghi'], \n", " label='Clear-Sky SW Down Radiation (pvlib-Ineichen)', color='green')\n", "\n", "# Format the x-axis to show time in HH:MM format\n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "\n", "# Set title and labels with larger font sizes\n", "plt.xlabel('Time', fontsize=16)\n", "plt.ylabel('SW Down Radiation (W/m²)', fontsize=16)\n", "plt.title('SW Down Radiation - May 23, 2024', fontsize=18)\n", "plt.legend(fontsize=14)\n", "plt.grid(True)\n", "plt.xticks(rotation=45, fontsize=12)\n", "plt.yticks(fontsize=12)\n", "\n", "plt.tight_layout()\n", "#plt.savefig(os.path.join(data_folder_sonic, 'combined_swr_plot_with_clear_sky_ghi.jpg'))\n", "plt.show()\n", "\n", "\n", "\n", "# Step 3: Filter non-zero clear sky GHI values\n", "non_zero_clear_sky = clear_sky_sw_down[clear_sky_sw_down['ghi'] > 0]\n", "\n", "\n", "\n", "merged_data_csi= pd.merge(non_zero_clear_sky, all_data_combined_swr, on='TIMESTAMP', how='inner')\n", "\n", "# Calculate Clear Sky Index (CSI) using 'SR15D1Dn_Irr' and 'clear_sky_ghi'\n", "merged_data_csi['CSI'] = merged_data_csi['SR15D1Dn_Irr'] / merged_data_csi['ghi']\n", "#print(merged_data)\n", "\n", "# Plotting CSI against time\n", "plt.figure(figsize=(10, 6))\n", "plt.plot(merged_data_csi['TIMESTAMP'], merged_data_csi['CSI'], \n", " marker='o', linestyle='-', color='b', markersize=1)\n", "\n", "# Set title and labels with larger font sizes\n", "plt.title('Clear Sky Index (CSI) vs Time', fontsize=18)\n", "plt.gca().xaxis.set_major_formatter(date_format)\n", "plt.xlabel('Time', fontsize=16)\n", "plt.ylabel('Clear Sky Index (CSI)', fontsize=16)\n", "plt.ylim(0, 1.5) # Set y-axis limits from 0 to 1\n", "\n", "plt.grid(True)\n", "plt.xticks(rotation=45, fontsize=12)\n", "plt.yticks(fontsize=12)\n", "\n", "plt.tight_layout()\n", "#plt.savefig(os.path.join(data_folder_sonic, 'csi_plot.jpg'))\n", "plt.show()\n", " #Assuming all_data_combined_lwr_net, all_data_combined_lwr_up, all_data_combined_lwr_down, \n", "# combined_data, all_data_combined_swr, merged_data_csi, and all_qv_data are appropriately defined\n", "\n", "# Create a figure with a 2x3 grid of plots\n", "fig, axs = plt.subplots(2, 3, figsize=(18, 12))\n", "\n", "# Plot 1: Combined Long Wave Radiation\n", "axs[1, 0].plot(all_data_combined_lwr_net['TIMESTAMP'], all_data_combined_lwr_net[lwr_net_column], label='Net Long Wave Radiation')\n", "axs[1, 0].set_xlabel('Time')\n", "axs[1, 0].set_ylabel('LW Radiation [W/m^2]')\n", "axs[1, 0].set_title('Net Long Wave Radiation')\n", "axs[1, 0].legend()\n", "axs[1, 0].grid(True)\n", "axs[1, 0].tick_params(axis='x', rotation=45)\n", "\n", "# Plot 2: Long Wave Radiation Up\n", "axs[0, 0].plot(all_data_combined_lwr_net['TIMESTAMP'], all_data_combined_lwr_net['IR20Up'], label='Long Wave Radiation Up, LW_up', color='orange')\n", "axs[0, 0].set_ylabel('LW Radiation Up [W/m^2]')\n", "axs[0, 0].set_title('Long Wave Radiation Up')\n", "axs[0, 0].legend()\n", "axs[0, 0].grid(True)\n", "axs[0, 0].tick_params(axis='x', rotation=45)\n", "# Plot 3: Temperature Series\n", "axs[0, 1].plot(combined_data['TIMESTAMP'], combined_data['AirTC_E5567_Avg'] + 273.15, label='T_2m (K) from mast')\n", "axs[0, 1].plot(resampled_surface_temp.index, resampled_surface_temp.values, label='T_srf from LW_up (K)') # Convert to °C\n", "axs[0, 1].set_xlabel('Time')\n", "axs[0, 1].set_ylabel('Temperature (K)')\n", "axs[0, 1].set_title('Time Series of T_2m and Surface Temperature')\n", "axs[0, 1].legend()\n", "axs[0, 1].grid(True)\n", "axs[0, 1].tick_params(axis='x', rotation=45)\n", "\n", "# Plot 4: SW Down Radiation Comparison\n", "axs[1, 1].plot(all_data_combined_swr['TIMESTAMP'], all_data_combined_swr[swr_column], label='Measured SW Down Radiation', color='blue')\n", "axs[1, 1].plot(clear_sky_sw_down['TIMESTAMP'], clear_sky_sw_down['ghi'], label='Clear-Sky SW Down Radiation (pvlib-Ineichen)', color='green')\n", "axs[1, 1].set_xlabel('Time')\n", "axs[1, 1].set_ylabel('SW Down Radiation (W/m²)')\n", "axs[1, 1].set_title('SW Down Radiation Comparison')\n", "axs[1, 1].legend()\n", "axs[1, 1].grid(True)\n", "axs[1, 1].tick_params(axis='x', rotation=45)\n", "\n", "# Plot 5: Clear Sky Index (CSI)\n", "axs[0, 2].plot(merged_data_csi['TIMESTAMP'], merged_data_csi['CSI'], marker='o', linestyle='-', color='b', markersize=1)\n", "axs[0, 2].set_title('Clear Sky Index (CSI) vs Time')\n", "axs[0, 2].set_xlabel('Time')\n", "axs[0, 2].set_ylabel('Clear Sky Index (CSI)')\n", "axs[0, 2].set_ylim(0, 1.5) # Set y-axis limits from 0 to 1.5\n", "axs[0, 2].grid(True)\n", "axs[0, 2].tick_params(axis='x', rotation=45)\n", "\n", "# Plot 6: Specific Humidity Variation\n", "axs[1, 2].set_title('Specific Humidity Variation Over the Whole Day')\n", "axs[1, 2].set_xlabel('Time')\n", "axs[1, 2].set_ylabel('Specific Humidity (g/kg)')\n", "axs[1, 2].grid(True)\n", "axs[1, 2].tick_params(axis='x', rotation=45)\n", "\n", "for height in heights:\n", " axs[1, 2].plot(all_qv_data['TIMESTAMP'], all_qv_data[f'qv_{height}m'], label=f'Specific Humidity - Height: {height}m')\n", "\n", "axs[1, 2].legend()\n", "\n", "plt.tight_layout()\n", "\n", "# Save the plot\n", "plt.savefig(os.path.join(data_folder_sonic, 'first_combined_plots_grid.png'),dpi=300)\n", "\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "2f728d8f", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:52:04.053057Z", "iopub.status.busy": "2025-06-13T18:52:04.053057Z", "iopub.status.idle": "2025-06-13T18:52:06.711971Z", "shell.execute_reply": "2025-06-13T18:52:06.711971Z" }, "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "# (Keep your existing data-loading steps here, e.g. all_data_combined_swr and clear_sky_sw_down.)\n", "\n", "# Create a larger figure\n", "plt.figure(figsize=(8, 6))\n", "\n", "# Plot measured SW down radiation with a thicker line\n", "plt.plot(\n", " all_data_combined_swr['TIMESTAMP'],\n", " all_data_combined_swr[swr_column],\n", " label='Measured SW_dn',\n", " color='blue',\n", " linewidth=2.0 # thicker line\n", ")\n", "\n", "# Plot clear-sky model prediction with a thicker line\n", "plt.plot(\n", " clear_sky_sw_down['TIMESTAMP'],\n", " clear_sky_sw_down['ghi'],\n", " label='Clear‐Sky SW_dn \\n(pvlib‐Ineichen)',\n", " color='green',\n", " linewidth=2.0 # thicker line\n", ")\n", "\n", "# Get current Axes object\n", "ax = plt.gca()\n", "\n", "# 1) Thicken all four spines (axis lines)\n", "for spine in ax.spines.values():\n", " spine.set_linewidth(1.5)\n", "\n", "# 2) Increase tick size and thickness\n", "ax.tick_params(axis='both', which='major', labelsize=12, width=1.5, length=6)\n", "ax.tick_params(axis='both', which='minor', width=1.0, length=4)\n", "\n", "# 3) Format x‐axis to show hours:minutes\n", "date_format = mdates.DateFormatter('%H:%M')\n", "ax.xaxis.set_major_formatter(date_format)\n", "\n", "# (Optional) If you want to have a major tick every hour:\n", "#ax.xaxis.set_major_locator(mdates.HourLocator(interval=1))\n", "# And maybe a minor tick every 30 minutes:\n", "#ax.xaxis.set_minor_locator(mdates.MinuteLocator(interval=30))\n", "# 4) Format x‐axis to show time in HH:MM\n", "#date_format = mdates.DateFormatter('%H:%M')\n", "ax.xaxis.set_major_formatter(date_format)\n", "ax.xaxis.set_major_locator(mdates.HourLocator(interval=2))\n", "ax.xaxis.set_minor_locator(mdates.MinuteLocator(interval=30))\n", "# 4) Labels and title with larger font\n", "ax.set_xlabel('Time', fontsize=16, fontweight='bold')\n", "ax.set_ylabel('SW Down Radiation (W/m²)', fontsize=16, fontweight='bold')\n", "ax.set_title('SW Down Radiation', fontsize=18, fontweight='bold')\n", "\n", "# 5) Increase legend font & frame linewidth\n", "legend = ax.legend(fontsize=14, frameon=True,loc='upper right')\n", "legend.get_frame().set_linewidth(1.5)\n", "\n", "# 6) Thicker grid lines (optional: use linestyle='–' for dashed)\n", "ax.grid(True, which='major', linestyle='--', linewidth=0.8, alpha=0.7)\n", "ax.grid(True, which='minor', linestyle=':', linewidth=0.5, alpha=0.5)\n", "\n", "# 7) Rotate x‐tick labels for readability\n", "plt.xticks(rotation=45)\n", "\n", "# 8) Tight layout and save\n", "plt.tight_layout()\n", "output_path = os.path.join(data_folder_sonic, 'combined_swr_plot_with_clear_sky_ghi.jpg')\n", "plt.savefig(output_path, dpi=300) # save at high resolution for report\n", "plt.show()\n" ] }, { "cell_type": "code", "execution_count": null, "id": "b3feedd7", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:52:07.009141Z", "iopub.status.busy": "2025-06-13T18:52:07.009141Z", "iopub.status.idle": "2025-06-13T18:52:08.355087Z", "shell.execute_reply": "2025-06-13T18:52:08.353943Z" }, "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "scrolled": false, "tags": [] }, "outputs": [], "source": [ "# (Assume merged_data_csi and date_format are already defined.)\n", "\n", "plt.figure(figsize=(8, 6))\n", "\n", "ax = plt.gca()\n", "\n", "# 1) Plot CSI with bolder line and slightly larger markers\n", "ax.plot(\n", " merged_data_csi['TIMESTAMP'],\n", " merged_data_csi['CSI'],\n", " #marker='o',\n", " linestyle='-',\n", " color='b',\n", " #markersize=4, # slightly larger markers\n", " linewidth=1.5, # bolder line\n", " label='CSI'\n", ")\n", "\n", "# 2) Thicken spines (axis lines)\n", "for spine in ax.spines.values():\n", " spine.set_linewidth(1.5)\n", "\n", "# 3) Tick parameters for major & minor ticks\n", "ax.tick_params(axis='both', which='major', labelsize=12, width=1.5, length=6)\n", "ax.tick_params(axis='both', which='minor', width=1.0, length=4)\n", "\n", "# 4) Format x‐axis with HH:MM\n", "ax.xaxis.set_major_formatter(date_format)\n", "\n", "# Optionally, place a major tick every hour and minor tick every 30 minutes:\n", "ax.xaxis.set_major_locator(mdates.HourLocator(interval=2))\n", "ax.xaxis.set_minor_locator(mdates.MinuteLocator(interval=30))\n", "\n", "# 5) Labels and title with bold font\n", "ax.set_title('Clear Sky Index (CSI) vs. Time', fontsize=18, fontweight='bold')\n", "ax.set_xlabel('Time', fontsize=16, fontweight='bold')\n", "ax.set_ylabel('Clear Sky Index (CSI)', fontsize=16, fontweight='bold')\n", "\n", "# 6) Y-axis limits and tick font size\n", "ax.set_ylim(0, 1.5)\n", "ax.set_yticks([0.0, 0.5, 1.0, 1.5]) # explicit ticks\n", "# (If you want minor ticks on the y-axis, you can add: ax.yaxis.set_minor_locator(...) )\n", "ax.tick_params(axis='y', labelsize=12, width=1.5, length=6)\n", "\n", "# 7) Grid styling: dashed major, dotted minor\n", "ax.grid(True, which='major', linestyle='--', linewidth=0.8, alpha=0.7)\n", "ax.grid(True, which='minor', linestyle='', linewidth=0.5, alpha=0.5)\n", "\n", "# 8) Rotate x‐tick labels for readability\n", "plt.xticks(rotation=45)\n", "\n", "# 9) (Optional) If you want a legend\n", "# legend = ax.legend(fontsize=14, frameon=True)\n", "# legend.get_frame().set_linewidth(1.5)\n", "\n", "# 10) Tight layout and save\n", "plt.tight_layout()\n", "output_path = os.path.join(data_folder_sonic, 'csi_plot_report_style.jpg')\n", "plt.savefig(output_path, dpi=300)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "c3d09444", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:52:08.655327Z", "iopub.status.busy": "2025-06-13T18:52:08.655327Z", "iopub.status.idle": "2025-06-13T18:52:08.712876Z", "shell.execute_reply": "2025-06-13T18:52:08.712876Z" }, "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "print(merged_data_csi)" ] }, { "cell_type": "code", "execution_count": null, "id": "40f13fef", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:52:08.987658Z", "iopub.status.busy": "2025-06-13T18:52:08.987658Z", "iopub.status.idle": "2025-06-13T18:52:09.167519Z", "shell.execute_reply": "2025-06-13T18:52:09.167519Z" }, "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "# Convert TIMESTAMP to datetime if it's not already\n", "all_data_combined_albedo['TIMESTAMP'] = pd.to_datetime(all_data_combined_albedo['TIMESTAMP'], errors='coerce')\n", "merged_data_csi['TIMESTAMP'] = pd.to_datetime(merged_data_csi['TIMESTAMP'], errors='coerce')\n", "\n", "# 2. Set the TIMESTAMP column as the index for both DataFrames\n", "all_data_combined_albedo.set_index('TIMESTAMP', inplace=True)\n", "merged_data_csi.set_index('TIMESTAMP', inplace=True)\n", "\n", "# 3. Extract only the CSI column from merged_data_csi\n", "csi_series = merged_data_csi[['CSI']]\n", "\n", "# 4. Merge the CSI series into all_data_combined_albedo on the timestamp index, using a left join\n", "all_data_combined_albedo = all_data_combined_albedo.merge(\n", " csi_series,\n", " left_index=True,\n", " right_index=True,\n", " how='left'\n", ")\n", "\n", "# Select only numeric columns for resampling\n", "numeric_columns = all_data_combined_albedo.select_dtypes(include='number').columns\n", "all_data_combined_albedo_numeric = all_data_combined_albedo[numeric_columns]\n", "\n", "# Resample the data to 10-minute intervals and calculate the mean for each interval\n", "all_data_combined_albedo_10min = all_data_combined_albedo_numeric.resample('10T').mean()\n", "\n", "# Calculate net radiation separately\n", "net_radiation_10min = (\n", " (all_data_combined_albedo[swr_up_column] - all_data_combined_albedo[swr_column])\n", " .resample('10T')\n", " .mean() +\n", " (all_data_combined_albedo[lwr_up_column] - all_data_combined_albedo[lwr_down_column])\n", " .resample('10T')\n", " .mean()\n", ")\n", "\n", "# Combine net radiation with the resampled data\n", "all_data_combined_albedo_10min['Net_Radiation_10min'] = net_radiation_10min\n", "\n", "# Reset the index to make TIMESTAMP a column again\n", "all_data_combined_albedo_10min.reset_index(inplace=True)\n", "\n", "# Display the resulting DataFrame\n", "print(all_data_combined_albedo_10min)\n" ] }, { "cell_type": "markdown", "id": "e9af7012", "metadata": { "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "tags": [] }, "source": [ "### Sonic" ] }, { "cell_type": "code", "execution_count": null, "id": "a90b8b70", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:52:09.732199Z", "iopub.status.busy": "2025-06-13T18:52:09.732199Z", "iopub.status.idle": "2025-06-13T18:52:36.641455Z", "shell.execute_reply": "2025-06-13T18:52:36.641455Z" }, "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "# Directory containing Sonic data\n", "# ⚠️ Edit this path to match your local setup before running.\n", "data_dir = f\"C:\\\\path\\\\to\\\\your\\\\Sonic\\\\{month_str}\\\\{date_str}\"\n", "\n", "\n", "# Initialize lists to store wT fluxes and corresponding timestamps\n", "timestamps_10min = []\n", "wT_fluxes_10min = []\n", "wrhoqv_fluxes_10min = []\n", "wrhoco2_fluxes_10min = []\n", "average_temperatures_10min = []\n", "average_temperatures_corr_10min = []\n", "average_h2o_density_10min = []\n", "average_co2_density_10min = []\n", "average_wind_ux_10min = []\n", "average_wind_uy_10min = []\n", "average_wind_uz_10min = []\n", "uw_flux_10min = []\n", "vw_flux_10min = []\n", "uv_flux_10min = []\n", "average_wind_ux_10min_corr=[]\n", "average_wind_uy_10min_corr=[]\n", "uw_flux_10min_corr=[]\n", "vw_flux_10min_corr=[]\n", "uv_flux_10min_corr=[]\n", "\n", "for filename in os.listdir(data_dir):\n", " if filename.endswith('.dat'):\n", " file_path = os.path.join(data_dir, filename)\n", " \n", " # Load data from .dat file, skipping the first row\n", " data = pd.read_csv(file_path, skiprows=1, delimiter=',', encoding='latin1', low_memory=False)\n", " \n", " # Drop any rows with missing values\n", " data.dropna(inplace=True)\n", " \n", " # Convert the timestamp to datetime format\n", " data['TIMESTAMP'] = pd.to_datetime(data['TIMESTAMP'], errors='coerce')\n", " data['Ux'] = pd.to_numeric(data['Ux'], errors='coerce')\n", " data['Uy'] = pd.to_numeric(data['Uy'], errors='coerce')\n", " data['Uz'] = pd.to_numeric(data['Uz'], errors='coerce')\n", " data['T_SONIC'] = pd.to_numeric(data['T_SONIC'], errors='coerce')\n", " data['T_SONIC_corr'] = pd.to_numeric(data['T_SONIC_corr'], errors='coerce')\n", " data['H2O_density'] = pd.to_numeric(data['H2O_density'], errors='coerce')\n", " data['CO2_density'] = pd.to_numeric(data['CO2_density'], errors='coerce')\n", "\n", " # Apply directional correction (rotate 116° clockwise)\n", " theta_deg = 116\n", " theta_rad = np.deg2rad(theta_deg)\n", "\n", " # Correct Ux and Uy (rotate into mast-aligned frame)\n", " #u_rot = data['Ux'] * np.cos(theta_rad) + data['Uy'] * np.sin(theta_rad)\n", " #v_rot = -data['Ux'] * np.sin(theta_rad) + data['Uy'] * np.cos(theta_rad)\n", "\n", " # Rotate Ux and Uy into corrected coordinate system\n", " data['Ux_corrected'] = data['Ux'] * np.cos(theta_rad) + data['Uy'] * np.sin(theta_rad)\n", " data['Uy_corrected'] = -data['Ux'] * np.sin(theta_rad) + data['Uy'] * np.cos(theta_rad)\n", "\n", "\n", "# Loop through each 10-minute chunk in the file\n", " for i in range(0, len(data), 12000): # 12000 samples in 10 minutes with 20 Hz frequency\n", " ten_minute_data = data.iloc[i:i+12000] # Get data for 10 minutes\n", " \n", " if not ten_minute_data.empty:\n", " # Calculate 10-minute means\n", " mean_wind_ux = ten_minute_data['Ux'].mean()\n", " mean_wind_uy = ten_minute_data['Uy'].mean()\n", " mean_wind_uz = ten_minute_data['Uz'].mean()\n", " mean_wind_ux_corr = ten_minute_data['Ux_corrected'].mean()\n", " mean_wind_uy_corr = ten_minute_data['Uy_corrected'].mean()\n", " # Calculate the mean temperature for the 10 minutes\n", " mean_T_sonic = ten_minute_data['T_SONIC'].mean()\n", " mean_T_sonic_corr = ten_minute_data['T_SONIC_corr'].mean()\n", " mean_h2o_density = ten_minute_data['H2O_density'].mean()\n", " mean_co2_density = ten_minute_data['CO2_density'].mean()\n", "\n", " # Calculate perturbations for wind components\n", " u_prime = ten_minute_data['Ux'] - mean_wind_ux\n", " v_prime = ten_minute_data['Uy'] - mean_wind_uy\n", " w_prime = ten_minute_data['Uz'] - mean_wind_uz\n", " u_prime_corr = ten_minute_data['Ux_corrected'] - mean_wind_ux_corr\n", " v_prime_corr = ten_minute_data['Uy_corrected'] - mean_wind_uy_corr\n", "\n", " # Calculate momentum fluxes\n", " uw_flux = u_prime * w_prime\n", " vw_flux = v_prime * w_prime\n", " uv_flux = u_prime * v_prime\n", " uw_flux_corr = u_prime_corr * w_prime\n", " vw_flux_corr = v_prime_corr * w_prime\n", " uv_flux_corr = u_prime_corr * v_prime_corr\n", " \n", " t_sonic_prime = ten_minute_data['T_SONIC_corr'] - mean_T_sonic_corr\n", " \n", " # Calculate wT flux (w'T' flux)\n", " wT_flux = (w_prime * (t_sonic_prime + 273.15)).mean()\n", " \n", " # w'qv'flux\n", " rhoqv_prime = ten_minute_data['H2O_density'] - mean_h2o_density\n", " wrhoqv_flux = (w_prime * rhoqv_prime).mean()\n", " \n", " # rhoCo2'\n", " rhoco2_prime = ten_minute_data['CO2_density'] - mean_co2_density\n", " # w'rhoco2'\n", " wrhoco2_flux = (rhoco2_prime * w_prime).mean()\n", " \n", " # Store data in lists\n", " timestamps_10min.append(ten_minute_data['TIMESTAMP'].iloc[0].floor('10T'))\n", " wT_fluxes_10min.append(wT_flux)\n", " wrhoqv_fluxes_10min.append(wrhoqv_flux)\n", " wrhoco2_fluxes_10min.append(wrhoco2_flux)\n", " average_temperatures_10min.append(mean_T_sonic + 273.15)\n", " average_temperatures_corr_10min.append(mean_T_sonic_corr + 273.15)\n", " average_h2o_density_10min.append(mean_h2o_density)\n", " average_co2_density_10min.append(mean_co2_density)\n", " average_wind_ux_10min.append(mean_wind_ux)\n", " average_wind_uy_10min.append(mean_wind_uy)\n", " average_wind_uz_10min.append(mean_wind_uz)\n", " average_wind_ux_10min_corr.append(mean_wind_ux_corr)\n", " average_wind_uy_10min_corr.append(mean_wind_uy_corr)\n", " uw_flux_10min.append(uw_flux.mean())\n", " vw_flux_10min.append(vw_flux.mean())\n", " uv_flux_10min.append(uv_flux.mean())\n", " uw_flux_10min_corr.append(uw_flux_corr.mean())\n", " vw_flux_10min_corr.append(vw_flux_corr.mean())\n", " uv_flux_10min_corr.append(uv_flux_corr.mean())\n", "# Create DataFrame for the collected data\n", "flux_data_10min = pd.DataFrame({\n", " 'TIMESTAMP': timestamps_10min,\n", " 'wT_Flux': wT_fluxes_10min,\n", " 'wrhoqv_Flux': wrhoqv_fluxes_10min,\n", " 'Average_Temperature': average_temperatures_10min,\n", " 'Average_Temperature_Corr': average_temperatures_corr_10min,\n", " 'Average_H2O_Density': average_h2o_density_10min,\n", " 'Average_CO2_Density': average_co2_density_10min,\n", " 'wrhoCO2_Flux': wrhoco2_fluxes_10min,\n", " 'Average_Wind_Ux': average_wind_ux_10min,\n", " 'Average_Wind_Uy': average_wind_uy_10min,\n", " 'Average_Wind_Uz': average_wind_uz_10min,\n", " 'Average_Wind_Ux_corrected': average_wind_ux_10min_corr,\n", " 'Average_Wind_Uy_corrected': average_wind_uy_10min_corr,\n", "\n", " 'uw_flux': uw_flux_10min,\n", " 'vw_flux': vw_flux_10min,\n", " 'uv_flux': uv_flux_10min,\n", " 'uw_flux_corr': uw_flux_10min_corr,\n", " 'vw_flux_corr': vw_flux_10min_corr,\n", " 'uv_flux_corr': uv_flux_10min_corr\n", "})\n", "\n", "print(flux_data_10min)" ] }, { "cell_type": "code", "execution_count": null, "id": "29db78d2", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:52:36.941312Z", "iopub.status.busy": "2025-06-13T18:52:36.941312Z", "iopub.status.idle": "2025-06-13T18:52:36.956579Z", "shell.execute_reply": "2025-06-13T18:52:36.956579Z" }, "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "scrolled": true, "tags": [] }, "outputs": [], "source": [ "print(flux_data_10min.columns)" ] }, { "cell_type": "code", "execution_count": null, "id": "eafb4783", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:52:37.243155Z", "iopub.status.busy": "2025-06-13T18:52:37.243155Z", "iopub.status.idle": "2025-06-13T18:52:37.346284Z", "shell.execute_reply": "2025-06-13T18:52:37.346284Z" }, "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "scrolled": false, "tags": [] }, "outputs": [], "source": [ "# Merge the resampled net radiation with flux_data_10min\n", "#merged_data_10min = pd.merge_asof(flux_data_10min.sort_values('TIMESTAMP'),\n", " # all_data_combined_albedo_10min.sort_values('TIMESTAMP'),\n", " # on='TIMESTAMP', direction='nearest')\n", "merged_data_10min = pd.merge(\n", " flux_data_10min,\n", " all_data_combined_albedo_10min,\n", " on='TIMESTAMP', # Only merge on exact timestamp matches\n", " how='inner' # Use 'inner' to keep only rows with matching timestamps in both dataframes\n", ")\n", "\n", "# Merge dataframes based on overlapping timestamps\n", "#merged_data_10min = pd.merge_asof(merged_data_10min.sort_values('TIMESTAMP'),\n", " # all_pressure_data_10min.sort_values('TIMESTAMP'),\n", " # on='TIMESTAMP',\n", " #direction='nearest')\n", " \n", "merged_data_10min = pd.merge(\n", " merged_data_10min,\n", " all_pressure_data_10min,\n", " on='TIMESTAMP', # Only merge on exact timestamp matches\n", " how='inner' # Use 'inner' to keep only rows with matching timestamps in both dataframes\n", ")\n", "\n", "#Calculations\n", "merged_data_10min['tau_xz'] = -merged_data_10min['rho_air_Tv'] * merged_data_10min['uw_flux']\n", "merged_data_10min['tau_yz'] = -merged_data_10min['rho_air_Tv'] * merged_data_10min['vw_flux']\n", "merged_data_10min['tau_xy'] = -merged_data_10min['rho_air_Tv'] * merged_data_10min['uv_flux']\n", "merged_data_10min['tau_xz_corr'] = -merged_data_10min['rho_air_Tv'] * merged_data_10min['uw_flux_corr']\n", "merged_data_10min['tau_yz_corr'] = -merged_data_10min['rho_air_Tv'] * merged_data_10min['vw_flux_corr']\n", "merged_data_10min['tau_xy_corr'] = -merged_data_10min['rho_air_Tv'] * merged_data_10min['uv_flux_corr']\n", "\n", "# Calculate sensible heat flux (SHF)\n", "merged_data_10min['SHF'] = merged_data_10min['wT_Flux'] * merged_data_10min['rho_air_Tv'] * Cp\n", "\n", "# Calculate latent heat flux (LHF)\n", "merged_data_10min['qv_sonic'] = merged_data_10min['Average_H2O_Density'] / merged_data_10min['rho_air_Tv']\n", "merged_data_10min['wqv_Flux'] = merged_data_10min['wrhoqv_Flux'] / merged_data_10min['rho_air_Tv']\n", "merged_data_10min['LHF'] = merged_data_10min['wqv_Flux'] * merged_data_10min['rho_air_Tv'] * Lv / 1000\n", "\n", "## Calculate ground heat flux (G)\n", "merged_data_10min['G'] = -(merged_data_10min['SHF'] + merged_data_10min['LHF'] + merged_data_10min['Net_Radiation_10min'])\n", "\n", "# Calculate wind speed\n", "merged_data_10min['Wind_Speed'] = np.sqrt(\n", " merged_data_10min['Average_Wind_Ux']**2 +\n", " merged_data_10min['Average_Wind_Uy']**2 +\n", " merged_data_10min['Average_Wind_Uz']**2\n", ")\n", "\n", "\n", "merged_data_10min['PPM_CO2'] = merged_data_10min['Average_CO2_Density'] * m_atm / (m_co2 * merged_data_10min['rho_air_Tv'])\n", "merged_data_10min['F_CO2'] = merged_data_10min['wrhoCO2_Flux'] * m_atm / (m_co2 * merged_data_10min['rho_air_Tv'])\n", "\n", "\n", "# Display the merged dataframe\n", "print(merged_data_10min)\n", "\n", "# Define the file path where you want to save the CSV within the data directory\n", "output_file = os.path.join(data_dir, 'merged_data_10min.csv')\n", "\n", "# Save merged_data_10min to CSV\n", "merged_data_10min.to_csv(output_file, index=False)\n", "\n", "print(f\"Saved merged data to {output_file}\")" ] }, { "cell_type": "code", "execution_count": null, "id": "c49a2051", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:52:37.624908Z", "iopub.status.busy": "2025-06-13T18:52:37.609239Z", "iopub.status.idle": "2025-06-13T18:52:37.632547Z", "shell.execute_reply": "2025-06-13T18:52:37.632547Z" }, "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "scrolled": true, "tags": [] }, "outputs": [], "source": [ "print(merged_data_10min.columns)" ] }, { "cell_type": "code", "execution_count": null, "id": "cf1578f1", "metadata": { "execution": { "iopub.execute_input": "2025-06-13T18:52:37.990373Z", "iopub.status.busy": "2025-06-13T18:52:37.990373Z", "iopub.status.idle": "2025-06-13T18:52:42.055080Z", "shell.execute_reply": "2025-06-13T18:52:42.055080Z" }, "papermill": { "duration": null, "end_time": null, "exception": null, "start_time": null, "status": "completed" }, "tags": [] }, "outputs": [], "source": [ "# Create a figure with 2x2 subplots\n", "fig, axs = plt.subplots(2, 2, figsize=(14, 10))\n", "\n", "# Plot 1: SHF and LHF vs Time\n", "axs[0, 0].scatter(merged_data_10min['TIMESTAMP'], merged_data_10min['SHF'], s=10, alpha=0.7, color='black', label='SHF')\n", "axs[0, 0].scatter(merged_data_10min['TIMESTAMP'], merged_data_10min['LHF'], s=10, alpha=0.7, color='green',label='LHF')\n", "axs[0, 0].xaxis.set_major_formatter(DateFormatter('%H:%M'))\n", "axs[0, 0].set_xlabel('Time')\n", "axs[0, 0].set_ylabel('Heat Flux (W/m^2)')\n", "axs[0, 0].set_title('SHF, LHF vs Time (averaged in 10 min intervals)')\n", "axs[0, 0].legend()\n", "axs[0, 0].grid(True)\n", "\n", "# Plot 2: G vs Time\n", "axs[0, 1].scatter(merged_data_10min['TIMESTAMP'], merged_data_10min['G'], color='red', s=10, alpha=0.7)\n", "axs[0, 1].xaxis.set_major_formatter(DateFormatter('%H:%M'))\n", "axs[0, 1].set_xlabel('Time')\n", "axs[0, 1].set_ylabel('G (W/m^2)')\n", "axs[0, 1].set_title('Heat flux into the ground vs Time (averaged in 10 min intervals)')\n", "axs[0, 1].grid(True)\n", "\n", "# Plot 3: F_CO2 vs Time\n", "axs[1, 0].scatter(merged_data_10min['TIMESTAMP'], merged_data_10min['F_CO2'], s=10, color='indigo', alpha=0.7)\n", "axs[1, 0].xaxis.set_major_formatter(DateFormatter('%H:%M'))\n", "axs[1, 0].set_xlabel('Time')\n", "axs[1, 0].set_ylabel('F_CO2 (ppm*ms^-1)')\n", "axs[1, 0].set_title('CO2 flux vs Time (averaged in 10 min intervals)')\n", "axs[1, 0].grid(True)\n", "axs[1, 0].set_ylim(-0.5, 0.5)\n", "\n", "# Plot 4: Net Radiation vs Time (assuming you have this data)\n", "axs[1, 1].scatter(merged_data_10min['TIMESTAMP'], merged_data_10min['Net_Radiation_10min'], color='blue', s=10, alpha=0.7)\n", "axs[1, 1].xaxis.set_major_formatter(DateFormatter('%H:%M'))\n", "axs[1, 1].set_xlabel('Time')\n", "axs[1, 1].set_ylabel('Net Radiation (W/m^2)')\n", "axs[1, 1].set_title('Net Radiation vs Time')\n", "axs[1, 1].grid(True)\n", "\n", "# Adjust layout to prevent overlap\n", "fig.tight_layout()\n", "\n", "# Save the figure as a PNG file in the data directory\n", "plt.savefig(os.path.join(data_dir, 'combined_plots_10min.png'), dpi=300)\n", "\n", "# Show the plot\n", "plt.show()" ] } ], "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" }, "papermill": { "default_parameters": {}, "duration": 29.18899, "end_time": "2025-06-13T18:53:12.882908", "environment_variables": {}, "exception": null, "input_path": "C:\\Users\\magda\\Master_Thesis\\PMT-profiles-fast.ipynb", "output_path": "C:\\Users\\magda\\Master_Thesis\\PMT-profiles-fast.ipynb", "parameters": { "date_str": "2024-03-15" }, "start_time": "2025-06-13T18:52:43.693918", "version": "2.6.0" } }, "nbformat": 4, "nbformat_minor": 5 }