From e5c920c2d3a9a3afc9f63423f30b3a958195298b Mon Sep 17 00:00:00 2001 From: patrick Date: Thu, 22 Sep 2022 16:21:09 +0200 Subject: [PATCH 1/6] Enh: added forecast comparaison parameters to env analysis --- rocketpy/EnvironmentAnalysis.py | 33 +++++++++++++++++++++++++++++++++ 1 file changed, 33 insertions(+) diff --git a/rocketpy/EnvironmentAnalysis.py b/rocketpy/EnvironmentAnalysis.py index 5a1447f5b..6e78c4f19 100644 --- a/rocketpy/EnvironmentAnalysis.py +++ b/rocketpy/EnvironmentAnalysis.py @@ -5,6 +5,7 @@ __license__ = "MIT" import bisect +import datetime import json import warnings from collections import defaultdict @@ -22,6 +23,7 @@ from scipy import stats from timezonefinder import TimezoneFinder from windrose import WindroseAxes +from rocketpy.Environment import Environment from rocketpy.Function import Function from rocketpy.units import convert_units @@ -68,6 +70,8 @@ def __init__( timezone=None, unit_system="metric", load_previous_data=None, + forecast_comparaison=False, + forecast_date=None, ): """Constructor for the EnvironmentAnalysis class. Parameters @@ -203,6 +207,28 @@ def __init__( # Run calculations self.process_data() + # Processing forecast + self.forecast = None + if forecast_comparaison: + self.forecast = {} + hours = list(self.pressureLevelDataDict.values())[0].keys() + for hour in hours: + hour_datetime = datetime.datetime( + year=forecast_date.year, month=forecast_date.month, day=forecast_date.day, + hour=int(hour)) + + Env = Environment( + railLength=5, + date=hour_datetime, + latitude=self.latitude, + longitude=self.longitude, + elevation=self.elevation, + ) + Env.setAtmosphericModel(type="Ensemble", file="GEFS") + Env.selectEnsembleMember(1) + self.forecast[hour] = Env + + def __bilinear_interpolation(self, x, y, x1, x2, y1, y2, z11, z12, z21, z22): """ Bilinear interpolation. @@ -2520,6 +2546,13 @@ def plot_wind_profile_over_average_day(self, SAcup_altitude_constraints=False): x_max = current_x_max if current_x_max > x_max else x_max y_max = current_y_max if current_y_max > y_max else y_max y_min = current_y_min if current_y_min < y_min else y_min + + if self.forecast: + forecast = self.forecast + x = self.average_wind_profile_at_given_hour[hour][0] + y = forecast[hour].windSpeed.getValue(x) + ax.plot(x, y, "b--") + ax.label_outer() ax.grid() # Set x and y limits for the last axis. Since axes are shared, set to all From 4e40c42398d52c756acd28db19c98da29968b25a Mon Sep 17 00:00:00 2001 From: Gui-FernandesBR Date: Tue, 27 Sep 2022 16:23:42 +0200 Subject: [PATCH 2/6] ENH: drafting compareEnvironments method --- rocketpy/EnvironmentAnalysis.py | 5 +- rocketpy/utilities.py | 134 ++++++++++++++++++++++++++++++++ 2 files changed, 136 insertions(+), 3 deletions(-) diff --git a/rocketpy/EnvironmentAnalysis.py b/rocketpy/EnvironmentAnalysis.py index e9bc97d35..c94196566 100644 --- a/rocketpy/EnvironmentAnalysis.py +++ b/rocketpy/EnvironmentAnalysis.py @@ -107,8 +107,7 @@ def __init__( pressureLevelDataFile=None, timezone=None, unit_system="metric", - load_previous_data=None, - forecast_comparaison=False, + forecast_comparison=False, forecast_date=None, maxExpectedAltitude=None, ): @@ -213,7 +212,7 @@ def __init__( # Processing forecast self.forecast = None - if forecast_comparaison: + if forecast_comparison: self.forecast = {} hours = list(self.pressureLevelDataDict.values())[0].keys() for hour in hours: diff --git a/rocketpy/utilities.py b/rocketpy/utilities.py index dc847ea3c..7bf510969 100644 --- a/rocketpy/utilities.py +++ b/rocketpy/utilities.py @@ -3,10 +3,12 @@ __copyright__ = "Copyright 20XX, RocketPy Team" __license__ = "MIT" +import matplotlib.pyplot as plt import numpy as np from scipy.integrate import solve_ivp from .Environment import Environment +from .EnvironmentAnalysis import EnvironmentAnalysis from .Function import Function @@ -197,3 +199,135 @@ def du(z, u): velocityFunction() return altitudeFunction, velocityFunction, final_sol + + +def compareEnvironments(env1, env2, names=["env1", "env2"]): + """Compares two environments and plots everything. + Useful when comparing Environment Analysis results against forecasts. + + Parameters + ---------- + env1 : Environment or EnvironmentAnalysis + Environment to compare with. + env2: Environment or EnvironmentAnalysis + Environment to compare with. + + Returns + ------- + diff: Dict + Dictionary with the differences. + """ + + # Raise TypeError if env1 or env2 are not Environment nor EnvironmentAnalysis + if not isinstance(env1, Environment) and not isinstance(env1, EnvironmentAnalysis): + raise TypeError("env1 must be an Environment or EnvironmentAnalysis object") + + if not isinstance(env2, Environment) and not isinstance(env2, EnvironmentAnalysis): + raise TypeError("env2 must be an Environment or EnvironmentAnalysis object") + + # If instances are Environment Analysis, get the equivalent environment + if isinstance(env1, EnvironmentAnalysis): + env1.process_temperature_profile_over_average_day() + env1.process_pressure_profile_over_average_day() + env1.process_wind_velocity_x_profile_over_average_day() + env1.process_wind_velocity_y_profile_over_average_day() + print() + if isinstance(env2, EnvironmentAnalysis): + env2.process_temperature_profile_over_average_day() + env2.process_pressure_profile_over_average_day() + env2.process_wind_velocity_x_profile_over_average_day() + env2.process_wind_velocity_y_profile_over_average_day() + print() + + # Plot graphs + print("\n\nAtmospheric Model Plots") + # Create height grid + grid = np.linspace(env1.elevation, env1.maxExpectedHeight) + + # Create figure + plt.figure(figsize=(9, 9)) + + # Create wind speed and wind direction subplot + ax1 = plt.subplot(221) + ax1.plot( + [env1.windSpeed(i) for i in grid], + grid, + "#ff7f0e", + label="Speed of Sound " + names[0], + ) + ax1.plot( + [env2.windSpeed(i) for i in grid], + grid, + "#ff7f0e", + label="Speed of Sound " + names[1], + ) + ax1.set_xlabel("Wind Speed (m/s)", color="#ff7f0e") + ax1.tick_params("x", colors="#ff7f0e") + # ax1up = ax1.twiny() + # ax1up.plot( + # [env1.windDirection(i) for i in grid], + # grid, + # color="#1f77b4", + # label="Density "+names[0], + # ) + # ax1up.plot( + # [env2.windDirection(i) for i in grid], + # grid, + # color="#1f77b4", + # label="Density "+names[1], + # ) + # ax1up.set_xlabel("Wind Direction (°)", color="#1f77b4") + # ax1up.tick_params("x", colors="#1f77b4") + # ax1up.set_xlim(0, 360) + ax1.set_ylabel("Height Above Sea Level (m)") + ax1.grid(True) + + # # Create density and speed of sound subplot + # ax2 = plt.subplot(222) + # ax2.plot( + # [env1.speedOfSound(i) for i in grid], + # grid, + # "#ff7f0e", + # label="Speed of Sound", + # ) + # ax2.set_xlabel("Speed of Sound (m/s)", color="#ff7f0e") + # ax2.tick_params("x", colors="#ff7f0e") + # ax2up = ax2.twiny() + # ax2up.plot([env1.density(i) for i in grid], grid, color="#1f77b4", label="Density") + # ax2up.set_xlabel("Density (kg/m³)", color="#1f77b4") + # ax2up.tick_params("x", colors="#1f77b4") + # ax2.set_ylabel("Height Above Sea Level (m)") + # ax2.grid(True) + + # Create wind u and wind v subplot + ax3 = plt.subplot(223) + ax3.plot([env1.windVelocityX(i) for i in grid], grid, label="Wind U " + names[0]) + ax3.plot([env1.windVelocityY(i) for i in grid], grid, label="Wind V " + names[0]) + ax3.plot([env2.windVelocityX(i) for i in grid], grid, label="Wind U " + names[1]) + ax3.plot([env2.windVelocityY(i) for i in grid], grid, label="Wind V " + names[1]) + ax3.legend(loc="best").set_draggable(True) + ax3.set_ylabel("Height Above Sea Level (m)") + ax3.set_xlabel("Wind Speed (m/s)") + ax3.grid(True) + + # Create pressure and temperature subplot + ax4 = plt.subplot(224) + ax4.plot([env1.pressure(i) / 100 for i in grid], grid, "#ff7f0e", label="Pressure") + ax4.set_xlabel("Pressure (hPa)", color="#ff7f0e") + ax4.tick_params("x", colors="#ff7f0e") + ax4up = ax4.twiny() + ax4up.plot( + [env1.temperature(i) for i in grid], + grid, + color="#1f77b4", + label="Temperature", + ) + ax4up.set_xlabel("Temperature (K)", color="#1f77b4") + ax4up.tick_params("x", colors="#1f77b4") + ax4.set_ylabel("Height Above Sea Level (m)") + ax4.grid(True) + + plt.subplots_adjust(wspace=0.5, hspace=0.3) + plt.show() + + return None # diff From 2ceca6a60464402b7c40f85178aedfc596125e22 Mon Sep 17 00:00:00 2001 From: patrick Date: Thu, 29 Sep 2022 21:09:21 +0200 Subject: [PATCH 3/6] Enh: Added windspeed forecast comparation plots and env configuration option --- rocketpy/EnvironmentAnalysis.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/rocketpy/EnvironmentAnalysis.py b/rocketpy/EnvironmentAnalysis.py index 6e78c4f19..595ed6049 100644 --- a/rocketpy/EnvironmentAnalysis.py +++ b/rocketpy/EnvironmentAnalysis.py @@ -72,6 +72,7 @@ def __init__( load_previous_data=None, forecast_comparaison=False, forecast_date=None, + forecast_args=None ): """Constructor for the EnvironmentAnalysis class. Parameters @@ -224,8 +225,8 @@ def __init__( longitude=self.longitude, elevation=self.elevation, ) - Env.setAtmosphericModel(type="Ensemble", file="GEFS") - Env.selectEnsembleMember(1) + forecast_args = forecast_args or {'type': 'Forecast', 'file': 'GFS'} + Env.setAtmosphericModel(**forecast_args) self.forecast[hour] = Env @@ -2549,8 +2550,8 @@ def plot_wind_profile_over_average_day(self, SAcup_altitude_constraints=False): if self.forecast: forecast = self.forecast - x = self.average_wind_profile_at_given_hour[hour][0] - y = forecast[hour].windSpeed.getValue(x) + y = self.average_wind_profile_at_given_hour[hour][1] + x = forecast[hour].windSpeed.getValue(y) * convert_units(1, "m/s", self.unit_system['wind_speed']) ax.plot(x, y, "b--") ax.label_outer() From 2002a8518d197436e9266519fa5feb0909f9e451 Mon Sep 17 00:00:00 2001 From: Lint Action Date: Thu, 29 Sep 2022 19:11:05 +0000 Subject: [PATCH 4/6] Fix code style issues with Black --- rocketpy/EnvironmentAnalysis.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/rocketpy/EnvironmentAnalysis.py b/rocketpy/EnvironmentAnalysis.py index 76cc7f794..29eaf9efc 100644 --- a/rocketpy/EnvironmentAnalysis.py +++ b/rocketpy/EnvironmentAnalysis.py @@ -231,7 +231,7 @@ def __init__( longitude=self.longitude, elevation=self.elevation, ) - forecast_args = forecast_args or {'type': 'Forecast', 'file': 'GFS'} + forecast_args = forecast_args or {"type": "Forecast", "file": "GFS"} Env.setAtmosphericModel(**forecast_args) self.forecast[hour] = Env @@ -2565,7 +2565,9 @@ def plot_wind_profile_over_average_day(self, clear_range_limits=False): if self.forecast: forecast = self.forecast y = self.average_wind_profile_at_given_hour[hour][1] - x = forecast[hour].windSpeed.getValue(y) * convert_units(1, "m/s", self.unit_system['wind_speed']) + x = forecast[hour].windSpeed.getValue(y) * convert_units( + 1, "m/s", self.unit_system["wind_speed"] + ) ax.plot(x, y, "b--") ax.label_outer() From 8b737b31cd598eb0ec08dcdad86db3d583f10d9e Mon Sep 17 00:00:00 2001 From: patrick Date: Thu, 29 Sep 2022 21:12:07 +0200 Subject: [PATCH 5/6] MAINT: removed unused func --- rocketpy/utilities.py | 133 ------------------------------------------ 1 file changed, 133 deletions(-) diff --git a/rocketpy/utilities.py b/rocketpy/utilities.py index 7bf510969..41fd251b8 100644 --- a/rocketpy/utilities.py +++ b/rocketpy/utilities.py @@ -8,7 +8,6 @@ from scipy.integrate import solve_ivp from .Environment import Environment -from .EnvironmentAnalysis import EnvironmentAnalysis from .Function import Function @@ -199,135 +198,3 @@ def du(z, u): velocityFunction() return altitudeFunction, velocityFunction, final_sol - - -def compareEnvironments(env1, env2, names=["env1", "env2"]): - """Compares two environments and plots everything. - Useful when comparing Environment Analysis results against forecasts. - - Parameters - ---------- - env1 : Environment or EnvironmentAnalysis - Environment to compare with. - env2: Environment or EnvironmentAnalysis - Environment to compare with. - - Returns - ------- - diff: Dict - Dictionary with the differences. - """ - - # Raise TypeError if env1 or env2 are not Environment nor EnvironmentAnalysis - if not isinstance(env1, Environment) and not isinstance(env1, EnvironmentAnalysis): - raise TypeError("env1 must be an Environment or EnvironmentAnalysis object") - - if not isinstance(env2, Environment) and not isinstance(env2, EnvironmentAnalysis): - raise TypeError("env2 must be an Environment or EnvironmentAnalysis object") - - # If instances are Environment Analysis, get the equivalent environment - if isinstance(env1, EnvironmentAnalysis): - env1.process_temperature_profile_over_average_day() - env1.process_pressure_profile_over_average_day() - env1.process_wind_velocity_x_profile_over_average_day() - env1.process_wind_velocity_y_profile_over_average_day() - print() - if isinstance(env2, EnvironmentAnalysis): - env2.process_temperature_profile_over_average_day() - env2.process_pressure_profile_over_average_day() - env2.process_wind_velocity_x_profile_over_average_day() - env2.process_wind_velocity_y_profile_over_average_day() - print() - - # Plot graphs - print("\n\nAtmospheric Model Plots") - # Create height grid - grid = np.linspace(env1.elevation, env1.maxExpectedHeight) - - # Create figure - plt.figure(figsize=(9, 9)) - - # Create wind speed and wind direction subplot - ax1 = plt.subplot(221) - ax1.plot( - [env1.windSpeed(i) for i in grid], - grid, - "#ff7f0e", - label="Speed of Sound " + names[0], - ) - ax1.plot( - [env2.windSpeed(i) for i in grid], - grid, - "#ff7f0e", - label="Speed of Sound " + names[1], - ) - ax1.set_xlabel("Wind Speed (m/s)", color="#ff7f0e") - ax1.tick_params("x", colors="#ff7f0e") - # ax1up = ax1.twiny() - # ax1up.plot( - # [env1.windDirection(i) for i in grid], - # grid, - # color="#1f77b4", - # label="Density "+names[0], - # ) - # ax1up.plot( - # [env2.windDirection(i) for i in grid], - # grid, - # color="#1f77b4", - # label="Density "+names[1], - # ) - # ax1up.set_xlabel("Wind Direction (°)", color="#1f77b4") - # ax1up.tick_params("x", colors="#1f77b4") - # ax1up.set_xlim(0, 360) - ax1.set_ylabel("Height Above Sea Level (m)") - ax1.grid(True) - - # # Create density and speed of sound subplot - # ax2 = plt.subplot(222) - # ax2.plot( - # [env1.speedOfSound(i) for i in grid], - # grid, - # "#ff7f0e", - # label="Speed of Sound", - # ) - # ax2.set_xlabel("Speed of Sound (m/s)", color="#ff7f0e") - # ax2.tick_params("x", colors="#ff7f0e") - # ax2up = ax2.twiny() - # ax2up.plot([env1.density(i) for i in grid], grid, color="#1f77b4", label="Density") - # ax2up.set_xlabel("Density (kg/m³)", color="#1f77b4") - # ax2up.tick_params("x", colors="#1f77b4") - # ax2.set_ylabel("Height Above Sea Level (m)") - # ax2.grid(True) - - # Create wind u and wind v subplot - ax3 = plt.subplot(223) - ax3.plot([env1.windVelocityX(i) for i in grid], grid, label="Wind U " + names[0]) - ax3.plot([env1.windVelocityY(i) for i in grid], grid, label="Wind V " + names[0]) - ax3.plot([env2.windVelocityX(i) for i in grid], grid, label="Wind U " + names[1]) - ax3.plot([env2.windVelocityY(i) for i in grid], grid, label="Wind V " + names[1]) - ax3.legend(loc="best").set_draggable(True) - ax3.set_ylabel("Height Above Sea Level (m)") - ax3.set_xlabel("Wind Speed (m/s)") - ax3.grid(True) - - # Create pressure and temperature subplot - ax4 = plt.subplot(224) - ax4.plot([env1.pressure(i) / 100 for i in grid], grid, "#ff7f0e", label="Pressure") - ax4.set_xlabel("Pressure (hPa)", color="#ff7f0e") - ax4.tick_params("x", colors="#ff7f0e") - ax4up = ax4.twiny() - ax4up.plot( - [env1.temperature(i) for i in grid], - grid, - color="#1f77b4", - label="Temperature", - ) - ax4up.set_xlabel("Temperature (K)", color="#1f77b4") - ax4up.tick_params("x", colors="#1f77b4") - ax4.set_ylabel("Height Above Sea Level (m)") - ax4.grid(True) - - plt.subplots_adjust(wspace=0.5, hspace=0.3) - plt.show() - - return None # diff From 0f070a4e940d1b310bf5e640ee7ec81ddb6ab0c4 Mon Sep 17 00:00:00 2001 From: patrick Date: Wed, 5 Oct 2022 21:52:28 +0200 Subject: [PATCH 6/6] Enh: Docs for new parameters --- rocketpy/EnvironmentAnalysis.py | 9 +++++++-- rocketpy/utilities.py | 1 - 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/rocketpy/EnvironmentAnalysis.py b/rocketpy/EnvironmentAnalysis.py index 29eaf9efc..b1f8af0ee 100644 --- a/rocketpy/EnvironmentAnalysis.py +++ b/rocketpy/EnvironmentAnalysis.py @@ -107,7 +107,6 @@ def __init__( pressureLevelDataFile=None, timezone=None, unit_system="metric", - forecast_comparison=False, forecast_date=None, forecast_args=None, maxExpectedAltitude=None, @@ -145,6 +144,12 @@ def __init__( unit_system : str, optional Unit system to be used when displaying results. Options are: SI, metric, imperial. Default is metric. + forecast_date : datetime.date, optional + Date for the forecast models. It will be requested the environment forecast + for multiple hours within that specified date. + forecast_args : dictionary, optional + Arguments for setting the forecast on the Environment class. With this argument + it is possible to change the forecast model being used. maxExpectedAltitude : float, optional Maximum expected altitude for your analysis. This is used to calculate plot limits from pressure level data profiles. If None is set, the @@ -213,7 +218,7 @@ def __init__( # Processing forecast self.forecast = None - if forecast_comparison: + if forecast_date: self.forecast = {} hours = list(self.pressureLevelDataDict.values())[0].keys() for hour in hours: diff --git a/rocketpy/utilities.py b/rocketpy/utilities.py index 41fd251b8..dc847ea3c 100644 --- a/rocketpy/utilities.py +++ b/rocketpy/utilities.py @@ -3,7 +3,6 @@ __copyright__ = "Copyright 20XX, RocketPy Team" __license__ = "MIT" -import matplotlib.pyplot as plt import numpy as np from scipy.integrate import solve_ivp