# pylint: disable=too-many-public-methods, too-many-instance-attributes
import bisect
import json
import re
import warnings
from collections import namedtuple
from datetime import datetime
import netCDF4
import numpy as np
import pytz
from rocketpy.environment.fetchers import (
fetch_atmospheric_data_from_windy,
fetch_gefs_ensemble,
fetch_gfs_file_return_dataset,
fetch_hiresw_file_return_dataset,
fetch_nam_file_return_dataset,
fetch_open_elevation,
fetch_rap_file_return_dataset,
fetch_wyoming_sounding,
)
from rocketpy.environment.tools import (
calculate_wind_heading,
calculate_wind_speed,
convert_wind_heading_to_direction,
find_latitude_index,
find_longitude_index,
find_time_index,
)
from rocketpy.environment.tools import geodesic_to_utm as geodesic_to_utm_tools
from rocketpy.environment.tools import (
get_elevation_data_from_dataset,
get_final_date_from_time_array,
get_initial_date_from_time_array,
get_interval_date_from_time_array,
get_pressure_levels_from_file,
mask_and_clean_dataset,
)
from rocketpy.environment.tools import utm_to_geodesic as utm_to_geodesic_tools
from rocketpy.environment.weather_model_mapping import WeatherModelMapping
from rocketpy.mathutils.function import NUMERICAL_TYPES, Function, funcify_method
from rocketpy.plots.environment_plots import _EnvironmentPlots
from rocketpy.prints.environment_prints import _EnvironmentPrints
from rocketpy.tools import (
bilinear_interpolation,
geopotential_height_to_geometric_height,
)
[docs]
class Environment:
"""Keeps all environment information stored, such as wind and temperature
conditions, as well as gravity.
Attributes
----------
Environment.earth_radius : float
Value of Earth's Radius as 6.3781e6 m.
Environment.air_gas_constant : float
Value of Air's Gas Constant as 287.05287 J/K/Kg
Environment.gravity : Function
Gravitational acceleration. Positive values point the
acceleration down. See :meth:`Environment.set_gravity_model` for more
information.
Environment.latitude : float
Launch site latitude.
Environment.longitude : float
Launch site longitude.
Environment.datum : string
The desired reference ellipsoid model, the following options are
available: ``SAD69``, ``WGS84``, ``NAD83``, and ``SIRGAS2000``.
Environment.initial_east : float
Launch site East UTM coordinate
Environment.initial_north : float
Launch site North UTM coordinate
Environment.initial_utm_zone : int
Launch site UTM zone number
Environment.initial_utm_letter : string
Launch site UTM letter, to keep the latitude band and describe the
UTM Zone
Environment.initial_hemisphere : string
Launch site South/North hemisphere
Environment.initial_ew : string
Launch site East/West hemisphere
Environment.elevation : float
Launch site elevation.
Environment.datetime_date : datetime
Date time of launch in UTC time zone using the ``datetime`` object.
Environment.local_date : datetime
Date time of launch in the local time zone, defined by the
``Environment.timezone`` parameter.
Environment.timezone : string
Local time zone specification. See `pytz`_. for time zone information.
.. _pytz: https://pytz.sourceforge.net/
Environment.elev_lon_array : array
Unidimensional array containing the longitude coordinates.
Environment.elev_lat_array : array
Unidimensional array containing the latitude coordinates.
Environment.elev_array : array
Two-dimensional Array containing the elevation information.
Environment.topographic_profile_activated : bool
True if the user already set a topographic profile. False otherwise.
Environment.max_expected_height : float
Maximum altitude in meters to keep weather data. The altitude must be
Above Sea Level (ASL). Especially useful for controlling plots.
Can be altered as desired by running ``max_expected_height = number``.
Environment.pressure_ISA : Function
Air pressure in Pa as a function of altitude as defined by the
International Standard Atmosphere ISO 2533.
Environment.temperature_ISA : Function
Air temperature in K as a function of altitude as defined by the
International Standard Atmosphere ISO 2533
Environment.pressure : Function
Air pressure in Pa as a function of altitude.
Environment.barometric_height : Function
Geometric height above sea level in m as a function of pressure.
Environment.temperature : Function
Air temperature in K as a function of altitude.
Environment.speed_of_sound : Function
Speed of sound in air in m/s as a function of altitude.
Environment.density : Function
Air density in kg/m³ as a function of altitude.
Environment.dynamic_viscosity : Function
Air dynamic viscosity in Pa*s as a function of altitude.
Environment.wind_speed : Function
Wind speed in m/s as a function of altitude.
Environment.wind_direction : Function
Wind direction (from which the wind blows) in degrees relative to north
(positive clockwise) as a function of altitude.
Environment.wind_heading : Function
Wind heading (direction towards which the wind blows) in degrees
relative to north (positive clockwise) as a function of altitude.
Environment.wind_velocity_x : Function
Wind U, or X (east) component of wind velocity in m/s as a function of
altitude.
Environment.wind_velocity_y : Function
Wind V, or Y (north) component of wind velocity in m/s as a function of
altitude.
Environment.atmospheric_model_type : string
Describes the atmospheric model which is being used. Can only assume the
following values: ``standard_atmosphere``, ``custom_atmosphere``,
``wyoming_sounding``, ``Forecast``, ``Reanalysis``,
``Ensemble``.
Environment.atmospheric_model_file : string
Address of the file used for the atmospheric model being used. Only
defined for ``wyoming_sounding``, ``Forecast``,
``Reanalysis``, ``Ensemble``
Environment.atmospheric_model_dict : dictionary
Dictionary used to properly interpret ``netCDF`` and ``OPeNDAP`` files.
Only defined for ``Forecast``, ``Reanalysis``, ``Ensemble``.
Environment.atmospheric_model_init_date : datetime
Datetime object instance of first available date in ``netCDF``
and ``OPeNDAP`` files when using ``Forecast``, ``Reanalysis`` or
``Ensemble``.
Environment.atmospheric_model_end_date : datetime
Datetime object instance of last available date in ``netCDF`` and
``OPeNDAP`` files when using ``Forecast``, ``Reanalysis`` or
``Ensemble``.
Environment.atmospheric_model_interval : int
Hour step between weather condition used in ``netCDF`` and
``OPeNDAP`` files when using ``Forecast``, ``Reanalysis`` or
``Ensemble``.
Environment.atmospheric_model_init_lat : float
Latitude of vertex just before the launch site in ``netCDF``
and ``OPeNDAP`` files when using ``Forecast``, ``Reanalysis`` or
``Ensemble``.
Environment.atmospheric_model_end_lat : float
Latitude of vertex just after the launch site in ``netCDF``
and ``OPeNDAP`` files when using ``Forecast``, ``Reanalysis`` or
``Ensemble``.
Environment.atmospheric_model_init_lon : float
Longitude of vertex just before the launch site in ``netCDF``
and ``OPeNDAP`` files when using ``Forecast``, ``Reanalysis`` or
``Ensemble``.
Environment.atmospheric_model_end_lon : float
Longitude of vertex just after the launch site in ``netCDF``
and ``OPeNDAP`` files when using ``Forecast``, ``Reanalysis`` or
``Ensemble``.
Environment.lat_array : array
Defined if ``netCDF`` or ``OPeNDAP`` file is used, for Forecasts,
Reanalysis and Ensembles. 2x2 matrix for each pressure level of
latitudes corresponding to the vertices of the grid cell which
surrounds the launch site.
Environment.lon_array : array
Defined if ``netCDF`` or ``OPeNDAP`` file is used, for Forecasts,
Reanalysis and Ensembles. 2x2 matrix for each pressure level of
longitudes corresponding to the vertices of the grid cell which
surrounds the launch site.
Environment.lon_index : int
Defined if ``netCDF`` or ``OPeNDAP`` file is used, for Forecasts,
Reanalysis and Ensembles. Index to a grid longitude which
is just over the launch site longitude, while ``lon_index`` - 1
points to a grid longitude which is just under the launch
site longitude.
Environment.lat_index : int
Defined if ``netCDF`` or ``OPeNDAP`` file is used, for Forecasts,
Reanalysis and Ensembles. Index to a grid latitude which
is just over the launch site latitude, while ``lat_index`` - 1
points to a grid latitude which is just under the launch
site latitude.
Environment.geopotentials : array
Defined if ``netCDF`` or ``OPeNDAP`` file is used, for Forecasts,
Reanalysis and Ensembles. 2x2 matrix for each pressure level of
geopotential heights corresponding to the vertices of the grid cell
which surrounds the launch site.
Environment.wind_us : array
Defined if ``netCDF`` or ``OPeNDAP`` file is used, for Forecasts,
Reanalysis and Ensembles. 2x2 matrix for each pressure level of
wind U (east) component corresponding to the vertices of the grid
cell which surrounds the launch site.
Environment.wind_vs : array
Defined if ``netCDF`` or ``OPeNDAP`` file is used, for Forecasts,
Reanalysis and Ensembles. 2x2 matrix for each pressure level of
wind V (north) component corresponding to the vertices of the grid
cell which surrounds the launch site.
Environment.levels : array
Defined if ``netCDF`` or ``OPeNDAP`` file is used, for Forecasts,
Reanalysis and Ensembles. List of pressure levels available in the file.
Environment.temperatures : array
Defined if ``netCDF`` or ``OPeNDAP`` file is used, for Forecasts,
Reanalysis and Ensembles. 2x2 matrix for each pressure level of
temperatures corresponding to the vertices of the grid cell which
surrounds the launch site.
Environment.time_array : array
Defined if ``netCDF`` or ``OPeNDAP`` file is used, for Forecasts,
Reanalysis and Ensembles. Array of dates available in the file.
Environment.height : array
Defined if ``netCDF`` or ``OPeNDAP`` file is used, for Forecasts,
Reanalysis and Ensembles. List of geometric height corresponding to
launch site location.
Environment.level_ensemble : array
Only defined when using Ensembles.
Environment.height_ensemble : array
Only defined when using Ensembles.
Environment.temperature_ensemble : array
Only defined when using Ensembles.
Environment.wind_u_ensemble : array
Only defined when using Ensembles.
Environment.wind_v_ensemble : array
Only defined when using Ensembles.
Environment.wind_heading_ensemble : array
Only defined when using Ensembles.
Environment.wind_direction_ensemble : array
Only defined when using Ensembles.
Environment.wind_speed_ensemble : array
Only defined when using Ensembles.
Environment.num_ensemble_members : int
Number of ensemble members. Only defined when using Ensembles.
Environment.ensemble_member : int
Current selected ensemble member. Only defined when using Ensembles.
Notes
-----
All the attributes listed as Function objects can be accessed as
regular arrays, or called as a Function. See :class:`rocketpy.Function`
for more information.
"""
[docs]
def __init__(
self,
gravity=None,
date=None,
latitude=0.0,
longitude=0.0,
elevation=0.0,
datum="SIRGAS2000",
timezone="UTC",
max_expected_height=80000.0,
):
"""Initializes the Environment class, capturing essential parameters of
the launch site, including the launch date, geographical coordinates,
and elevation. This class is designed to calculate crucial variables
for the Flight simulation, such as atmospheric air pressure, density,
and gravitational acceleration.
Note that the default atmospheric model is the International Standard
Atmosphere as defined by ISO 2533 unless specified otherwise in
:meth:`Environment.set_atmospheric_model`.
Parameters
----------
gravity : int, float, callable, string, array, optional
Surface gravitational acceleration. Positive values point the
acceleration down. If None, the Somigliana formula is used.
See :meth:`Environment.set_gravity_model` for more information.
date : list or tuple, optional
List or tuple of length 4, stating (year, month, day, hour) in the
time zone of the parameter ``timezone``.
Alternatively, can be a ``datetime`` object specifying launch
date and time. The dates are stored as follows:
- :attr:`Environment.local_date`: Local time of launch in
the time zone specified by the parameter ``timezone``.
- :attr:`Environment.datetime_date`: UTC time of launch.
Must be given if a Forecast, Reanalysis
or Ensemble, will be set as an atmospheric model.
Default is None.
See :meth:`Environment.set_date` for more information.
latitude : float, optional
Latitude in degrees (ranging from -90 to 90) of rocket
launch location. Must be given if a Forecast, Reanalysis
or Ensemble will be used as an atmospheric model or if
Open-Elevation will be used to compute elevation. Positive
values correspond to the North. Default value is 0, which
corresponds to the equator.
longitude : float, optional
Longitude in degrees (ranging from -180 to 180) of rocket
launch location. Must be given if a Forecast, Reanalysis
or Ensemble will be used as an atmospheric model or if
Open-Elevation will be used to compute elevation. Positive
values correspond to the East. Default value is 0, which
corresponds to the Greenwich Meridian.
elevation : float, optional
Elevation of launch site measured as height above sea
level in meters. Alternatively, can be set as
``Open-Elevation`` which uses the Open-Elevation API to
find elevation data. For this option, latitude and
longitude must also be specified. Default value is 0.
datum : string, optional
The desired reference ellipsoidal model, the following options are
available: "SAD69", "WGS84", "NAD83", and "SIRGAS2000". The default
is "SIRGAS2000".
timezone : string, optional
Name of the time zone. To see all time zones, import pytz and run
``print(pytz.all_timezones)``. Default time zone is "UTC".
max_expected_height : float, optional
Maximum altitude in meters to keep weather data. The altitude must
be above sea level (ASL). Especially useful for visualization. Can
be altered as desired by running ``max_expected_height = number``.
Depending on the atmospheric model, this value may be automatically
modified.
Returns
-------
None
"""
# Initialize constants and atmospheric variables
self.__initialize_empty_variables()
self.__initialize_constants()
self.__initialize_elevation_and_max_height(elevation, max_expected_height)
# Initialize plots and prints objects
self.prints = _EnvironmentPrints(self)
self.plots = _EnvironmentPlots(self)
# Set the atmosphere model to the standard atmosphere
self.set_atmospheric_model("standard_atmosphere")
# Initialize date, latitude, longitude, and Earth geometry
self.__initialize_date(date, timezone)
self.set_location(latitude, longitude)
self.__initialize_earth_geometry(datum)
self.__initialize_utm_coordinates()
# Set the gravity model
self.gravity = self.set_gravity_model(gravity)
def __initialize_constants(self):
"""Sets some important constants and atmospheric variables."""
self.earth_radius = 6.3781 * (10**6)
self.air_gas_constant = 287.05287 # in J/K/kg
self.standard_g = 9.80665
self.__weather_model_map = WeatherModelMapping()
self.__atm_type_file_to_function_map = {
"forecast": {
"GFS": fetch_gfs_file_return_dataset,
"NAM": fetch_nam_file_return_dataset,
"RAP": fetch_rap_file_return_dataset,
"HIRESW": fetch_hiresw_file_return_dataset,
},
"ensemble": {
"GEFS": fetch_gefs_ensemble,
},
}
self.__standard_atmosphere_layers = {
"geopotential_height": [ # in geopotential m
-2e3,
0,
11e3,
20e3,
32e3,
47e3,
51e3,
71e3,
80e3,
],
"temperature": [ # in K
301.15,
288.15,
216.65,
216.65,
228.65,
270.65,
270.65,
214.65,
196.65,
],
"beta": [-6.5e-3, -6.5e-3, 0, 1e-3, 2.8e-3, 0, -2.8e-3, -2e-3, 0], # in K/m
"pressure": [ # in Pa
1.27774e5,
1.01325e5,
2.26320e4,
5.47487e3,
8.680164e2,
1.10906e2,
6.69384e1,
3.95639e0,
8.86272e-2,
],
}
def __initialize_empty_variables(self):
self.atmospheric_model_file = str()
self.atmospheric_model_dict = {}
def __initialize_elevation_and_max_height(self, elevation, max_expected_height):
"""Saves the elevation and the maximum expected height."""
self.elevation = elevation
self.set_elevation(elevation)
self._max_expected_height = max_expected_height
def __initialize_date(self, date, timezone):
"""Saves the date and configure timezone."""
if date is not None:
self.set_date(date, timezone)
else:
self.date = None
self.datetime_date = None
self.local_date = None
self.timezone = None
def __initialize_earth_geometry(self, datum):
"""Initialize Earth geometry, save datum and Recalculate Earth Radius"""
self.datum = datum
self.ellipsoid = self.set_earth_geometry(datum)
self.earth_radius = self.calculate_earth_radius(
lat=self.latitude,
semi_major_axis=self.ellipsoid.semi_major_axis,
flattening=self.ellipsoid.flattening,
)
def __initialize_utm_coordinates(self):
"""Store launch site coordinates referenced to UTM projection system."""
if -80 < self.latitude < 84:
(
self.initial_east,
self.initial_north,
self.initial_utm_zone,
self.initial_utm_letter,
self.initial_hemisphere,
self.initial_ew,
) = self.geodesic_to_utm(
lat=self.latitude,
lon=self.longitude,
flattening=self.ellipsoid.flattening,
semi_major_axis=self.ellipsoid.semi_major_axis,
)
else: # pragma: no cover
warnings.warn(
"UTM coordinates are not available for latitudes "
"above 84 or below -80 degrees. The UTM conversions will fail."
)
self.initial_north = None
self.initial_east = None
self.initial_utm_zone = None
self.initial_utm_letter = None
self.initial_hemisphere = None
self.initial_ew = None
# Auxiliary private setters.
def __set_pressure_function(self, source):
self.pressure = Function(
source,
inputs="Height Above Sea Level (m)",
outputs="Pressure (Pa)",
interpolation="linear",
)
def __set_barometric_height_function(self, source):
self.barometric_height = Function(
source,
inputs="Pressure (Pa)",
outputs="Height Above Sea Level (m)",
interpolation="linear",
extrapolation="natural",
)
if callable(self.barometric_height.source):
# discretize to speed up flight simulation
self.barometric_height.set_discrete(
0,
self.max_expected_height,
100,
extrapolation="constant",
mutate_self=True,
)
def __set_temperature_function(self, source):
self.temperature = Function(
source,
inputs="Height Above Sea Level (m)",
outputs="Temperature (K)",
interpolation="linear",
)
def __set_wind_velocity_x_function(self, source):
self.wind_velocity_x = Function(
source,
inputs="Height Above Sea Level (m)",
outputs="Wind Velocity X (m/s)",
interpolation="linear",
)
def __set_wind_velocity_y_function(self, source):
self.wind_velocity_y = Function(
source,
inputs="Height Above Sea Level (m)",
outputs="Wind Velocity Y (m/s)",
interpolation="linear",
)
def __set_wind_speed_function(self, source):
self.wind_speed = Function(
source,
inputs="Height Above Sea Level (m)",
outputs="Wind Speed (m/s)",
interpolation="linear",
)
def __set_wind_direction_function(self, source):
self.wind_direction = Function(
source,
inputs="Height Above Sea Level (m)",
outputs="Wind Direction (Deg True)",
interpolation="linear",
)
def __set_wind_heading_function(self, source):
self.wind_heading = Function(
source,
inputs="Height Above Sea Level (m)",
outputs="Wind Heading (Deg True)",
interpolation="linear",
)
def __reset_barometric_height_function(self):
# NOTE: this assumes self.pressure and max_expected_height are already set.
self.barometric_height = self.pressure.inverse_function()
if callable(self.barometric_height.source):
# discretize to speed up flight simulation
self.barometric_height.set_discrete(
0,
self.max_expected_height,
100,
extrapolation="constant",
mutate_self=True,
)
self.barometric_height.set_inputs("Pressure (Pa)")
self.barometric_height.set_outputs("Height Above Sea Level (m)")
def __reset_wind_speed_function(self):
# NOTE: assume wind_velocity_x and wind_velocity_y as Function objects
self.wind_speed = (self.wind_velocity_x**2 + self.wind_velocity_y**2) ** 0.5
self.wind_speed.set_inputs("Height Above Sea Level (m)")
self.wind_speed.set_outputs("Wind Speed (m/s)")
self.wind_speed.set_title("Wind Speed Profile")
# commented because I never finished, leave it for future implementation
# def __reset_wind_heading_function(self):
# NOTE: this assumes wind_u and wind_v as numpy arrays with same length.
# TODO: should we implement arctan2 in the Function class?
# self.wind_heading = calculate_wind_heading(
# self.wind_velocity_x, self.wind_velocity_y
# )
# self.wind_heading.set_inputs("Height Above Sea Level (m)")
# self.wind_heading.set_outputs("Wind Heading (Deg True)")
# self.wind_heading.set_title("Wind Heading Profile")
def __reset_wind_direction_function(self):
self.wind_direction = convert_wind_heading_to_direction(self.wind_heading)
self.wind_direction.set_inputs("Height Above Sea Level (m)")
self.wind_direction.set_outputs("Wind Direction (Deg True)")
self.wind_direction.set_title("Wind Direction Profile")
# Validators (used to verify an attribute is being set correctly.)
def __validate_dictionary(self, file, dictionary):
# removed CMC until it is fixed.
available_models = ["GFS", "NAM", "RAP", "HIRESW", "GEFS", "ERA5"]
if isinstance(dictionary, str):
dictionary = self.__weather_model_map.get(dictionary)
elif file in available_models:
dictionary = self.__weather_model_map.get(file)
if not isinstance(dictionary, dict):
raise TypeError(
"Please specify a dictionary or choose a valid model from the "
f"following list: {available_models}"
)
return dictionary
def __validate_datetime(self):
if self.datetime_date is None:
raise ValueError(
"Please specify the launch date and time using the "
"Environment.set_date() method."
)
# Define setters
[docs]
def set_date(self, date, timezone="UTC"):
"""Set date and time of launch and update weather conditions if
date dependent atmospheric model is used.
Parameters
----------
date : list, tuple, datetime
List or tuple of length 4, stating (year, month, day, hour) in the
time zone of the parameter ``timezone``. See Notes for more
information. Alternatively, can be a ``datetime`` object specifying
launch date and time.
timezone : string, optional
Name of the time zone. To see all time zones, import pytz and run
``print(pytz.all_timezones)``. Default time zone is "UTC".
Returns
-------
None
Notes
-----
- If the ``date`` is given as a list or tuple, it should be in the same
time zone as specified by the ``timezone`` parameter. This local
time will be available in the attribute :attr:`Environment.local_date`
while the UTC time will be available in the attribute
:attr:`Environment.datetime_date`.
- If the ``date`` is given as a ``datetime`` object without a time zone,
it will be assumed to be in the same time zone as specified by the
``timezone`` parameter. However, if the ``datetime`` object has a time
zone specified in its ``tzinfo`` attribute, the ``timezone``
parameter will be ignored.
Examples
--------
Let's set the launch date as an list:
>>> date = [2000, 1, 1, 13] # January 1st, 2000 at 13:00 UTC+1
>>> env = Environment()
>>> env.set_date(date, timezone="Europe/Rome")
>>> print(env.datetime_date) # Get UTC time
2000-01-01 12:00:00+00:00
>>> print(env.local_date)
2000-01-01 13:00:00+01:00
Now let's set the launch date as a ``datetime`` object:
>>> from datetime import datetime
>>> date = datetime(2000, 1, 1, 13, 0, 0)
>>> env = Environment()
>>> env.set_date(date, timezone="Europe/Rome")
>>> print(env.datetime_date) # Get UTC time
2000-01-01 12:00:00+00:00
>>> print(env.local_date)
2000-01-01 13:00:00+01:00
"""
# Store date and configure time zone
self.timezone = timezone
tz = pytz.timezone(self.timezone)
if not isinstance(date, datetime):
local_date = datetime(*date)
else:
local_date = date
if local_date.tzinfo is None:
local_date = tz.localize(local_date)
self.date = date
self.local_date = local_date
self.datetime_date = self.local_date.astimezone(pytz.UTC)
# Update atmospheric conditions if atmosphere type is Forecast,
# Reanalysis or Ensemble
if hasattr(self, "atmospheric_model_type") and self.atmospheric_model_type in [
"Forecast",
"Reanalysis",
"Ensemble",
]:
self.set_atmospheric_model(
type=self.atmospheric_model_type,
file=self.atmospheric_model_file,
dictionary=self.atmospheric_model_dict,
)
[docs]
def set_location(self, latitude, longitude):
"""Set latitude and longitude of launch and update atmospheric
conditions if location dependent model is being used.
Parameters
----------
latitude : float
Latitude of launch site. May range from -90 to 90 degrees.
longitude : float
Longitude of launch site. Either from 0 to 360 degrees or from -180
to 180 degrees.
Returns
-------
None
"""
if not isinstance(latitude, NUMERICAL_TYPES) and isinstance(
longitude, NUMERICAL_TYPES
): # pragma: no cover
raise TypeError("Latitude and Longitude must be numbers!")
# Store latitude and longitude
self.latitude = latitude
self.longitude = longitude
# Update atmospheric conditions if atmosphere type is Forecast,
# Reanalysis or Ensemble
if hasattr(self, "atmospheric_model_type") and self.atmospheric_model_type in [
"Forecast",
"Reanalysis",
"Ensemble",
]:
self.set_atmospheric_model(
type=self.atmospheric_model_type,
file=self.atmospheric_model_file,
dictionary=self.atmospheric_model_dict,
)
[docs]
def set_gravity_model(self, gravity=None):
"""Defines the gravity model based on the given user input to the
gravity parameter. The gravity model is responsible for computing the
gravity acceleration at a given height above sea level in meters.
Parameters
----------
gravity : int, float, callable, string, list, optional
The gravitational acceleration in m/s² to be used in the
simulation, this value is positive when pointing downwards.
The input type can be one of the following:
- ``int`` or ``float``: The gravity acceleration is set as a\
constant function with respect to height;
- ``callable``: This callable should receive the height above\
sea level in meters and return the gravity acceleration;
- ``list``: The datapoints should be structured as\
``[(h_i,g_i), ...]`` where ``h_i`` is the height above sea\
level in meters and ``g_i`` is the gravity acceleration in m/s²;
- ``string``: The string should correspond to a path to a CSV file\
containing the gravity acceleration data;
- ``None``: The Somigliana formula is used to compute the gravity\
acceleration.
This parameter is used as a :class:`Function` object source, check\
out the available input types for a more detailed explanation.
Returns
-------
Function
Function object representing the gravity model.
Notes
-----
This method **does not** set the gravity acceleration, it only returns
a :class:`Function` object representing the gravity model.
Examples
--------
Let's prepare a `Environment` object with a constant gravity
acceleration:
>>> g_0 = 9.80665
>>> env_cte_g = Environment(gravity=g_0)
>>> env_cte_g.gravity([0, 100, 1000])
[np.float64(9.80665), np.float64(9.80665), np.float64(9.80665)]
It's also possible to variate the gravity acceleration by defining
its function of height:
>>> R_t = 6371000
>>> g_func = lambda h : g_0 * (R_t / (R_t + h))**2
>>> env_var_g = Environment(gravity=g_func)
>>> g = env_var_g.gravity(1000)
>>> print(f"{g:.6f}")
9.803572
"""
if gravity is None:
return self.somigliana_gravity.set_discrete(
0, self.max_expected_height, 100
)
else:
return Function(gravity, "height (m)", "gravity (m/s²)").set_discrete(
0, self.max_expected_height, 100
)
@property
def max_expected_height(self):
return self._max_expected_height
@max_expected_height.setter
def max_expected_height(self, value):
if value < self.elevation: # pragma: no cover
raise ValueError(
"Max expected height cannot be lower than the surface elevation"
)
self._max_expected_height = value
self.plots.grid = np.linspace(self.elevation, self.max_expected_height)
@funcify_method("height (m)", "gravity (m/s²)")
def somigliana_gravity(self, height):
"""Computes the gravity acceleration with the Somigliana formula [1]_.
An height correction is applied to the normal gravity that is
accurate for heights used in aviation. The formula is based on the
WGS84 ellipsoid, but is accurate for other reference ellipsoids.
Parameters
----------
height : float
Height above the reference ellipsoid in meters.
Returns
-------
Function
Function object representing the gravity model.
References
----------
.. [1] https://en.wikipedia.org/wiki/Theoretical_gravity#Somigliana_equation
"""
a = 6378137.0 # semi_major_axis
f = 1 / 298.257223563 # flattening_factor
m_rot = 3.449786506841e-3 # rotation_factor
g_e = 9.7803253359 # normal gravity at equator
k_somgl = 1.931852652458e-3 # normal gravity formula const.
first_ecc_sqrd = 6.694379990141e-3 # square of first eccentricity
# Compute quantities
sin_lat_sqrd = (np.sin(self.latitude * np.pi / 180)) ** 2
gravity_somgl = g_e * (
(1 + k_somgl * sin_lat_sqrd) / (np.sqrt(1 - first_ecc_sqrd * sin_lat_sqrd))
)
height_correction = (
1
- height * 2 / a * (1 + f + m_rot - 2 * f * sin_lat_sqrd)
+ 3 * height**2 / a**2
)
return height_correction * gravity_somgl
[docs]
def set_elevation(self, elevation="Open-Elevation"):
"""Set elevation of launch site given user input or using the
Open-Elevation API.
Parameters
----------
elevation : float, string, optional
Elevation of launch site measured as height above sea level in
meters. Alternatively, can be set as ``Open-Elevation`` which uses
the Open-Elevation API to find elevation data. For this option,
latitude and longitude must have already been specified.
See Also
--------
:meth:`rocketpy.Environment.set_location`
Returns
-------
None
"""
if elevation not in ["Open-Elevation", "SRTM"]:
# NOTE: this is assuming the elevation is a number (i.e. float, int, etc.)
self.elevation = elevation
else:
self.elevation = fetch_open_elevation(self.latitude, self.longitude)
print(f"Elevation received: {self.elevation} m")
[docs]
def set_topographic_profile( # pylint: disable=redefined-builtin, unused-argument
self, type, file, dictionary="netCDF4", crs=None
):
"""[UNDER CONSTRUCTION] Defines the Topographic profile, importing data
from previous downloaded files. Mainly data from the Shuttle Radar
Topography Mission (SRTM) and NASA Digital Elevation Model will be used
but other models and methods can be implemented in the future.
So far, this function can only handle data from NASADEM, available at:
https://cmr.earthdata.nasa.gov/search/concepts/C1546314436-LPDAAC_ECS.html
Parameters
----------
type : string
Defines the topographic model to be used, usually 'NASADEM Merged
DEM Global 1 arc second nc' can be used. To download this kind of
data, access 'https://search.earthdata.nasa.gov/search'.
NASADEM data products were derived from original telemetry data from
the Shuttle Radar Topography Mission (SRTM).
file : string
The path/name of the topographic file. Usually .nc provided by
dictionary : string, optional
Dictionary which helps to read the specified file. By default
'netCDF4' which works well with .nc files will be used.
crs : string, optional
Coordinate reference system, by default None, which will use the crs
provided by the file.
"""
if type == "NASADEM_HGT":
if dictionary == "netCDF4":
nasa_dem = netCDF4.Dataset(file, "r", format="NETCDF4")
self.elev_lon_array = nasa_dem.variables["lon"][:].tolist()
self.elev_lat_array = nasa_dem.variables["lat"][:].tolist()
self.elev_array = nasa_dem.variables["NASADEM_HGT"][:].tolist()
# crsArray = nasa_dem.variables['crs'][:].tolist().
self.topographic_profile_activated = True
print("Region covered by the Topographical file: ")
print(
f"Latitude from {self.elev_lat_array[-1]:.6f}° to "
f"{self.elev_lat_array[0]:.6f}°"
)
print(
f"Longitude from {self.elev_lon_array[0]:.6f}° to "
f"{self.elev_lon_array[-1]:.6f}°"
)
[docs]
def get_elevation_from_topographic_profile(self, lat, lon):
"""Function which receives as inputs the coordinates of a point and
finds its elevation in the provided Topographic Profile.
Parameters
----------
lat : float
latitude of the point.
lon : float
longitude of the point.
Returns
-------
elevation : float | int
Elevation provided by the topographic data, in meters.
"""
# TODO: refactor this method. pylint: disable=too-many-statements
if self.topographic_profile_activated is False: # pragma: no cover
raise ValueError(
"You must define a Topographic profile first, please use the "
"Environment.set_topographic_profile() method first."
)
# Find latitude index
# Check if reversed or sorted
if self.elev_lat_array[0] < self.elev_lat_array[-1]:
# Deal with sorted self.elev_lat_array
lat_index = bisect.bisect(self.elev_lat_array, lat)
else:
# Deal with reversed self.elev_lat_array
self.elev_lat_array.reverse()
lat_index = len(self.elev_lat_array) - bisect.bisect_left(
self.elev_lat_array, lat
)
self.elev_lat_array.reverse()
# Take care of latitude value equal to maximum longitude in the grid
if (
lat_index == len(self.elev_lat_array)
and self.elev_lat_array[lat_index - 1] == lat
):
lat_index = lat_index - 1
# Check if latitude value is inside the grid
if lat_index == 0 or lat_index == len(self.elev_lat_array):
raise ValueError(
f"Latitude {lat} not inside region covered by file, which is from "
f"{self.elev_lat_array[0]} to {self.elev_lat_array[-1]}."
)
# Find longitude index
# Determine if file uses -180 to 180 or 0 to 360
if self.elev_lon_array[0] < 0 or self.elev_lon_array[-1] < 0:
# Convert input to -180 - 180
lon = lon if lon < 180 else -180 + lon % 180
else:
# Convert input to 0 - 360
lon = lon % 360
# Check if reversed or sorted
if self.elev_lon_array[0] < self.elev_lon_array[-1]:
# Deal with sorted self.elev_lon_array
lon_index = bisect.bisect(self.elev_lon_array, lon)
else:
# Deal with reversed self.elev_lon_array
self.elev_lon_array.reverse()
lon_index = len(self.elev_lon_array) - bisect.bisect_left(
self.elev_lon_array, lon
)
self.elev_lon_array.reverse()
# Take care of longitude value equal to maximum longitude in the grid
if (
lon_index == len(self.elev_lon_array)
and self.elev_lon_array[lon_index - 1] == lon
):
lon_index = lon_index - 1
# Check if longitude value is inside the grid
if lon_index == 0 or lon_index == len(self.elev_lon_array):
raise ValueError(
f"Longitude {lon} not inside region covered by file, which is from "
f"{self.elev_lon_array[0]} to {self.elev_lon_array[-1]}."
)
# Get the elevation
elevation = self.elev_array[lat_index][lon_index]
return elevation
[docs]
def set_atmospheric_model( # pylint: disable=too-many-statements
self,
type, # pylint: disable=redefined-builtin
file=None,
dictionary=None,
pressure=None,
temperature=None,
wind_u=0,
wind_v=0,
):
"""Defines an atmospheric model for the Environment. Supported
functionality includes using data from the `International Standard
Atmosphere`, importing data from weather reanalysis, forecasts and
ensemble forecasts, importing data from upper air soundings and
inputting data as custom functions, arrays or csv files.
Parameters
----------
type : string
One of the following options:
- ``standard_atmosphere``: sets pressure and temperature profiles
corresponding to the International Standard Atmosphere defined by
ISO 2533 and ranging from -2 km to 80 km of altitude above sea
level. Note that the wind profiles are set to zero when this type
is chosen.
- ``wyoming_sounding``: sets pressure, temperature, wind-u
and wind-v profiles and surface elevation obtained from
an upper air sounding given by the file parameter through
an URL. This URL should point to a data webpage given by
selecting plot type as text: list, a station and a time at
`weather.uwyo`_.
An example of a valid link would be:
http://weather.uwyo.edu/cgi-bin/sounding?region=samer&TYPE=TEXT%3ALIST&YEAR=2019&MONTH=02&FROM=0200&TO=0200&STNM=82599
.. _weather.uwyo: http://weather.uwyo.edu/upperair/sounding.html
- ``windy_atmosphere``: sets pressure, temperature, wind-u and
wind-v profiles and surface elevation obtained from the Windy API.
See file argument to specify the model as either ``ECMWF``,
``GFS`` or ``ICON``.
- ``Forecast``: sets pressure, temperature, wind-u and wind-v
profiles and surface elevation obtained from a weather forecast
file in ``netCDF`` format or from an ``OPeNDAP`` URL, both given
through the file parameter. When this type is chosen, the date
and location of the launch should already have been set through
the date and location parameters when initializing the
Environment. The ``netCDF`` and ``OPeNDAP`` datasets must contain
at least geopotential height or geopotential, temperature, wind-u
and wind-v profiles as a function of pressure levels. If surface
geopotential or geopotential height is given, elevation is also
set. Otherwise, elevation is not changed. Profiles are
interpolated bi-linearly using supplied latitude and longitude.
The date used is the nearest one to the date supplied.
Furthermore, a dictionary must be supplied through the dictionary
parameter in order for the dataset to be accurately read. Lastly,
the dataset must use a rectangular grid sorted in either ascending
or descending order of latitude and longitude.
- ``Reanalysis``: sets pressure, temperature, wind-u and wind-v
profiles and surface elevation obtained from a weather forecast
file in ``netCDF`` format or from an ``OPeNDAP`` URL, both given
through the file parameter. When this type is chosen, the date and
location of the launch should already have been set through the
date and location parameters when initializing the Environment.
The ``netCDF`` and ``OPeNDAP`` datasets must contain at least
geopotential height or geopotential, temperature, wind-u and
wind-v profiles as a function of pressure levels. If surface
geopotential or geopotential height is given, elevation is also
set. Otherwise, elevation is not changed. Profiles are
interpolated bi-linearly using supplied latitude and longitude.
The date used is the nearest one to the date supplied.
Furthermore, a dictionary must be supplied through the dictionary
parameter in order for the dataset to be accurately read. Lastly,
the dataset must use a rectangular grid sorted in either ascending
or descending order of latitude and longitude.
- ``Ensemble``: sets pressure, temperature, wind-u and wind-v
profiles and surface elevation obtained from a weather forecast
file in ``netCDF`` format or from an ``OPeNDAP`` URL, both given
through the file parameter. When this type is chosen, the date and
location of the launch should already have been set through the
date and location parameters when initializing the Environment.
The ``netCDF`` and ``OPeNDAP`` datasets must contain at least
geopotential height or geopotential, temperature, wind-u and
wind-v profiles as a function of pressure levels. If surface
geopotential or geopotential height is given, elevation is also
set. Otherwise, elevation is not changed. Profiles are
interpolated bi-linearly using supplied latitude and longitude.
The date used is the nearest one to the date supplied.
Furthermore, a dictionary must be supplied through the dictionary
parameter in order for the dataset to be accurately read. Lastly,
the dataset must use a rectangular grid sorted in either ascending
or descending order of latitude and longitude. By default the
first ensemble forecast is activated.
.. seealso::
To activate other ensemble forecasts see
:meth:`rocketpy.Environment.select_ensemble_member`.
- ``custom_atmosphere``: sets pressure, temperature, wind-u and
wind-v profiles given though the pressure, temperature, wind-u and
wind-v parameters of this method. If pressure or temperature is
not given, it will default to the `International Standard
Atmosphere`. If the wind components are not given, it will default
to 0.
file : string, optional
String that must be given when type is either ``wyoming_sounding``,
``Forecast``, ``Reanalysis``, ``Ensemble`` or ``Windy``. It
specifies the location of the data given, either through a local
file address or a URL. If type is ``Forecast``, this parameter can
also be either ``GFS``, ``FV3``, ``RAP`` or ``NAM`` for latest of
these forecasts.
.. note::
Time reference for the Forecasts are:
- ``GFS``: `Global` - 0.25deg resolution - Updates every 6
hours, forecast for 81 points spaced by 3 hours
- ``RAP``: `Regional USA` - 0.19deg resolution - Updates hourly,
forecast for 40 points spaced hourly
- ``NAM``: `Regional CONUS Nest` - 5 km resolution - Updates
every 6 hours, forecast for 21 points spaced by 3 hours
If type is ``Ensemble``, this parameter can also be ``GEFS``
for the latest of this ensemble.
.. note::
Time referece for the Ensembles are:
- GEFS: Global, bias-corrected, 0.5deg resolution, 21 forecast
members, Updates every 6 hours, forecast for 65 points spaced
by 4 hours
- CMC (currently not available): Global, 0.5deg resolution, 21 \
forecast members, Updates every 12 hours, forecast for 65 \
points spaced by 4 hours
If type is ``Windy``, this parameter can be either ``GFS``,
``ECMWF``, ``ICON`` or ``ICONEU``. Default in this case is ``ECMWF``.
dictionary : dictionary, string, optional
Dictionary that must be given when type is either ``Forecast``,
``Reanalysis`` or ``Ensemble``. It specifies the dictionary to be
used when reading ``netCDF`` and ``OPeNDAP`` files, allowing the
correct retrieval of data. Acceptable values include ``ECMWF``,
``NOAA`` and ``UCAR`` for default dictionaries which can generally
be used to read datasets from these institutes. Alternatively, a
dictionary structure can also be given, specifying the short names
used for time, latitude, longitude, pressure levels, temperature
profile, geopotential or geopotential height profile, wind-u and
wind-v profiles in the dataset given in the file parameter.
Additionally, ensemble dictionaries must have the ensemble as well.
An example is the following dictionary, used for ``NOAA``:
.. code-block:: python
dictionary = {
"time": "time",
"latitude": "lat",
"longitude": "lon",
"level": "lev",
"ensemble": "ens",
"temperature": "tmpprs",
"surface_geopotential_height": "hgtsfc",
"geopotential_height": "hgtprs",
"geopotential": None,
"u_wind": "ugrdprs",
"v_wind": "vgrdprs",
}
pressure : float, string, array, callable, optional
This defines the atmospheric pressure profile.
Should be given if the type parameter is ``custom_atmosphere``. If not,
than the the ``Standard Atmosphere`` pressure will be used.
If a float is given, it will define a constant pressure
profile. The float should be in units of Pa.
If a string is given, it should point to a `.CSV` file
containing at most one header line and two columns of data.
The first column must be the geometric height above sea level in
meters while the second column must be the pressure in Pa.
If an array is given, it is expected to be a list or array
of coordinates (height in meters, pressure in Pa).
Finally, a callable or function is also accepted. The
function should take one argument, the height above sea
level in meters and return a corresponding pressure in Pa.
temperature : float, string, array, callable, optional
This defines the atmospheric temperature profile. Should be given
if the type parameter is ``custom_atmosphere``. If not, than the the
``Standard Atmosphere`` temperature will be used. If a float is
given, it will define a constant temperature profile. The float
should be in units of K. If a string is given, it should point to a
`.CSV` file containing at most one header line and two columns of
data. The first column must be the geometric height above sea level
in meters while the second column must be the temperature in K.
If an array is given, it is expected to be a list or array of
coordinates (height in meters, temperature in K). Finally, a
callable or function is also accepted. The function should take one
argument, the height above sea level in meters and return a
corresponding temperature in K.
wind_u : float, string, array, callable, optional
This defines the atmospheric wind-u profile, corresponding the
magnitude of the wind speed heading East. Should be given if the
type parameter is ``custom_atmosphere``. If not, it will be assumed
to be constant and equal to 0. If a float is given, it will define
a constant wind-u profile. The float should be in units of m/s. If a
string is given, it should point to a .CSV file containing at most
one header line and two columns of data. The first column must be
the geometric height above sea level in meters while the second
column must be the wind-u in m/s. If an array is given, it is
expected to be an array of coordinates (height in meters,
wind-u in m/s). Finally, a callable or function is also accepted.
The function should take one argument, the height above sea level in
meters and return a corresponding wind-u in m/s.
wind_v : float, string, array, callable, optional
This defines the atmospheric wind-v profile, corresponding the
magnitude of the wind speed heading North. Should be given if the
type parameter is ``custom_atmosphere``. If not, it will be assumed
to be constant and equal to 0. If a float is given, it will define a
constant wind-v profile. The float should be in units of m/s. If a
string is given, it should point to a .CSV file containing at most
one header line and two columns of data. The first column must be
the geometric height above sea level in meters while the second
column must be the wind-v in m/s. If an array is given, it is
expected to be an array of coordinates (height in meters, wind-v in
m/s). Finally, a callable or function is also accepted. The function
should take one argument, the height above sea level in meters and
return a corresponding wind-v in m/s.
Returns
-------
None
"""
# Save atmospheric model type
self.atmospheric_model_type = type
type = type.lower()
# Handle each case # TODO: use match case when python 3.9 is no longer supported
if type == "standard_atmosphere":
self.process_standard_atmosphere()
elif type == "wyoming_sounding":
self.process_wyoming_sounding(file)
elif type == "custom_atmosphere":
self.process_custom_atmosphere(pressure, temperature, wind_u, wind_v)
elif type == "windy":
self.process_windy_atmosphere(file)
elif type in ["forecast", "reanalysis", "ensemble"]:
dictionary = self.__validate_dictionary(file, dictionary)
try:
fetch_function = self.__atm_type_file_to_function_map[type][file]
except KeyError:
fetch_function = None
# Fetches the dataset using OpenDAP protocol or uses the file path
dataset = fetch_function() if fetch_function is not None else file
if type in ["forecast", "reanalysis"]:
self.process_forecast_reanalysis(dataset, dictionary)
else:
self.process_ensemble(dataset, dictionary)
else: # pragma: no cover
raise ValueError(f"Unknown model type '{type}'.")
if type not in ["ensemble"]:
# Ensemble already computed these values
self.calculate_density_profile()
self.calculate_speed_of_sound_profile()
self.calculate_dynamic_viscosity()
# Save dictionary and file
self.atmospheric_model_file = file
self.atmospheric_model_dict = dictionary
# Atmospheric model processing methods
[docs]
def process_standard_atmosphere(self):
"""Sets pressure and temperature profiles corresponding to the
International Standard Atmosphere defined by ISO 2533 and
ranging from -2 km to 80 km of altitude above sea level. Note
that the wind profiles are set to zero.
Returns
-------
None
"""
# Save temperature, pressure and wind profiles
self.pressure = self.pressure_ISA
self.barometric_height = self.barometric_height_ISA
self.temperature = self.temperature_ISA
# Set wind profiles to zero
self.__set_wind_direction_function(0)
self.__set_wind_heading_function(0)
self.__set_wind_velocity_x_function(0)
self.__set_wind_velocity_y_function(0)
self.__set_wind_speed_function(0)
# 80k meters is the limit of the standard atmosphere
self._max_expected_height = 80000
[docs]
def process_custom_atmosphere(
self, pressure=None, temperature=None, wind_u=0, wind_v=0
):
"""Import pressure, temperature and wind profile given by user.
Parameters
----------
pressure : float, string, array, callable, optional
This defines the atmospheric pressure profile.
Should be given if the type parameter is ``custom_atmosphere``.
If not, than the the Standard Atmosphere pressure will be used.
If a float is given, it will define a constant pressure
profile. The float should be in units of Pa.
If a string is given, it should point to a .CSV file
containing at most one header line and two columns of data.
The first column must be the geometric height above sea level in
meters while the second column must be the pressure in Pa.
If an array is given, it is expected to be a list or array
of coordinates (height in meters, pressure in Pa).
Finally, a callable or function is also accepted. The
function should take one argument, the height above sea
level in meters and return a corresponding pressure in Pa.
temperature : float, string, array, callable, optional
This defines the atmospheric temperature profile.
Should be given if the type parameter is ``custom_atmosphere``.
If not, than the the Standard Atmosphere temperature will be used.
If a float is given, it will define a constant temperature
profile. The float should be in units of K.
If a string is given, it should point to a .CSV file
containing at most one header line and two columns of data.
The first column must be the geometric height above sea level in
meters while the second column must be the temperature in K.
If an array is given, it is expected to be a list or array
of coordinates (height in meters, temperature in K).
Finally, a callable or function is also accepted. The
function should take one argument, the height above sea
level in meters and return a corresponding temperature in K.
wind_u : float, string, array, callable, optional
This defines the atmospheric wind-u profile, corresponding
the the magnitude of the wind speed heading East.
Should be given if the type parameter is ``custom_atmosphere``.
If not, it will be assumed constant and 0.
If a float is given, it will define a constant wind-u
profile. The float should be in units of m/s.
If a string is given, it should point to a .CSV file
containing at most one header line and two columns of data.
The first column must be the geometric height above sea level in
meters while the second column must be the wind-u in m/s.
If an array is given, it is expected to be an array of
coordinates (height in meters, wind-u in m/s).
Finally, a callable or function is also accepted. The
function should take one argument, the height above sea
level in meters and return a corresponding wind-u in m/s.
wind_v : float, string, array, callable, optional
This defines the atmospheric wind-v profile, corresponding
the the magnitude of the wind speed heading North.
Should be given if the type parameter is ``custom_atmosphere``.
If not, it will be assumed constant and 0.
If a float is given, it will define a constant wind-v
profile. The float should be in units of m/s.
If a string is given, it should point to a .CSV file
containing at most one header line and two columns of data.
The first column must be the geometric height above sea level in
meters while the second column must be the wind-v in m/s.
If an array is given, it is expected to be an array of
coordinates (height in meters, wind-v in m/s).
Finally, a callable or function is also accepted. The
function should take one argument, the height above sea
level in meters and return a corresponding wind-v in m/s.
Return
------
None
"""
# Initialize an estimate of the maximum expected atmospheric model height
max_expected_height = self.max_expected_height or 1000
# Save pressure profile
if pressure is None:
# Use standard atmosphere
self.pressure = self.pressure_ISA
self.barometric_height = self.barometric_height_ISA
else:
# Use custom input
self.__set_pressure_function(pressure)
self.__reset_barometric_height_function()
# Check maximum height of custom pressure input
if not callable(self.pressure.source):
max_expected_height = max(self.pressure[-1, 0], max_expected_height)
# Save temperature profile
if temperature is None:
# Use standard atmosphere
self.temperature = self.temperature_ISA
else:
self.__set_temperature_function(temperature)
# Check maximum height of custom temperature input
if not callable(self.temperature.source):
max_expected_height = max(self.temperature[-1, 0], max_expected_height)
# Save wind profile
self.__set_wind_velocity_x_function(wind_u)
self.__set_wind_velocity_y_function(wind_v)
# Check maximum height of custom wind input
if not callable(self.wind_velocity_x.source):
max_expected_height = max(self.wind_velocity_x[-1, 0], max_expected_height)
def wind_heading_func(h): # TODO: create another custom reset for heading
return calculate_wind_heading(
self.wind_velocity_x.get_value_opt(h),
self.wind_velocity_y.get_value_opt(h),
)
self.__set_wind_heading_function(wind_heading_func)
self.__reset_wind_direction_function()
self.__reset_wind_speed_function()
self._max_expected_height = max_expected_height
[docs]
def process_windy_atmosphere(
self, model="ECMWF"
): # pylint: disable=too-many-statements
"""Process data from Windy.com to retrieve atmospheric forecast data.
Parameters
----------
model : string, optional
The atmospheric model to use. Default is ``ECMWF``. Options are:
``ECMWF`` for the `ECMWF-HRES` model, ``GFS`` for the `GFS` model,
``ICON`` for the `ICON-Global` model or ``ICONEU`` for the `ICON-EU`
model.
"""
if model.lower() not in ["ecmwf", "gfs", "icon", "iconeu"]:
raise ValueError(
f"Invalid model '{model}'. "
"Valid options are 'ECMWF', 'GFS', 'ICON' or 'ICONEU'."
)
response = fetch_atmospheric_data_from_windy(
self.latitude, self.longitude, model
)
# Determine time index from model
time_array = np.array(response["data"]["hours"])
time_units = "milliseconds since 1970-01-01 00:00:00"
launch_time_in_units = netCDF4.date2num(self.datetime_date, time_units)
# Find the index of the closest time in time_array to the launch time
time_index = (np.abs(time_array - launch_time_in_units)).argmin()
# Define available pressure levels
pressure_levels = np.array(
[1000, 950, 925, 900, 850, 800, 700, 600, 500, 400, 300, 250, 200, 150]
)
# Process geopotential height array
(
geopotential_height_array,
altitude_array,
temperature_array,
wind_u_array,
wind_v_array,
) = self.__parse_windy_file(response, time_index, pressure_levels)
# Determine wind speed, heading and direction
wind_speed_array = calculate_wind_speed(wind_u_array, wind_v_array)
wind_heading_array = calculate_wind_heading(wind_u_array, wind_v_array)
wind_direction_array = convert_wind_heading_to_direction(wind_heading_array)
# Combine all data into big array
data_array = mask_and_clean_dataset(
100 * pressure_levels, # Convert hPa to Pa
altitude_array,
temperature_array,
wind_u_array,
wind_v_array,
wind_heading_array,
wind_direction_array,
wind_speed_array,
)
# Save atmospheric data
self.__set_pressure_function(data_array[:, (1, 0)])
self.__set_barometric_height_function(data_array[:, (0, 1)])
self.__set_temperature_function(data_array[:, (1, 2)])
self.__set_wind_velocity_x_function(data_array[:, (1, 3)])
self.__set_wind_velocity_y_function(data_array[:, (1, 4)])
self.__set_wind_heading_function(data_array[:, (1, 5)])
self.__set_wind_direction_function(data_array[:, (1, 6)])
self.__set_wind_speed_function(data_array[:, (1, 7)])
# Save maximum expected height
self._max_expected_height = max(altitude_array[0], altitude_array[-1])
# Get elevation data from file
self.elevation = float(response["header"]["elevation"])
# Compute info data
self.atmospheric_model_init_date = get_initial_date_from_time_array(
time_array, time_units
)
self.atmospheric_model_end_date = get_final_date_from_time_array(
time_array, time_units
)
self.atmospheric_model_interval = get_interval_date_from_time_array(
time_array, time_units
)
self.atmospheric_model_init_lat = self.latitude
self.atmospheric_model_end_lat = self.latitude
self.atmospheric_model_init_lon = self.longitude
self.atmospheric_model_end_lon = self.longitude
# Save debugging data
self.geopotentials = geopotential_height_array
self.wind_us = wind_u_array
self.wind_vs = wind_v_array
self.levels = pressure_levels
self.temperatures = temperature_array
self.time_array = time_array
self.height = altitude_array
def __parse_windy_file(self, response, time_index, pressure_levels):
geopotential_height_array = np.array(
[response["data"][f"gh-{pL}h"][time_index] for pL in pressure_levels]
)
# Convert geopotential height to geometric altitude (ASL)
altitude_array = geopotential_height_to_geometric_height(
geopotential_height_array, self.earth_radius
)
# Process temperature array (in Kelvin)
temperature_array = np.array(
[response["data"][f"temp-{pL}h"][time_index] for pL in pressure_levels]
)
# Process wind-u and wind-v array (in m/s)
wind_u_array = np.array(
[response["data"][f"wind_u-{pL}h"][time_index] for pL in pressure_levels]
)
wind_v_array = np.array(
[response["data"][f"wind_v-{pL}h"][time_index] for pL in pressure_levels]
)
return (
geopotential_height_array,
altitude_array,
temperature_array,
wind_u_array,
wind_v_array,
)
[docs]
def process_wyoming_sounding(self, file): # pylint: disable=too-many-statements
"""Import and process the upper air sounding data from `Wyoming
Upper Air Soundings` database given by the url in file. Sets
pressure, temperature, wind-u, wind-v profiles and surface elevation.
Parameters
----------
file : string
URL of an upper air sounding data output from `Wyoming
Upper Air Soundings` database.
Example:
http://weather.uwyo.edu/cgi-bin/sounding?region=samer&TYPE=TEXT%3ALIST&YEAR=2019&MONTH=02&FROM=0200&TO=0200&STNM=82599
Notes
-----
More can be found at: http://weather.uwyo.edu/upperair/sounding.html.
Returns
-------
None
"""
# Request Wyoming Sounding from file url
response = fetch_wyoming_sounding(file)
# Process Wyoming Sounding by finding data table and station info
response_split_text = re.split("(<.{0,1}PRE>)", response.text)
data_table = response_split_text[2]
station_info = response_split_text[6]
# Transform data table into np array
data_array = []
for line in data_table.split("\n")[5:-1]:
# Split data table into lines and remove header and footer
columns = re.split(" +", line) # Split line into columns
# 12 is the number of column entries when all entries are given
if len(columns) == 12:
data_array.append(columns[1:])
data_array = np.array(data_array, dtype=float)
# Retrieve pressure from data array
data_array[:, 0] = 100 * data_array[:, 0] # Converts hPa to Pa
self.__set_pressure_function(data_array[:, (1, 0)])
self.__set_barometric_height_function(data_array[:, (0, 1)])
# Retrieve temperature from data array
data_array[:, 2] = data_array[:, 2] + 273.15 # Converts C to K
self.__set_temperature_function(data_array[:, (1, 2)])
# Retrieve wind-u and wind-v from data array
## Converts Knots to m/s
data_array[:, 7] = data_array[:, 7] * 1.852 / 3.6
## Convert wind direction to wind heading
data_array[:, 5] = (data_array[:, 6] + 180) % 360
data_array[:, 3] = data_array[:, 7] * np.sin(data_array[:, 5] * np.pi / 180)
data_array[:, 4] = data_array[:, 7] * np.cos(data_array[:, 5] * np.pi / 180)
# Convert geopotential height to geometric height
data_array[:, 1] = geopotential_height_to_geometric_height(
data_array[:, 1], self.earth_radius
)
# Save atmospheric data
self.__set_wind_velocity_x_function(data_array[:, (1, 3)])
self.__set_wind_velocity_y_function(data_array[:, (1, 4)])
self.__set_wind_heading_function(data_array[:, (1, 5)])
self.__set_wind_direction_function(data_array[:, (1, 6)])
self.__set_wind_speed_function(data_array[:, (1, 7)])
# Retrieve station elevation from station info
station_elevation_text = station_info.split("\n")[6]
# Convert station elevation text into float value
self.elevation = float(
re.findall(r"[0-9]+\.[0-9]+|[0-9]+", station_elevation_text)[0]
)
# Save maximum expected height
self._max_expected_height = data_array[-1, 1]
[docs]
def process_noaaruc_sounding(self, file): # pylint: disable=too-many-statements
"""Import and process the upper air sounding data from `NOAA
Ruc Soundings` database (https://rucsoundings.noaa.gov/) given as
ASCII GSD format pages passed by its url to the file parameter. Sets
pressure, temperature, wind-u, wind-v profiles and surface elevation.
Parameters
----------
file : string
URL of an upper air sounding data output from `NOAA Ruc Soundings`
in ASCII GSD format.
Example:
https://rucsoundings.noaa.gov/get_raobs.cgi?data_source=RAOB&latest=latest&start_year=2019&start_month_name=Feb&start_mday=5&start_hour=12&start_min=0&n_hrs=1.0&fcst_len=shortest&airport=83779&text=Ascii%20text%20%28GSD%20format%29&hydrometeors=false&start=latest
See also
--------
This method is deprecated and will be fully deleted in version 1.8.0.
Returns
-------
None
"""
warnings.warn(
"NOAA RUC models are no longer available. "
"This method is deprecated and will be fully deleted in version 1.8.0.",
DeprecationWarning,
)
return file
[docs]
def process_forecast_reanalysis(
self, file, dictionary
): # pylint: disable=too-many-locals,too-many-statements
"""Import and process atmospheric data from weather forecasts
and reanalysis given as ``netCDF`` or ``OPeNDAP`` files.
Sets pressure, temperature, wind-u and wind-v
profiles and surface elevation obtained from a weather
file in ``netCDF`` format or from an ``OPeNDAP`` URL, both
given through the file parameter. The date and location of the launch
should already have been set through the date and
location parameters when initializing the Environment.
The ``netCDF`` and ``OPeNDAP`` datasets must contain at least
geopotential height or geopotential, temperature,
wind-u and wind-v profiles as a function of pressure levels.
If surface geopotential or geopotential height is given,
elevation is also set. Otherwise, elevation is not changed.
Profiles are interpolated bi-linearly using supplied
latitude and longitude. The date used is the nearest one
to the date supplied. Furthermore, a dictionary must be
supplied through the dictionary parameter in order for the
dataset to be accurately read. Lastly, the dataset must use
a rectangular grid sorted in either ascending or descending
order of latitude and longitude.
Parameters
----------
file : string
String containing path to local ``netCDF`` file or URL of an
``OPeNDAP`` file, such as NOAA's NOMAD or UCAR TRHEDDS server.
dictionary : dictionary
Specifies the dictionary to be used when reading ``netCDF`` and
``OPeNDAP`` files, allowing for the correct retrieval of data.
The dictionary structure should specify the short names
used for time, latitude, longitude, pressure levels,
temperature profile, geopotential or geopotential height
profile, wind-u and wind-v profiles in the dataset given in
the file parameter. An example is the following dictionary,
generally used to read ``OPeNDAP`` files from NOAA's NOMAD
server:
.. code-block:: python
dictionary = {
"time": "time",
"latitude": "lat",
"longitude": "lon",
"level": "lev",
"temperature": "tmpprs",
"surface_geopotential_height": "hgtsfc",
"geopotential_height": "hgtprs",
"geopotential": None,
"u_wind": "ugrdprs",
"v_wind": "vgrdprs",
}
Returns
-------
None
"""
# Check if date, lat and lon are known
self.__validate_datetime()
# Read weather file
if isinstance(file, str):
data = netCDF4.Dataset(file)
if dictionary["time"] not in data.variables.keys():
dictionary = self.__weather_model_map.get("ECMWF_v0")
else:
data = file
# Get time, latitude and longitude data from file
time_array = data.variables[dictionary["time"]]
lon_list = data.variables[dictionary["longitude"]][:].tolist()
lat_list = data.variables[dictionary["latitude"]][:].tolist()
# Find time, latitude and longitude indexes
time_index = find_time_index(self.datetime_date, time_array)
lon, lon_index = find_longitude_index(self.longitude, lon_list)
_, lat_index = find_latitude_index(self.latitude, lat_list)
# Get pressure level data from file
levels = get_pressure_levels_from_file(data, dictionary)
# Get geopotential data from file
try:
geopotentials = data.variables[dictionary["geopotential_height"]][
time_index, :, (lat_index - 1, lat_index), (lon_index - 1, lon_index)
]
except KeyError:
try:
geopotentials = (
data.variables[dictionary["geopotential"]][
time_index,
:,
(lat_index - 1, lat_index),
(lon_index - 1, lon_index),
]
/ self.standard_g
)
except KeyError as e:
raise ValueError(
"Unable to read geopotential height"
" nor geopotential from file. At least"
" one of them is necessary. Check "
" file and dictionary."
) from e
# Get temperature from file
try:
temperatures = data.variables[dictionary["temperature"]][
time_index, :, (lat_index - 1, lat_index), (lon_index - 1, lon_index)
]
except Exception as e:
raise ValueError(
"Unable to read temperature from file. Check file and dictionary."
) from e
# Get wind data from file
try:
wind_us = data.variables[dictionary["u_wind"]][
time_index, :, (lat_index - 1, lat_index), (lon_index - 1, lon_index)
]
except KeyError as e:
raise ValueError(
"Unable to read wind-u component. Check file and dictionary."
) from e
try:
wind_vs = data.variables[dictionary["v_wind"]][
time_index, :, (lat_index - 1, lat_index), (lon_index - 1, lon_index)
]
except KeyError as e:
raise ValueError(
"Unable to read wind-v component. Check file and dictionary."
) from e
# Prepare for bilinear interpolation
x, y = self.latitude, lon
x1, y1 = lat_list[lat_index - 1], lon_list[lon_index - 1]
x2, y2 = lat_list[lat_index], lon_list[lon_index]
# Determine properties in lat, lon
height = bilinear_interpolation(
x,
y,
x1,
x2,
y1,
y2,
geopotentials[:, 0, 0],
geopotentials[:, 0, 1],
geopotentials[:, 1, 0],
geopotentials[:, 1, 1],
)
temper = bilinear_interpolation(
x,
y,
x1,
x2,
y1,
y2,
temperatures[:, 0, 0],
temperatures[:, 0, 1],
temperatures[:, 1, 0],
temperatures[:, 1, 1],
)
wind_u = bilinear_interpolation(
x,
y,
x1,
x2,
y1,
y2,
wind_us[:, 0, 0],
wind_us[:, 0, 1],
wind_us[:, 1, 0],
wind_us[:, 1, 1],
)
wind_v = bilinear_interpolation(
x,
y,
x1,
x2,
y1,
y2,
wind_vs[:, 0, 0],
wind_vs[:, 0, 1],
wind_vs[:, 1, 0],
wind_vs[:, 1, 1],
)
# Determine wind speed, heading and direction
wind_speed = calculate_wind_speed(wind_u, wind_v)
wind_heading = calculate_wind_heading(wind_u, wind_v)
wind_direction = convert_wind_heading_to_direction(wind_heading)
# Convert geopotential height to geometric height
height = geopotential_height_to_geometric_height(height, self.earth_radius)
# Combine all data into big array
data_array = mask_and_clean_dataset(
levels,
height,
temper,
wind_u,
wind_v,
wind_heading,
wind_direction,
wind_speed,
)
# Save atmospheric data
self.__set_pressure_function(data_array[:, (1, 0)])
self.__set_barometric_height_function(data_array[:, (0, 1)])
self.__set_temperature_function(data_array[:, (1, 2)])
self.__set_wind_velocity_x_function(data_array[:, (1, 3)])
self.__set_wind_velocity_y_function(data_array[:, (1, 4)])
self.__set_wind_heading_function(data_array[:, (1, 5)])
self.__set_wind_direction_function(data_array[:, (1, 6)])
self.__set_wind_speed_function(data_array[:, (1, 7)])
# Save maximum expected height
self._max_expected_height = max(height[0], height[-1])
# Get elevation data from file
if dictionary["surface_geopotential_height"] is not None:
self.elevation = get_elevation_data_from_dataset(
dictionary, data, time_index, lat_index, lon_index, x, y, x1, x2, y1, y2
)
# Compute info data
self.atmospheric_model_init_date = get_initial_date_from_time_array(time_array)
self.atmospheric_model_end_date = get_final_date_from_time_array(time_array)
if self.atmospheric_model_init_date != self.atmospheric_model_end_date:
self.atmospheric_model_interval = get_interval_date_from_time_array(
time_array
)
else:
self.atmospheric_model_interval = 0
self.atmospheric_model_init_lat = lat_list[0]
self.atmospheric_model_end_lat = lat_list[-1]
self.atmospheric_model_init_lon = lon_list[0]
self.atmospheric_model_end_lon = lon_list[-1]
# Save debugging data
self.lat_array = lat_list
self.lon_array = lon_list
self.lon_index = lon_index
self.lat_index = lat_index
self.geopotentials = geopotentials
self.wind_us = wind_us
self.wind_vs = wind_vs
self.levels = levels
self.temperatures = temperatures
self.time_array = time_array[:].tolist()
self.height = height
# Close weather data
data.close()
[docs]
def process_ensemble(
self, file, dictionary
): # pylint: disable=too-many-locals,too-many-statements
"""Import and process atmospheric data from weather ensembles
given as ``netCDF`` or ``OPeNDAP`` files. Sets pressure, temperature,
wind-u and wind-v profiles and surface elevation obtained from a weather
ensemble file in ``netCDF`` format or from an ``OPeNDAP`` URL, both
given through the file parameter. The date and location of the launch
should already have been set through the date and location parameters
when initializing the Environment. The ``netCDF`` and ``OPeNDAP``
datasets must contain at least geopotential height or geopotential,
temperature, wind-u and wind-v profiles as a function of pressure
levels. If surface geopotential or geopotential height is given,
elevation is also set. Otherwise, elevation is not changed. Profiles are
interpolated bi-linearly using supplied latitude and longitude. The date
used is the nearest one to the date supplied. Furthermore, a dictionary
must be supplied through the dictionary parameter in order for the
dataset to be accurately read. Lastly, the dataset must use a
rectangular grid sorted in either ascending or descending order of
latitude and longitude. By default the first ensemble forecast is
activated. To activate other ensemble forecasts see
:meth:`Environment.select_ensemble_member()`.
Parameters
----------
file : string
String containing path to local ``.nc`` file or URL of an
``OPeNDAP`` file, such as NOAA's NOMAD or UCAR TRHEDDS server.
dictionary : dictionary
Specifies the dictionary to be used when reading ``netCDF`` and
``OPeNDAP`` files, allowing for the correct retrieval of data.
The dictionary structure should specify the short names
used for time, latitude, longitude, pressure levels,
temperature profile, geopotential or geopotential height
profile, wind-u and wind-v profiles in the dataset given in
the file parameter. An example is the following dictionary,
generally used to read ``OPeNDAP`` files from NOAA's NOMAD
server:
.. code-block:: python
dictionary = {
"time": "time",
"latitude": "lat",
"longitude": "lon",
"level": "lev",
"ensemble": "ens",
"surface_geopotential_height": "hgtsfc",
"geopotential_height": "hgtprs",
"geopotential": None,
"u_wind": "ugrdprs",
"v_wind": "vgrdprs",
}
See also
--------
See the :class:``rocketpy.environment.weather_model_mapping`` for some
dictionary examples.
"""
# Check if date, lat and lon are known
self.__validate_datetime()
# Read weather file
if isinstance(file, str):
data = netCDF4.Dataset(file)
else:
data = file
# Get time, latitude and longitude data from file
time_array = data.variables[dictionary["time"]]
lon_list = data.variables[dictionary["longitude"]][:].tolist()
lat_list = data.variables[dictionary["latitude"]][:].tolist()
# Find time, latitude and longitude indexes
time_index = find_time_index(self.datetime_date, time_array)
lon, lon_index = find_longitude_index(self.longitude, lon_list)
_, lat_index = find_latitude_index(self.latitude, lat_list)
# Get ensemble data from file
try:
num_members = len(data.variables[dictionary["ensemble"]][:])
except KeyError as e:
raise ValueError(
"Unable to read ensemble data from file. Check file and dictionary."
) from e
# Get pressure level data from file
levels = get_pressure_levels_from_file(data, dictionary)
inverse_dictionary = {v: k for k, v in dictionary.items()}
param_dictionary = {
"time": time_index,
"ensemble": range(num_members),
"level": range(len(levels)),
"latitude": (lat_index - 1, lat_index),
"longitude": (lon_index - 1, lon_index),
}
# Get dimensions
try:
dimensions = data.variables[dictionary["geopotential_height"]].dimensions[:]
except KeyError:
dimensions = data.variables[dictionary["geopotential"]].dimensions[:]
# Get params
params = tuple(param_dictionary[inverse_dictionary[dim]] for dim in dimensions)
# Get geopotential data from file
try:
geopotentials = data.variables[dictionary["geopotential_height"]][params]
except KeyError:
try:
geopotentials = (
data.variables[dictionary["geopotential"]][params] / self.standard_g
)
except KeyError as e:
raise ValueError(
"Unable to read geopotential height nor geopotential from file. "
"At least one of them is necessary. Check file and dictionary."
) from e
# Get temperature from file
try:
temperatures = data.variables[dictionary["temperature"]][params]
except KeyError as e:
raise ValueError(
"Unable to read temperature from file. Check file and dictionary."
) from e
# Get wind data from file
try:
wind_us = data.variables[dictionary["u_wind"]][params]
except KeyError as e:
raise ValueError(
"Unable to read wind-u component. Check file and dictionary."
) from e
try:
wind_vs = data.variables[dictionary["v_wind"]][params]
except KeyError as e:
raise ValueError(
"Unable to read wind-v component. Check file and dictionary."
) from e
# Prepare for bilinear interpolation
x, y = self.latitude, lon
x1, y1 = lat_list[lat_index - 1], lon_list[lon_index - 1]
x2, y2 = lat_list[lat_index], lon_list[lon_index]
# Determine properties in lat, lon
height = bilinear_interpolation(
x,
y,
x1,
x2,
y1,
y2,
geopotentials[:, :, 0, 0],
geopotentials[:, :, 0, 1],
geopotentials[:, :, 1, 0],
geopotentials[:, :, 1, 1],
)
temper = bilinear_interpolation(
x,
y,
x1,
x2,
y1,
y2,
temperatures[:, :, 0, 0],
temperatures[:, :, 0, 1],
temperatures[:, :, 1, 0],
temperatures[:, :, 1, 1],
)
wind_u = bilinear_interpolation(
x,
y,
x1,
x2,
y1,
y2,
wind_us[:, :, 0, 0],
wind_us[:, :, 0, 1],
wind_us[:, :, 1, 0],
wind_us[:, :, 1, 1],
)
wind_v = bilinear_interpolation(
x,
y,
x1,
x2,
y1,
y2,
wind_vs[:, :, 0, 0],
wind_vs[:, :, 0, 1],
wind_vs[:, :, 1, 0],
wind_vs[:, :, 1, 1],
)
# Determine wind speed, heading and direction
wind_speed = calculate_wind_speed(wind_u, wind_v)
wind_heading = calculate_wind_heading(wind_u, wind_v)
wind_direction = convert_wind_heading_to_direction(wind_heading)
# Convert geopotential height to geometric height
height = geopotential_height_to_geometric_height(height, self.earth_radius)
# Save ensemble data
self.level_ensemble = levels
self.height_ensemble = height
self.temperature_ensemble = temper
self.wind_u_ensemble = wind_u
self.wind_v_ensemble = wind_v
self.wind_heading_ensemble = wind_heading
self.wind_direction_ensemble = wind_direction
self.wind_speed_ensemble = wind_speed
self.num_ensemble_members = num_members
# Activate default ensemble
self.select_ensemble_member()
# Get elevation data from file
if dictionary["surface_geopotential_height"] is not None:
self.elevation = get_elevation_data_from_dataset(
dictionary, data, time_index, lat_index, lon_index, x, y, x1, x2, y1, y2
)
# Compute info data
self.atmospheric_model_init_date = get_initial_date_from_time_array(time_array)
self.atmospheric_model_end_date = get_final_date_from_time_array(time_array)
self.atmospheric_model_interval = get_interval_date_from_time_array(time_array)
self.atmospheric_model_init_lat = lat_list[0]
self.atmospheric_model_end_lat = lat_list[-1]
self.atmospheric_model_init_lon = lon_list[0]
self.atmospheric_model_end_lon = lon_list[-1]
# Save debugging data
self.lat_array = lat_list
self.lon_array = lon_list
self.lon_index = lon_index
self.lat_index = lat_index
self.geopotentials = geopotentials
self.wind_us = wind_us
self.wind_vs = wind_vs
self.levels = levels
self.temperatures = temperatures
self.time_array = time_array[:].tolist()
self.height = height
# Close weather data
data.close()
[docs]
def select_ensemble_member(self, member=0):
"""Activates the specified ensemble member, ensuring all atmospheric
variables read from the Environment instance correspond to the selected
ensemble member.
Parameters
----------
member : int, optional
The ensemble member to activate. Index starts from 0. Default is 0.
Raises
------
ValueError
If the specified ensemble member index is out of range.
Notes
-----
The first ensemble member (index 0) is activated by default when loading
an ensemble model. This member typically represents a control member
that is generated without perturbations. Other ensemble members are
generated by perturbing the control member.
"""
# Verify ensemble member
if member >= self.num_ensemble_members:
raise ValueError(
f"Please choose member from 0 to {self.num_ensemble_members - 1}"
)
# Read ensemble member
levels = self.level_ensemble[:]
height = self.height_ensemble[member, :]
temperature = self.temperature_ensemble[member, :]
wind_u = self.wind_u_ensemble[member, :]
wind_v = self.wind_v_ensemble[member, :]
wind_heading = self.wind_heading_ensemble[member, :]
wind_direction = self.wind_direction_ensemble[member, :]
wind_speed = self.wind_speed_ensemble[member, :]
# Combine all data into big array
data_array = mask_and_clean_dataset(
levels,
height,
temperature,
wind_u,
wind_v,
wind_heading,
wind_direction,
wind_speed,
)
# Save atmospheric data
self.__set_pressure_function(data_array[:, (1, 0)])
self.__set_barometric_height_function(data_array[:, (0, 1)])
self.__set_temperature_function(data_array[:, (1, 2)])
self.__set_wind_velocity_x_function(data_array[:, (1, 3)])
self.__set_wind_velocity_y_function(data_array[:, (1, 4)])
self.__set_wind_heading_function(data_array[:, (1, 5)])
self.__set_wind_direction_function(data_array[:, (1, 6)])
self.__set_wind_speed_function(data_array[:, (1, 7)])
# Save other attributes
self._max_expected_height = max(height[0], height[-1])
self.ensemble_member = member
# Update air density, speed of sound and dynamic viscosity
self.calculate_density_profile()
self.calculate_speed_of_sound_profile()
self.calculate_dynamic_viscosity()
[docs]
def load_international_standard_atmosphere(self): # pragma: no cover
"""Defines the pressure and temperature profile functions set
by `ISO 2533` for the International Standard atmosphere and saves
them as ``Environment.pressure_ISA`` and ``Environment.temperature_ISA``.
Notes
-----
This method is **deprecated** and will be removed in version 1.6.0. You
can access :meth:`rocketpy.Environment.pressure_ISA` and
:meth:`rocketpy.Environment.temperature_ISA` directly without the need
to call this method.
"""
warnings.warn(
"load_international_standard_atmosphere() is deprecated in version "
"1.5.0 and will be removed in version 1.7.0. This method is no longer "
"needed as the International Standard Atmosphere is already calculated "
"when the Environment object is created.",
DeprecationWarning,
)
@funcify_method("Height Above Sea Level (m)", "Pressure (Pa)", "spline", "natural")
def pressure_ISA(self):
"""Pressure, in Pa, as a function of height above sea level as defined
by the `International Standard Atmosphere ISO 2533`."""
# Retrieve lists
pressure = self.__standard_atmosphere_layers["pressure"]
geopotential_height = self.__standard_atmosphere_layers["geopotential_height"]
temperature = self.__standard_atmosphere_layers["temperature"]
beta = self.__standard_atmosphere_layers["beta"]
# Get constants
earth_radius = self.earth_radius
g = self.standard_g
R = self.air_gas_constant
# Create function to compute pressure at a given geometric height
def pressure_function(h):
"""Computes the pressure at a given geometric height h using the
International Standard Atmosphere model."""
# Convert geometric to geopotential height
H = earth_radius * h / (earth_radius + h)
# Check if height is within bounds, return extrapolated value if not
if H < -2000:
return pressure[0]
elif H > 80000:
return pressure[-1]
# Find layer that contains height h
layer = bisect.bisect(geopotential_height, H) - 1
# Retrieve layer base geopotential height, temp, beta and pressure
base_geopotential_height = geopotential_height[layer]
base_temperature = temperature[layer]
base_pressure = pressure[layer]
B = beta[layer]
# Compute pressure
if B != 0:
P = base_pressure * (
1 + (B / base_temperature) * (H - base_geopotential_height)
) ** (-g / (B * R))
else:
T = base_temperature + B * (H - base_geopotential_height)
P = base_pressure * np.exp(
-(H - base_geopotential_height) * (g / (R * T))
)
return P
# Discretize this Function to speed up the trajectory simulation
altitudes = np.linspace(0, 80000, 100) # TODO: should be -2k instead of 0
pressures = [pressure_function(h) for h in altitudes]
return np.column_stack([altitudes, pressures])
@funcify_method("Pressure (Pa)", "Height Above Sea Level (m)")
def barometric_height_ISA(self):
"""Returns the inverse function of the ``pressure_ISA`` function."""
return self.pressure_ISA.inverse_function()
@funcify_method("Height Above Sea Level (m)", "Temperature (K)", "linear")
def temperature_ISA(self):
"""Air temperature, in K, as a function of altitude as defined by the
`International Standard Atmosphere ISO 2533`."""
temperature = self.__standard_atmosphere_layers["temperature"]
geopotential_height = self.__standard_atmosphere_layers["geopotential_height"]
altitude_asl = [
geopotential_height_to_geometric_height(h, self.earth_radius)
for h in geopotential_height
]
return np.column_stack([altitude_asl, temperature])
[docs]
def calculate_density_profile(self):
r"""Compute the density of the atmosphere as a function of
height. This function is automatically called whenever a new atmospheric
model is set.
Notes
-----
1. The density is calculated as:
.. math:: \rho = \frac{P}{RT}
Examples
--------
Creating an Environment object and calculating the density
at Sea Level:
>>> env = Environment()
>>> env.calculate_density_profile()
>>> float(env.density(0))
1.225000018124288
Creating an Environment object and calculating the density
at 1000m above Sea Level:
>>> env = Environment()
>>> env.calculate_density_profile()
>>> float(env.density(1000))
1.1115112430077818
"""
# Retrieve pressure P, gas constant R and temperature T
P = self.pressure
R = self.air_gas_constant
T = self.temperature
# Compute density using P/RT
D = P / (R * T)
# Set new output for the calculated density
D.set_outputs("Air Density (kg/m³)")
# Save calculated density
self.density = D
[docs]
def calculate_speed_of_sound_profile(self):
r"""Compute the speed of sound in the atmosphere as a function
of height. This function is automatically called whenever a new
atmospheric model is set.
Notes
-----
1. The speed of sound is calculated as:
.. math:: a = \sqrt{\gamma \cdot R \cdot T}
"""
# Retrieve gas constant R and temperature T
R = self.air_gas_constant
T = self.temperature
G = 1.4
# Compute speed of sound using sqrt(gamma*R*T)
a = (G * R * T) ** 0.5
# Set new output for the calculated speed of sound
a.set_outputs("Speed of Sound (m/s)")
# Save calculated speed of sound
self.speed_of_sound = a
[docs]
def calculate_dynamic_viscosity(self):
r"""Compute the dynamic viscosity of the atmosphere as a function of
height by using the formula given in ISO 2533. This function is
automatically called whenever a new atmospheric model is set.
Notes
-----
1. The dynamic viscosity is calculated as:
.. math::
\mu = \frac{B \cdot T^{1.5}}{(T + S)}
where `B` and `S` are constants, and `T` is the temperature.
2. This equation is invalid for very high or very low temperatures.
3. Also invalid under conditions occurring at altitudes above 90 km.
"""
# Retrieve temperature T and set constants
T = self.temperature
B = 1.458e-6 # Kg/m/s/K^0.5
S = 110.4 # K
# Compute dynamic viscosity using u = B*T^(1.4)/(T+S) (See ISO2533)
u = (B * T ** (1.5)) / (T + S)
# Set new output for the calculated density
u.set_outputs("Dynamic Viscosity (Pa s)")
# Save calculated density
self.dynamic_viscosity = u
[docs]
def add_wind_gust(self, wind_gust_x, wind_gust_y):
"""Adds a function to the current stored wind profile, in order to
simulate a wind gust.
Parameters
----------
wind_gust_x : float, callable
Callable, function of altitude, which will be added to the
x velocity of the current stored wind profile. If float is given,
it will be considered as a constant function in altitude.
wind_gust_y : float, callable
Callable, function of altitude, which will be added to the
y velocity of the current stored wind profile. If float is given,
it will be considered as a constant function in altitude.
"""
# Recalculate wind_velocity_x and wind_velocity_y
self.__set_wind_velocity_x_function(self.wind_velocity_x + wind_gust_x)
self.__set_wind_velocity_y_function(self.wind_velocity_y + wind_gust_y)
# Reset wind heading and velocity magnitude
self.wind_heading = Function(
lambda h: (180 / np.pi)
* np.arctan2(
self.wind_velocity_x.get_value_opt(h),
self.wind_velocity_y.get_value_opt(h),
)
% 360,
"Height (m)",
"Wind Heading (degrees)",
extrapolation="constant",
)
self.wind_speed = Function(
lambda h: (
self.wind_velocity_x.get_value_opt(h) ** 2
+ self.wind_velocity_y.get_value_opt(h) ** 2
)
** 0.5,
"Height (m)",
"Wind Speed (m/s)",
extrapolation="constant",
)
[docs]
def info(self):
"""Prints important data and graphs available about the Environment."""
self.prints.all()
self.plots.info()
[docs]
def all_info(self):
"""Prints out all data and graphs available about the Environment."""
self.prints.all()
self.plots.all()
# TODO: Create a better .json format and allow loading a class from it.
[docs]
def export_environment(self, filename="environment"):
"""Export important attributes of Environment class to a ``.json`` file,
saving the information needed to recreate the same environment using
the ``custom_atmosphere`` model.
Parameters
----------
filename : string
The name of the file to be saved, without the extension.
"""
pressure = self.pressure.source
temperature = self.temperature.source
wind_x = self.wind_velocity_x.source
wind_y = self.wind_velocity_y.source
export_env_dictionary = {
"gravity": self.gravity(self.elevation),
"date": [
self.datetime_date.year,
self.datetime_date.month,
self.datetime_date.day,
self.datetime_date.hour,
],
"latitude": self.latitude,
"longitude": self.longitude,
"elevation": self.elevation,
"datum": self.datum,
"timezone": self.timezone,
"max_expected_height": float(self.max_expected_height),
"atmospheric_model_type": self.atmospheric_model_type,
"atmospheric_model_file": self.atmospheric_model_file,
"atmospheric_model_dict": self.atmospheric_model_dict,
"atmospheric_model_pressure_profile": pressure,
"atmospheric_model_temperature_profile": temperature,
"atmospheric_model_wind_velocity_x_profile": wind_x,
"atmospheric_model_wind_velocity_y_profile": wind_y,
}
with open(filename + ".json", "w") as f:
json.dump(export_env_dictionary, f, sort_keys=False, indent=4, default=str)
print(
f"Your Environment file was saved at '{filename}.json'. You can use "
"it in the future by using the custom_atmosphere atmospheric model."
)
[docs]
def set_earth_geometry(self, datum):
"""Sets the Earth geometry for the ``Environment`` class based on the
provided datum.
Parameters
----------
datum : str
The datum to be used for the Earth geometry. The following options
are supported: 'SIRGAS2000', 'SAD69', 'NAD83', 'WGS84'.
Returns
-------
earth_geometry : namedtuple
The namedtuple containing the Earth geometry.
"""
geodesy = namedtuple("earth_geometry", "semi_major_axis flattening")
ellipsoid = {
"SIRGAS2000": geodesy(6378137.0, 1 / 298.257223563),
"SAD69": geodesy(6378160.0, 1 / 298.25),
"NAD83": geodesy(6378137.0, 1 / 298.257024899),
"WGS84": geodesy(6378137.0, 1 / 298.257223563),
}
try:
return ellipsoid[datum]
except KeyError as e: # pragma: no cover
available_datums = ', '.join(ellipsoid.keys())
raise AttributeError(
f"The reference system '{datum}' is not recognized. Please use one of "
f"the following recognized datum: {available_datums}"
) from e
# Auxiliary functions - Geodesic Coordinates
[docs]
@staticmethod
def geodesic_to_utm(
lat, lon, semi_major_axis=6378137.0, flattening=1 / 298.257223563
):
"""Function which converts geodetic coordinates, i.e. lat/lon, to UTM
projection coordinates. Can be used only for latitudes between -80.00°
and 84.00°
Parameters
----------
lat : float
The latitude coordinates of the point of analysis, must be contained
between -80.00° and 84.00°
lon : float
The longitude coordinates of the point of analysis, must be
contained between -180.00° and 180.00°
semi_major_axis : float
The semi-major axis of the ellipsoid used to represent the Earth,
must be given in meters (default is 6,378,137.0 m, which corresponds
to the WGS84 ellipsoid)
flattening : float
The flattening of the ellipsoid used to represent the Earth, usually
between 1/250 and 1/150 (default is 1/298.257223563, which
corresponds to the WGS84 ellipsoid)
Returns
-------
x : float
East coordinate, always positive
y : float
North coordinate, always positive
utm_zone : int
The number of the UTM zone of the point of analysis, can vary
between 1 and 60
utm_letter : string
The letter of the UTM zone of the point of analysis, can vary
between C and X, omitting the letters "I" and "O"
hemis : string
Returns "S" for southern hemisphere and "N" for Northern hemisphere
EW : string
Returns "W" for western hemisphere and "E" for eastern hemisphere
"""
return geodesic_to_utm_tools(lat, lon, semi_major_axis, flattening)
[docs]
@staticmethod
def utm_to_geodesic(
x, y, utm_zone, hemis, semi_major_axis=6378137.0, flattening=1 / 298.257223563
):
"""Function to convert UTM coordinates to geodesic coordinates
(i.e. latitude and longitude).
Parameters
----------
x : float
East UTM coordinate in meters
y : float
North UTM coordinate in meters
utm_zone : int
The number of the UTM zone of the point of analysis, can vary
between 1 and 60
hemis : string
Equals to "S" for southern hemisphere and "N" for Northern
hemisphere
semi_major_axis : float
The semi-major axis of the ellipsoid used to represent the Earth,
must be given in meters (default is 6,378,137.0 m, which corresponds
to the WGS84 ellipsoid)
flattening : float
The flattening of the ellipsoid used to represent the Earth, usually
between 1/250 and 1/150 (default is 1/298.257223563, which
corresponds to the WGS84 ellipsoid)
Returns
-------
lat : float
latitude of the analyzed point
lon : float
latitude of the analyzed point
"""
return utm_to_geodesic_tools(x, y, utm_zone, hemis, semi_major_axis, flattening)
[docs]
@staticmethod
def calculate_earth_radius(
lat, semi_major_axis=6378137.0, flattening=1 / 298.257223563
):
"""Function to calculate the Earth's radius at a specific latitude
based on ellipsoidal reference model. The Earth radius here is
assumed as the distance between the ellipsoid's center of gravity and a
point on ellipsoid surface at the desired latitude.
Parameters
----------
lat : float
latitude at which the Earth radius will be calculated
semi_major_axis : float
The semi-major axis of the ellipsoid used to represent the Earth,
must be given in meters (default is 6,378,137.0 m, which corresponds
to the WGS84 ellipsoid)
flattening : float
The flattening of the ellipsoid used to represent the Earth, usually
between 1/250 and 1/150 (default is 1/298.257223563, which
corresponds to the WGS84 ellipsoid)
Returns
-------
radius : float
Earth radius at the desired latitude, in meters
Notes
-----
The ellipsoid is an approximation for the Earth model and
will result in an estimate of the perfect distance between
Earth's relief and its center of gravity.
"""
semi_minor_axis = semi_major_axis * (1 - flattening)
# Numpy trigonometric functions work with radians, so convert to radians
lat = lat * np.pi / 180
# Calculate the Earth Radius in meters
e_radius = np.sqrt(
(
(np.cos(lat) * (semi_major_axis**2)) ** 2
+ (np.sin(lat) * (semi_minor_axis**2)) ** 2
)
/ (
(np.cos(lat) * semi_major_axis) ** 2
+ (np.sin(lat) * semi_minor_axis) ** 2
)
)
return e_radius
[docs]
@staticmethod
def decimal_degrees_to_arc_seconds(angle):
"""Function to convert an angle in decimal degrees to degrees, arc
minutes and arc seconds.
Parameters
----------
angle : float
The angle that you need convert. Must be given in decimal degrees.
Returns
-------
degrees : int
The degrees.
arc_minutes : int
The arc minutes. 1 arc-minute = (1/60)*degree
arc_seconds : float
The arc Seconds. 1 arc-second = (1/3600)*degree
Examples
--------
Convert 45.5 degrees to degrees, arc minutes and arc seconds:
>>> from rocketpy import Environment
>>> Environment.decimal_degrees_to_arc_seconds(45.5)
(45, 30, 0.0)
"""
sign = -1 if angle < 0 else 1
degrees = int(abs(angle)) * sign
remainder = abs(angle) - abs(degrees)
arc_minutes = int(remainder * 60)
arc_seconds = (remainder * 60 - arc_minutes) * 60
return degrees, arc_minutes, arc_seconds
def to_dict(self, include_outputs=False):
env_dict = {
"gravity": self.gravity,
"date": self.date,
"latitude": self.latitude,
"longitude": self.longitude,
"elevation": self.elevation,
"datum": self.datum,
"timezone": self.timezone,
"max_expected_height": self.max_expected_height,
"atmospheric_model_type": self.atmospheric_model_type,
"pressure": self.pressure,
"temperature": self.temperature,
"wind_velocity_x": self.wind_velocity_x,
"wind_velocity_y": self.wind_velocity_y,
"wind_heading": self.wind_heading,
"wind_direction": self.wind_direction,
"wind_speed": self.wind_speed,
}
if include_outputs:
env_dict["density"] = self.density
env_dict["barometric_height"] = self.barometric_height
env_dict["speed_of_sound"] = self.speed_of_sound
env_dict["dynamic_viscosity"] = self.dynamic_viscosity
return env_dict
@classmethod
def from_dict(cls, data): # pylint: disable=too-many-statements
env = cls(
gravity=data["gravity"],
date=data["date"],
latitude=data["latitude"],
longitude=data["longitude"],
elevation=data["elevation"],
datum=data["datum"],
timezone=data["timezone"],
max_expected_height=data["max_expected_height"],
)
atmospheric_model = data["atmospheric_model_type"]
if atmospheric_model == "standard_atmosphere":
env.set_atmospheric_model("standard_atmosphere")
elif atmospheric_model == "custom_atmosphere":
env.set_atmospheric_model(
type="custom_atmosphere",
pressure=data["pressure"],
temperature=data["temperature"],
wind_u=data["wind_velocity_x"],
wind_v=data["wind_velocity_y"],
)
else:
env.__set_pressure_function(data["pressure"])
env.__set_temperature_function(data["temperature"])
env.__set_wind_velocity_x_function(data["wind_velocity_x"])
env.__set_wind_velocity_y_function(data["wind_velocity_y"])
env.__set_wind_heading_function(data["wind_heading"])
env.__set_wind_direction_function(data["wind_direction"])
env.__set_wind_speed_function(data["wind_speed"])
env.elevation = data["elevation"]
env.max_expected_height = data["max_expected_height"]
if atmospheric_model in ("windy", "forecast", "reanalysis", "ensemble"):
env.atmospheric_model_init_date = data["atmospheric_model_init_date"]
env.atmospheric_model_end_date = data["atmospheric_model_end_date"]
env.atmospheric_model_interval = data["atmospheric_model_interval"]
env.atmospheric_model_init_lat = data["atmospheric_model_init_lat"]
env.atmospheric_model_end_lat = data["atmospheric_model_end_lat"]
env.atmospheric_model_init_lon = data["atmospheric_model_init_lon"]
env.atmospheric_model_end_lon = data["atmospheric_model_end_lon"]
if atmospheric_model == "ensemble":
env.level_ensemble = data["level_ensemble"]
env.height_ensemble = data["height_ensemble"]
env.temperature_ensemble = data["temperature_ensemble"]
env.wind_u_ensemble = data["wind_u_ensemble"]
env.wind_v_ensemble = data["wind_v_ensemble"]
env.wind_heading_ensemble = data["wind_heading_ensemble"]
env.wind_direction_ensemble = data["wind_direction_ensemble"]
env.wind_speed_ensemble = data["wind_speed_ensemble"]
env.num_ensemble_members = data["num_ensemble_members"]
env.__reset_barometric_height_function()
env.calculate_density_profile()
env.calculate_speed_of_sound_profile()
env.calculate_dynamic_viscosity()
return env
if __name__ == "__main__": # pragma: no cover
import doctest
results = doctest.testmod()
if results.failed < 1:
print(f"All the {results.attempted} tests passed!")
else:
print(f"{results.failed} out of {results.attempted} tests failed.")