diff --git a/CHANGELOG.md b/CHANGELOG.md index cc5bc989e..d71dbb101 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -32,7 +32,20 @@ Attention: The newest changes should be on top --> ### Added -- ENH: Air brakes controller functions now support 8-parameter signature [#854](https://github.com/RocketPy-Team/RocketPy/pull/854) +- + +### Changed + +- + +### Fixed + +- + +## [v1.12.0] - 2026-03-08 + +### Added + - TST: Add acceptance tests for 3DOF flight simulation based on Bella Lui rocket [#914] (https://github.com/RocketPy-Team/RocketPy/pull/914_ - ENH: Add background map auto download functionality to Monte Carlo plots [#896](https://github.com/RocketPy-Team/RocketPy/pull/896) - MNT: net thrust addition to 3 dof in flight class [#907] (https://github.com/RocketPy-Team/RocketPy/pull/907) @@ -58,6 +71,8 @@ Attention: The newest changes should be on top --> ### Fixed +- BUG: Migrate Forecasts to UCAR THREDDS [#943](https://github.com/RocketPy-Team/RocketPy/pull/943) +- BUG: Fix hard-coded radius value for parachute added mass calculation [#889](https://github.com/RocketPy-Team/RocketPy/pull/889) - DOC: Fix documentation build [#908](https://github.com/RocketPy-Team/RocketPy/pull/908) - BUG: energy_data plot not working for 3 dof sims [[#906](https://github.com/RocketPy-Team/RocketPy/issues/906)] - BUG: Fix CSV column header spacing in FlightDataExporter [#864](https://github.com/RocketPy-Team/RocketPy/issues/864) diff --git a/docs/conf.py b/docs/conf.py index ae8a4b17d..e535082e7 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -22,12 +22,12 @@ # -- Project information ----------------------------------------------------- project = "RocketPy" -copyright = "2025, RocketPy Team" +copyright = "2026, RocketPy Team" author = "RocketPy Team" # The full version, including alpha/beta/rc tags -release = "1.11.0" +release = "1.12.0" # -- General configuration --------------------------------------------------- diff --git a/docs/user/environment/1-atm-models/ensemble.rst b/docs/user/environment/1-atm-models/ensemble.rst index 97c247f68..504cbfe60 100644 --- a/docs/user/environment/1-atm-models/ensemble.rst +++ b/docs/user/environment/1-atm-models/ensemble.rst @@ -1,3 +1,5 @@ +.. _ensemble_atmosphere: + Ensemble ======== @@ -21,7 +23,21 @@ Ensemble Forecast Global Ensemble Forecast System (GEFS) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ -The ``GEFS`` model is a global ensemble forecast model ... +.. danger:: + + **GEFS shortcut unavailable**: ``file="GEFS"`` is currently disabled in + RocketPy because NOMADS OPeNDAP is deactivated for this endpoint. + +.. note:: + + If you have a GEFS-compatible NetCDF or OPeNDAP dataset from another + provider (or a local copy), you can still load it explicitly by passing the + dataset path/URL in ``file`` and a compatible mapping in ``dictionary``. + + +The ``GEFS`` model is a global ensemble forecast system useful for uncertainty +analysis, but RocketPy's automatic ``file="GEFS"`` shortcut is temporarily +disabled. .. code-block:: python @@ -71,20 +87,16 @@ CMC Ensemble resulted in a change of the model's endpoint. Efforts are underway to \ restore access to the CMC Ensemble model as swiftly as possible. -.. code-block:: python +At the moment, there is no built-in ``file="CMC"`` shortcut in +``Environment.set_atmospheric_model``. - env_cmc = Environment( - date=date_info, - latitude=-21.960641, - longitude=-47.482122, - elevation=640, - ) - env_cmc.set_atmospheric_model(type="Ensemble", file="CMC") - env_cmc.all_info() +If you have a CMC-compatible NetCDF or OPeNDAP dataset, load it explicitly by +passing the dataset path/URL in ``file`` and a matching mapping dictionary in +``dictionary``. Ensemble Reanalysis ------------------- Ensemble reanalyses are also possible with RocketPy. See the -:ref:`reanalysis_ensemble` section for more information. +:ref:`reanalysis_ensemble` section for more information. \ No newline at end of file diff --git a/docs/user/environment/1-atm-models/forecast.rst b/docs/user/environment/1-atm-models/forecast.rst index c88c71ff2..ac91504e0 100644 --- a/docs/user/environment/1-atm-models/forecast.rst +++ b/docs/user/environment/1-atm-models/forecast.rst @@ -24,7 +24,7 @@ Global Forecast System (GFS) Using the latest forecast from GFS is simple. Set the atmospheric model to ``forecast`` and specify that GFS is the file you want. -Note that since data is downloaded from the NOMADS server, this line of code can +Note that since data is downloaded from a remote OPeNDAP server, this line of code can take longer than usual. .. jupyter-execute:: @@ -111,36 +111,15 @@ The same coordinates for SpacePort America will be used. High Resolution Window (HIRESW) ------------------------------- -The High Resolution Window (HIRESW) model is a sophisticated weather forecasting -system that operates at a high spatial resolution of approximately 3 km. -It utilizes two main dynamical cores: the Advanced Research WRF (WRF-ARW) and -the Finite Volume Cubed Sphere (FV3), each designed to enhance the accuracy of -weather predictions. +.. danger:: -You can easily set up HIRESW in RocketPy by specifying the date, latitude, and -longitude of your location. Let's use SpacePort America as an example. + **HIRESW shortcut unavailable**: ``file="HIRESW"`` is currently disabled in + RocketPy because NOMADS OPeNDAP is deactivated for this endpoint. -.. jupyter-execute:: - - env_hiresw = Environment( - date=tomorrow, - latitude=32.988528, - longitude=-106.975056, - ) +If you have a HIRESW-compatible dataset from another provider (or a local copy), +you can still load it explicitly by passing the path/URL in ``file`` and an +appropriate mapping in ``dictionary``. - env_hiresw.set_atmospheric_model( - type="Forecast", - file="HIRESW", - dictionary="HIRESW", - ) - - env_hiresw.plots.atmospheric_model() - -.. note:: - - The HRES model is updated every 12 hours, providing forecasts with a \ - resolution of 3 km. The model can predict weather conditions up to 48 hours \ - in advance. RocketPy uses the CONUS domain with ARW core. Using Windy Atmosphere @@ -248,6 +227,5 @@ Also, the servers may be down or may face high traffic. .. seealso:: - To see a complete list of available models on the NOAA's NOMADS server, visit - `NOMADS `_. - + To browse available NCEP model collections on UCAR THREDDS, visit + `THREDDS NCEP Catalog `_. \ No newline at end of file diff --git a/docs/user/environment/1-atm-models/soundings.rst b/docs/user/environment/1-atm-models/soundings.rst index 9a276477e..279750df5 100644 --- a/docs/user/environment/1-atm-models/soundings.rst +++ b/docs/user/environment/1-atm-models/soundings.rst @@ -57,31 +57,22 @@ This service allows users to download virtual soundings from numerical weather prediction models such as GFS, RAP, and NAM, and also real soundings from the Integrated Global Radiosonde Archive (IGRA). -These options can be retrieved as a text file in GSD format. -By generating such a file through the link above, the file's URL can be used to -import the atmospheric data into RocketPy. - -We will use the same sounding station as we did for the Wyoming Soundings. +These options can be retrieved as a text file in GSD format. However, +RocketPy no longer provides a dedicated ``set_atmospheric_model`` type for +NOAA RUC Soundings. .. note:: Select ROABs as the initial data source, specify the station through its \ WMO-ID, and opt for the ASCII (GSD format) button. -Initialize a new Environment instance: - -.. code-block:: python +If you need to use RUC-sounding-like data in RocketPy, convert it to one of the +supported workflows: - url = r"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" - - env = Environment() - env.set_atmospheric_model(type="NOAARucSounding", file=url) - env.plots.atmospheric_model() +- Use :ref:`custom_atmosphere` after parsing the text data. +- Use :ref:`reanalysis` or :ref:`forecast` with NetCDF/OPeNDAP sources. .. note:: The leading `r` in the URL string is used to indicate a raw string, which \ - is useful when dealing with backslashes in URLs. - - - + is useful when dealing with backslashes in URLs. \ No newline at end of file diff --git a/docs/user/environment/1-atm-models/standard_atmosphere.rst b/docs/user/environment/1-atm-models/standard_atmosphere.rst index 0c125dfd8..d6c1de782 100644 --- a/docs/user/environment/1-atm-models/standard_atmosphere.rst +++ b/docs/user/environment/1-atm-models/standard_atmosphere.rst @@ -1,3 +1,5 @@ +.. _standard_atmosphere: + Standard Atmosphere =================== @@ -29,4 +31,4 @@ The International Standard Atmosphere can also be reset at any time by using the .. jupyter-execute:: - env.set_atmospheric_model(type="standard_atmosphere") + env.set_atmospheric_model(type="standard_atmosphere") \ No newline at end of file diff --git a/docs/user/environment/3-further/other_apis.rst b/docs/user/environment/3-further/other_apis.rst index c70fd58f7..01d4b9a30 100644 --- a/docs/user/environment/3-further/other_apis.rst +++ b/docs/user/environment/3-further/other_apis.rst @@ -1,3 +1,5 @@ +.. _environment_other_apis: + Connecting to other APIs ======================== @@ -25,14 +27,19 @@ the following dimensions and variables: - Latitude - Longitude - Pressure Levels +- Temperature (as a function of Time, Pressure Levels, Latitude and Longitude) - Geopotential Height (as a function of Time, Pressure Levels, Latitude and Longitude) +- or Geopotential (as a function of Time, Pressure Levels, Latitude and Longitude) - Surface Geopotential Height (as a function of Time, Latitude and Longitude) + (optional) - Wind - U Component (as a function of Time, Pressure Levels, Latitude and Longitude) - Wind - V Component (as a function of Time, Pressure Levels, Latitude and Longitude) +Some projected grids also require a ``projection`` key in the mapping. + -For example, let's imagine we want to use the HIRESW model from this endpoint: -`https://nomads.ncep.noaa.gov/dods/hiresw/ `_ +For example, let's imagine we want to use a forecast model available via an +OPeNDAP endpoint. Looking through the variable list in the link above, we find the following correspondence: @@ -72,15 +79,85 @@ Therefore, we can create an environment like this: dictionary=name_mapping, ) +Built-in mapping dictionaries +----------------------------- + +Instead of a custom dictionary, you can pass a built-in mapping name in the +``dictionary`` argument. Common options include: + +- ``"ECMWF"`` +- ``"ECMWF_v0"`` +- ``"NOAA"`` +- ``"GFS"`` +- ``"NAM"`` +- ``"RAP"`` +- ``"HIRESW"`` (mapping available; latest-model shortcut currently disabled) +- ``"GEFS"`` (mapping available; latest-model shortcut currently disabled) +- ``"MERRA2"`` +- ``"CMC"`` (for compatible datasets loaded explicitly) + +What a mapping name means +^^^^^^^^^^^^^^^^^^^^^^^^^ + +- Base mapping names (for example ``"GFS"``, ``"NAM"`` and ``"RAP"``) map + RocketPy weather keys to the current default variable naming used by the + corresponding provider datasets. +- These defaults are aligned with current shortcut workflows (for example, + THREDDS-backed latest model sources) and may use projected coordinates + (``x``/``y`` plus ``projection``) depending on the model. + +Legacy mapping names +^^^^^^^^^^^^^^^^^^^^ + +If you are loading archived or older NOMADS-style datasets, use the explicit +legacy aliases: + +- ``"GFS_LEGACY"`` +- ``"NAM_LEGACY"`` +- ``"NOAA_LEGACY"`` +- ``"RAP_LEGACY"`` +- ``"CMC_LEGACY"`` +- ``"GEFS_LEGACY"`` +- ``"HIRESW_LEGACY"`` +- ``"MERRA2_LEGACY"`` + +Legacy aliases primarily cover older variable naming patterns such as +``lev``, ``tmpprs``, ``hgtprs``, ``ugrdprs`` and ``vgrdprs``. + +.. note:: + + Mapping names are case-insensitive. For example, + ``"gfs_legacy"`` and ``"GFS_LEGACY"`` are equivalent. + +For custom dictionaries, the canonical structure is: + +.. code-block:: python + + mapping = { + "time": "time", + "latitude": "lat", + "longitude": "lon", + "level": "lev", + "temperature": "tmpprs", + "surface_geopotential_height": "hgtsfc", # optional + "geopotential_height": "hgtprs", # or geopotential + "geopotential": None, + "u_wind": "ugrdprs", + "v_wind": "vgrdprs", + } + +.. important:: + + Ensemble datasets require an additional key for member selection: + ``"ensemble": ""``. + .. caution:: - Notice the ``file`` argument were suppressed in the code above. This is because \ - the URL depends on the date you are running the simulation. For example, as \ - it for now, a possible link could be: https://nomads.ncep.noaa.gov/dods/hiresw/hiresw20240803/hiresw_conusfv3_12z \ - (for the 3rd of August, 2024, at 12:00 UTC). \ - You should replace the date in the URL with the date you are running the simulation. \ - Different models may have different URL structures, so be sure to check the \ - documentation of the model you are using. + The ``file`` argument was intentionally omitted in the example above. This is + because the URL depends on the provider, dataset, and date you are running + the simulation. Build the endpoint according to the provider specification + and always validate that the target service is active before running your + simulation workflow. Without OPeNDAP protocol @@ -94,4 +171,3 @@ Environment class, for example: - `Meteomatics `_: `#545 `_ - `Open-Meteo `_: `#520 `_ - diff --git a/pyproject.toml b/pyproject.toml index 35ea34382..b9433c6d3 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "rocketpy" -version = "1.11.0" +version = "1.12.0" description="Advanced 6-DOF trajectory simulation for High-Power Rocketry." dynamic = ["dependencies"] readme = "README.md" diff --git a/rocketpy/environment/environment.py b/rocketpy/environment/environment.py index 6743b06ae..39441ecae 100644 --- a/rocketpy/environment/environment.py +++ b/rocketpy/environment/environment.py @@ -27,6 +27,7 @@ find_latitude_index, find_longitude_index, find_time_index, + geodesic_to_lambert_conformal, geodesic_to_utm, get_elevation_data_from_dataset, get_final_date_from_time_array, @@ -138,15 +139,15 @@ class Environment: 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``. + ``wyoming_sounding``, ``windy``, ``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`` + defined for ``wyoming_sounding``, ``windy``, ``forecast``, + ``reanalysis``, ``ensemble`` Environment.atmospheric_model_dict : dictionary Dictionary used to properly interpret ``netCDF`` and ``OPeNDAP`` files. - Only defined for ``Forecast``, ``Reanalysis``, ``Ensemble``. + 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 @@ -295,21 +296,21 @@ def __init__( - :attr:`Environment.datetime_date`: UTC time of launch. - Must be given if a Forecast, Reanalysis - or Ensemble, will be set as an atmospheric model. + Must be given if a ``windy``, ``forecast``, ``reanalysis`` + or ``ensemble`` atmospheric model will be used. 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 + launch location. Must be given if a ``windy``, ``forecast``, + ``reanalysis`` or ``ensemble`` atmospheric model will be used 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 + launch location. Must be given if a ``windy``, ``forecast``, + ``reanalysis`` or ``ensemble`` atmospheric model will be used 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. @@ -605,13 +606,81 @@ def __set_earth_rotation_vector(self): # Validators (used to verify an attribute is being set correctly.) + @staticmethod + def __dictionary_matches_dataset(dictionary, dataset): + """Check whether a mapping dictionary is compatible with a dataset.""" + variables = dataset.variables + required_keys = ( + "time", + "latitude", + "longitude", + "level", + "temperature", + "u_wind", + "v_wind", + ) + + for key in required_keys: + variable_name = dictionary.get(key) + if variable_name is None or variable_name not in variables: + return False + + projection_name = dictionary.get("projection") + if projection_name is not None and projection_name not in variables: + return False + + geopotential_height_name = dictionary.get("geopotential_height") + geopotential_name = dictionary.get("geopotential") + has_geopotential_height = ( + geopotential_height_name is not None + and geopotential_height_name in variables + ) + has_geopotential = ( + geopotential_name is not None and geopotential_name in variables + ) + + return has_geopotential_height or has_geopotential + + def __resolve_dictionary_for_dataset(self, dictionary, dataset): + """Resolve a compatible mapping dictionary for the loaded dataset. + + If the provided mapping is incompatible with the dataset variables, + this method tries built-in mappings and falls back to the first + compatible one. + """ + if self.__dictionary_matches_dataset(dictionary, dataset): + return dictionary + + for model_name, candidate in self.__weather_model_map.all_dictionaries.items(): + if self.__dictionary_matches_dataset(candidate, dataset): + warnings.warn( + "Provided weather mapping does not match dataset variables. " + f"Falling back to built-in mapping '{model_name}'." + ) + return candidate + + return dictionary + def __validate_dictionary(self, file, dictionary): # removed CMC until it is fixed. - available_models = ["GFS", "NAM", "RAP", "HIRESW", "GEFS", "ERA5", "MERRA2"] + available_models = [ + "GFS", + "NAM", + "RAP", + "HIRESW", + "GEFS", + "ERA5", + "MERRA2", + ] if isinstance(dictionary, str): dictionary = self.__weather_model_map.get(dictionary) - elif file in available_models: - dictionary = self.__weather_model_map.get(file) + elif isinstance(file, str): + matching_model = next( + (model for model in available_models if model.lower() == file.lower()), + None, + ) + if matching_model is not None: + dictionary = self.__weather_model_map.get(matching_model) if not isinstance(dictionary, dict): raise TypeError( "Please specify a dictionary or choose a valid model from the " @@ -1045,171 +1114,41 @@ def set_atmospheric_model( # pylint: disable=too-many-statements 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. + """Define the atmospheric model for this Environment. 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``, ``UCAR`` and ``MERRA2`` 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", - } + Atmospheric model selector (case-insensitive). Accepted values are + ``"standard_atmosphere"``, ``"wyoming_sounding"``, ``"windy"``, + ``"forecast"``, ``"reanalysis"``, ``"ensemble"`` and + ``"custom_atmosphere"``. + file : string | netCDF4.Dataset, optional + Data source or model shortcut. Meaning depends on ``type``: + + - ``"standard_atmosphere"`` and ``"custom_atmosphere"``: ignored. + - ``"wyoming_sounding"``: URL of the sounding text page. + - ``"windy"``: one of ``"ECMWF"``, ``"GFS"``, ``"ICON"`` or + ``"ICONEU"``. + - ``"forecast"``: local path, OPeNDAP URL, open + ``netCDF4.Dataset``, or one of ``"GFS"``, ``"NAM"`` or ``"RAP"`` + for the latest available forecast. + - ``"reanalysis"``: local path, OPeNDAP URL, or open + ``netCDF4.Dataset``. + - ``"ensemble"``: local path, OPeNDAP URL, open + ``netCDF4.Dataset``, or ``"GEFS"`` for the latest available + forecast. + dictionary : dict | str, optional + Variable-name mapping for ``"forecast"``, ``"reanalysis"`` and + ``"ensemble"``. It may be a custom dictionary or a built-in + mapping name (for example: ``"ECMWF"``, ``"ECMWF_v0"``, + ``"NOAA"``, ``"GFS"``, ``"NAM"``, ``"RAP"``, ``"HIRESW"``, + ``"GEFS"``, ``"MERRA2"`` or ``"CMC"``). + + If ``dictionary`` is omitted and ``file`` is one of RocketPy's + latest-model shortcuts, the matching built-in mapping is selected + automatically. For ensemble datasets, the mapping must include the + ensemble dimension key (typically ``"ensemble"``). pressure : float, string, array, callable, optional This defines the atmospheric pressure profile. @@ -1272,6 +1211,36 @@ def set_atmospheric_model( # pylint: disable=too-many-statements Returns ------- None + + Raises + ------ + ValueError + If ``type`` is unknown, if required launch date/time information is + missing for date-dependent models, if Windy model names are invalid, + or if required atmospheric variables cannot be read from the input + dataset. + TypeError + If ``dictionary`` is invalid for ``"forecast"``, ``"reanalysis"`` + or ``"ensemble"``. + KeyError + If a built-in mapping name passed in ``dictionary`` is unknown. + + See Also + -------- + :ref:`atmospheric_models` + Overview of all atmospheric-model workflows in the user guide. + :ref:`forecast` + Forecast and Windy usage details, including latest-model shortcuts. + :ref:`reanalysis` + Reanalysis and MERRA-2 examples. + :ref:`soundings` + Wyoming sounding workflow and RUC migration notes. + :ref:`custom_atmosphere` + Defining pressure, temperature and wind profiles directly. + :ref:`ensemble_atmosphere` + Ensemble forecasts and member-selection workflow. + :ref:`environment_other_apis` + Building custom mapping dictionaries for NetCDF/OPeNDAP APIs. """ # Save atmospheric model type self.atmospheric_model_type = type @@ -1287,6 +1256,36 @@ def set_atmospheric_model( # pylint: disable=too-many-statements case "windy": self.process_windy_atmosphere(file) case "forecast" | "reanalysis" | "ensemble": + if isinstance(file, str): + shortcut_map = self.__atm_type_file_to_function_map.get(type, {}) + matching_shortcut = next( + ( + shortcut + for shortcut in shortcut_map + if shortcut.lower() == file.lower() + ), + None, + ) + if matching_shortcut is not None: + file = matching_shortcut + + if isinstance(file, str): + file_upper = file.upper() + if type == "forecast" and file_upper == "HIRESW": + raise ValueError( + "The HIRESW latest-model shortcut is currently " + "unavailable because NOMADS OPeNDAP is deactivated. " + "Please use another forecast source or provide a " + "compatible dataset path/URL explicitly." + ) + if type == "ensemble" and file_upper == "GEFS": + raise ValueError( + "The GEFS latest-model shortcut is currently " + "unavailable because NOMADS OPeNDAP is deactivated. " + "Please use another ensemble source or provide a " + "compatible dataset path/URL explicitly." + ) + dictionary = self.__validate_dictionary(file, dictionary) try: fetch_function = self.__atm_type_file_to_function_map[type][file] @@ -1471,6 +1470,12 @@ def process_windy_atmosphere(self, model="ECMWF"): # pylint: disable=too-many-s ``ECMWF`` for the `ECMWF-HRES` model, ``GFS`` for the `GFS` model, ``ICON`` for the `ICON-Global` model or ``ICONEU`` for the `ICON-EU` model. + + Raises + ------ + ValueError + If ``model`` is not one of ``ECMWF``, ``GFS``, ``ICON`` or + ``ICONEU``. """ if model.lower() not in ["ecmwf", "gfs", "icon", "iconeu"]: @@ -1728,6 +1733,13 @@ def process_forecast_reanalysis(self, file, dictionary): # pylint: disable=too- Returns ------- None + + Raises + ------ + ValueError + If launch date/time was not set before loading date-dependent data, + or if required geopotential/geopotential-height, temperature, + wind-u, or wind-v variables cannot be read from the dataset. """ # Check if date, lat and lon are known self.__validate_datetime() @@ -1735,20 +1747,34 @@ def process_forecast_reanalysis(self, file, dictionary): # pylint: disable=too- # 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 + dictionary = self.__resolve_dictionary_for_dataset(dictionary, data) + # 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() + lon_array = data.variables[dictionary["longitude"]] + lat_array = data.variables[dictionary["latitude"]] + + # Some THREDDS datasets use projected x/y coordinates. + if dictionary.get("projection") is not None: + projection_variable = data.variables[dictionary["projection"]] + x_units = getattr(lon_array, "units", "m") + target_lon, target_lat = geodesic_to_lambert_conformal( + self.latitude, + self.longitude, + projection_variable, + x_units=x_units, + ) + else: + target_lon = self.longitude + target_lat = self.latitude # 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) + lon, lon_index = find_longitude_index(target_lon, lon_array) + _, lat_index = find_latitude_index(target_lat, lat_array) # Get pressure level data from file levels = get_pressure_levels_from_file(data, dictionary) @@ -1806,9 +1832,9 @@ def process_forecast_reanalysis(self, file, dictionary): # pylint: disable=too- ) 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] + x, y = target_lat, lon + x1, y1 = float(lat_array[lat_index - 1]), float(lon_array[lon_index - 1]) + x2, y2 = float(lat_array[lat_index]), float(lon_array[lon_index]) # Determine properties in lat, lon height = bilinear_interpolation( @@ -1860,6 +1886,17 @@ def process_forecast_reanalysis(self, file, dictionary): # pylint: disable=too- wind_vs[:, 1, 1], ) + # Some datasets expose different level counts between fields + # (e.g., temperature on isobaric1 and geopotential on isobaric). + min_profile_length = min( + len(levels), len(height), len(temper), len(wind_u), len(wind_v) + ) + levels = levels[:min_profile_length] + height = height[:min_profile_length] + temper = temper[:min_profile_length] + wind_u = wind_u[:min_profile_length] + wind_v = wind_v[:min_profile_length] + # Determine wind speed, heading and direction wind_speed = calculate_wind_speed(wind_u, wind_v) wind_heading = calculate_wind_heading(wind_u, wind_v) @@ -1917,14 +1954,14 @@ def process_forecast_reanalysis(self, file, dictionary): # pylint: disable=too- ) 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] + self.atmospheric_model_init_lat = float(lat_array[0]) + self.atmospheric_model_end_lat = float(lat_array[len(lat_array) - 1]) + self.atmospheric_model_init_lon = float(lon_array[0]) + self.atmospheric_model_end_lon = float(lon_array[len(lon_array) - 1]) # Save debugging data - self.lat_array = lat_list - self.lon_array = lon_list + self.lat_array = [x1, x2] + self.lon_array = [y1, y2] self.lon_index = lon_index self.lat_index = lat_index self.geopotentials = geopotentials @@ -1932,7 +1969,10 @@ def process_forecast_reanalysis(self, file, dictionary): # pylint: disable=too- self.wind_vs = wind_vs self.levels = levels self.temperatures = temperatures - self.time_array = time_array[:].tolist() + self.time_array = [ + float(time_array[0]), + float(time_array[time_array.shape[0] - 1]), + ] self.height = height # Close weather data @@ -1994,6 +2034,13 @@ def process_ensemble(self, file, dictionary): # pylint: disable=too-many-locals -------- See the :class:``rocketpy.environment.weather_model_mapping`` for some dictionary examples. + + Raises + ------ + ValueError + If launch date/time was not set before loading date-dependent data, + or if required geopotential/geopotential-height, temperature, + wind-u, or wind-v variables cannot be read from the dataset. """ # Check if date, lat and lon are known self.__validate_datetime() @@ -2004,23 +2051,40 @@ def process_ensemble(self, file, dictionary): # pylint: disable=too-many-locals else: data = file + dictionary = self.__resolve_dictionary_for_dataset(dictionary, data) + # 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() + lon_array = data.variables[dictionary["longitude"]] + lat_array = data.variables[dictionary["latitude"]] + + # Some THREDDS datasets use projected x/y coordinates. + # TODO CHECK THIS I AM NOT SURE????? + if dictionary.get("projection") is not None: + projection_variable = data.variables[dictionary["projection"]] + x_units = getattr(lon_array, "units", "m") + target_lon, target_lat = geodesic_to_lambert_conformal( + self.latitude, + self.longitude, + projection_variable, + x_units=x_units, + ) + else: + target_lon = self.longitude + target_lat = self.latitude # 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) + lon, lon_index = find_longitude_index(target_lon, lon_array) + _, lat_index = find_latitude_index(target_lat, lat_array) # Get ensemble data from file + has_ensemble_dimension = True 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 + except KeyError: + has_ensemble_dimension = False + num_members = 1 # Get pressure level data from file levels = get_pressure_levels_from_file(data, dictionary) @@ -2079,10 +2143,16 @@ def process_ensemble(self, file, dictionary): # pylint: disable=too-many-locals "Unable to read wind-v component. Check file and dictionary." ) from e + if not has_ensemble_dimension: + geopotentials = np.expand_dims(geopotentials, axis=0) + temperatures = np.expand_dims(temperatures, axis=0) + wind_us = np.expand_dims(wind_us, axis=0) + wind_vs = np.expand_dims(wind_vs, axis=0) + # 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] + x, y = target_lat, lon + x1, y1 = float(lat_array[lat_index - 1]), float(lon_array[lon_index - 1]) + x2, y2 = float(lat_array[lat_index]), float(lon_array[lon_index]) # Determine properties in lat, lon height = bilinear_interpolation( @@ -2134,6 +2204,19 @@ def process_ensemble(self, file, dictionary): # pylint: disable=too-many-locals wind_vs[:, :, 1, 1], ) + min_profile_length = min( + len(levels), + height.shape[1], + temper.shape[1], + wind_u.shape[1], + wind_v.shape[1], + ) + levels = levels[:min_profile_length] + height = height[:, :min_profile_length] + temper = temper[:, :min_profile_length] + wind_u = wind_u[:, :min_profile_length] + wind_v = wind_v[:, :min_profile_length] + # Determine wind speed, heading and direction wind_speed = calculate_wind_speed(wind_u, wind_v) wind_heading = calculate_wind_heading(wind_u, wind_v) @@ -2166,14 +2249,14 @@ def process_ensemble(self, file, dictionary): # pylint: disable=too-many-locals 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] + self.atmospheric_model_init_lat = float(lat_array[0]) + self.atmospheric_model_end_lat = float(lat_array[len(lat_array) - 1]) + self.atmospheric_model_init_lon = float(lon_array[0]) + self.atmospheric_model_end_lon = float(lon_array[len(lon_array) - 1]) # Save debugging data - self.lat_array = lat_list - self.lon_array = lon_list + self.lat_array = [x1, x2] + self.lon_array = [y1, y2] self.lon_index = lon_index self.lat_index = lat_index self.geopotentials = geopotentials @@ -2181,7 +2264,10 @@ def process_ensemble(self, file, dictionary): # pylint: disable=too-many-locals self.wind_vs = wind_vs self.levels = levels self.temperatures = temperatures - self.time_array = time_array[:].tolist() + self.time_array = [ + float(time_array[0]), + float(time_array[time_array.shape[0] - 1]), + ] self.height = height # Close weather data diff --git a/rocketpy/environment/fetchers.py b/rocketpy/environment/fetchers.py index d5ac2a1df..589159f1c 100644 --- a/rocketpy/environment/fetchers.py +++ b/rocketpy/environment/fetchers.py @@ -113,33 +113,18 @@ def fetch_gfs_file_return_dataset(max_attempts=10, base_delay=2): RuntimeError If unable to load the latest weather data for GFS. """ - time_attempt = datetime.now(tz=timezone.utc) + file_url = ( + "https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/GFS/Global_0p25deg/Best" + ) attempt_count = 0 - dataset = None - - # TODO: the code below is trying to determine the hour of the latest available - # forecast by trial and error. This is not the best way to do it. We should - # actually check the NOAA website for the latest forecast time. Refactor needed. while attempt_count < max_attempts: - time_attempt -= timedelta(hours=6) # GFS updates every 6 hours - file_url = ( - f"https://nomads.ncep.noaa.gov/dods/gfs_0p25/gfs" - f"{time_attempt.year:04d}{time_attempt.month:02d}" - f"{time_attempt.day:02d}/" - f"gfs_0p25_{6 * (time_attempt.hour // 6):02d}z" - ) try: - # Attempts to create a dataset from the file using OpenDAP protocol. - dataset = netCDF4.Dataset(file_url) - return dataset + return netCDF4.Dataset(file_url) except OSError: attempt_count += 1 time.sleep(base_delay**attempt_count) - if dataset is None: - raise RuntimeError( - "Unable to load latest weather data for GFS through " + file_url - ) + raise RuntimeError("Unable to load latest weather data for GFS through " + file_url) def fetch_nam_file_return_dataset(max_attempts=10, base_delay=2): @@ -163,28 +148,16 @@ def fetch_nam_file_return_dataset(max_attempts=10, base_delay=2): RuntimeError If unable to load the latest weather data for NAM. """ - # Attempt to get latest forecast - time_attempt = datetime.now(tz=timezone.utc) + file_url = "https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/NAM/CONUS_12km/Best" attempt_count = 0 - dataset = None - while attempt_count < max_attempts: - time_attempt -= timedelta(hours=6) # NAM updates every 6 hours - file = ( - f"https://nomads.ncep.noaa.gov/dods/nam/nam{time_attempt.year:04d}" - f"{time_attempt.month:02d}{time_attempt.day:02d}/" - f"nam_conusnest_{6 * (time_attempt.hour // 6):02d}z" - ) try: - # Attempts to create a dataset from the file using OpenDAP protocol. - dataset = netCDF4.Dataset(file) - return dataset + return netCDF4.Dataset(file_url) except OSError: attempt_count += 1 time.sleep(base_delay**attempt_count) - if dataset is None: - raise RuntimeError("Unable to load latest weather data for NAM through " + file) + raise RuntimeError("Unable to load latest weather data for NAM through " + file_url) def fetch_rap_file_return_dataset(max_attempts=10, base_delay=2): @@ -208,28 +181,88 @@ def fetch_rap_file_return_dataset(max_attempts=10, base_delay=2): RuntimeError If unable to load the latest weather data for RAP. """ - # Attempt to get latest forecast - time_attempt = datetime.now(tz=timezone.utc) + file_url = "https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/RAP/CONUS_13km/Best" attempt_count = 0 - dataset = None + while attempt_count < max_attempts: + try: + return netCDF4.Dataset(file_url) + except OSError: + attempt_count += 1 + time.sleep(base_delay**attempt_count) + + raise RuntimeError("Unable to load latest weather data for RAP through " + file_url) + +def fetch_hrrr_file_return_dataset(max_attempts=10, base_delay=2): + """Fetches the latest HRRR (High-Resolution Rapid Refresh) dataset from + the NOAA's GrADS data server using the OpenDAP protocol. + + Parameters + ---------- + max_attempts : int, optional + The maximum number of attempts to fetch the dataset. Default is 10. + base_delay : int, optional + The base delay in seconds between attempts. Default is 2. + + Returns + ------- + netCDF4.Dataset + The HRRR dataset. + + Raises + ------ + RuntimeError + If unable to load the latest weather data for HRRR. + """ + file_url = "https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/HRRR/CONUS_2p5km/Best" + attempt_count = 0 while attempt_count < max_attempts: - time_attempt -= timedelta(hours=1) # RAP updates every hour - file = ( - f"https://nomads.ncep.noaa.gov/dods/rap/rap{time_attempt.year:04d}" - f"{time_attempt.month:02d}{time_attempt.day:02d}/" - f"rap_{time_attempt.hour:02d}z" - ) try: - # Attempts to create a dataset from the file using OpenDAP protocol. - dataset = netCDF4.Dataset(file) - return dataset + return netCDF4.Dataset(file_url) except OSError: attempt_count += 1 time.sleep(base_delay**attempt_count) - if dataset is None: - raise RuntimeError("Unable to load latest weather data for RAP through " + file) + raise RuntimeError( + "Unable to load latest weather data for HRRR through " + file_url + ) + + +def fetch_aigfs_file_return_dataset(max_attempts=10, base_delay=2): + """Fetches the latest AIGFS (Artificial Intelligence GFS) dataset from + the NOAA's GrADS data server using the OpenDAP protocol. + + Parameters + ---------- + max_attempts : int, optional + The maximum number of attempts to fetch the dataset. Default is 10. + base_delay : int, optional + The base delay in seconds between attempts. Default is 2. + + Returns + ------- + netCDF4.Dataset + The AIGFS dataset. + + Raises + ------ + RuntimeError + If unable to load the latest weather data for AIGFS. + """ + file_url = ( + "https://thredds.ucar.edu/thredds/dodsC/grib/NCEP/AIGFS/Global_0p25deg/Best" + ) + attempt_count = 0 + while attempt_count < max_attempts: + try: + return netCDF4.Dataset(file_url) + except OSError: + attempt_count += 1 + time.sleep(base_delay**attempt_count) + + raise RuntimeError( + "Unable to load latest weather data for AIGFS through " + file_url + ) def fetch_hiresw_file_return_dataset(max_attempts=10, base_delay=2): diff --git a/rocketpy/environment/tools.py b/rocketpy/environment/tools.py index 1239ee6b9..fb0179c9e 100644 --- a/rocketpy/environment/tools.py +++ b/rocketpy/environment/tools.py @@ -5,7 +5,7 @@ future to improve their performance and usability. """ -import bisect +import math import warnings import netCDF4 @@ -109,6 +109,63 @@ def calculate_wind_speed(u, v, w=0.0): return np.sqrt(u**2 + v**2 + w**2) +def geodesic_to_lambert_conformal(lat, lon, projection_variable, x_units="m"): + """Convert geodesic coordinates to Lambert conformal projected coordinates. + + Parameters + ---------- + lat : float + Latitude in degrees. + lon : float + Longitude in degrees, ranging from -180 to 180. + projection_variable : netCDF4.Variable + Projection variable containing Lambert conformal metadata. + x_units : str, optional + Units used by the dataset x coordinate. Supported values are meters + and kilometers. Default is "m". + + Returns + ------- + tuple[float, float] + Projected coordinates ``(x, y)`` in the same units as ``x_units``. + """ + lat_radians = math.radians(lat) + lon_radians = math.radians(lon % 360) + + lat_origin = math.radians(float(projection_variable.latitude_of_projection_origin)) + lon_origin = math.radians(float(projection_variable.longitude_of_central_meridian)) + + standard_parallel = projection_variable.standard_parallel + if np.ndim(standard_parallel) == 0: + standard_parallels = [float(standard_parallel)] + else: + standard_parallels = np.asarray(standard_parallel, dtype=float).tolist() + + if len(standard_parallels) >= 2: + phi_1 = math.radians(standard_parallels[0]) + phi_2 = math.radians(standard_parallels[1]) + n = math.log(math.cos(phi_1) / math.cos(phi_2)) / math.log( + math.tan(math.pi / 4 + phi_2 / 2) / math.tan(math.pi / 4 + phi_1 / 2) + ) + else: + phi_1 = math.radians(standard_parallels[0]) + n = math.sin(phi_1) + + earth_radius = float(getattr(projection_variable, "earth_radius", 6371229.0)) + f_const = (math.cos(phi_1) * math.tan(math.pi / 4 + phi_1 / 2) ** n) / n + + rho = earth_radius * f_const / (math.tan(math.pi / 4 + lat_radians / 2) ** n) + rho_origin = earth_radius * f_const / (math.tan(math.pi / 4 + lat_origin / 2) ** n) + theta = n * (lon_radians - lon_origin) + + x_meters = rho * math.sin(theta) + y_meters = rho_origin - rho * math.cos(theta) + + if str(x_units).lower().startswith("km"): + return x_meters / 1000.0, y_meters / 1000.0 + return x_meters, y_meters + + ## These functions are meant to be used with netcdf4 datasets @@ -168,7 +225,7 @@ def mask_and_clean_dataset(*args): return data_array -def find_longitude_index(longitude, lon_list): +def find_longitude_index(longitude, lon_list): # pylint: disable=too-many-statements """Finds the index of the given longitude in a list of longitudes. Parameters @@ -188,30 +245,48 @@ def find_longitude_index(longitude, lon_list): ValueError If the longitude is not within the range covered by the list. """ - # Determine if file uses -180 to 180 or 0 to 360 - if lon_list[0] < 0 or lon_list[-1] < 0: - # Convert input to -180 - 180 - lon = longitude if longitude < 180 else -180 + longitude % 180 - else: - # Convert input to 0 - 360 - lon = longitude % 360 - # Check if reversed or sorted - if lon_list[0] < lon_list[-1]: - # Deal with sorted lon_list - lon_index = bisect.bisect(lon_list, lon) + + def _coord_value(source, index): + return float(source[index]) + + lon_len = len(lon_list) + lon_start = _coord_value(lon_list, 0) + lon_end = _coord_value(lon_list, lon_len - 1) + + # Determine if file uses geographic longitudes in [-180, 180] or [0, 360]. + # Do not remap projected x coordinates. + is_geographic_longitude = abs(lon_start) <= 360 and abs(lon_end) <= 360 + if is_geographic_longitude: + if lon_start < 0 or lon_end < 0: + lon = longitude if longitude < 180 else -180 + longitude % 180 + else: + lon = longitude % 360 else: - # Deal with reversed lon_list - lon_list.reverse() - lon_index = len(lon_list) - bisect.bisect_left(lon_list, lon) - lon_list.reverse() + lon = longitude + + is_ascending = lon_start < lon_end + + # Binary search to find the insertion index such that index-1 and index + # bracket the requested longitude. + low = 0 + high = lon_len + while low < high: + mid = (low + high) // 2 + mid_value = _coord_value(lon_list, mid) + if (mid_value < lon) if is_ascending else (mid_value > lon): + low = mid + 1 + else: + high = mid + lon_index = low + # Take care of longitude value equal to maximum longitude in the grid - if lon_index == len(lon_list) and lon_list[lon_index - 1] == lon: - lon_index = lon_index - 1 + if lon_index == lon_len and _coord_value(lon_list, lon_index - 1) == lon: + lon_index -= 1 # Check if longitude value is inside the grid - if lon_index == 0 or lon_index == len(lon_list): + if lon_index in (0, lon_len): raise ValueError( f"Longitude {lon} not inside region covered by file, which is " - f"from {lon_list[0]} to {lon_list[-1]}." + f"from {lon_start} to {lon_end}." ) return lon, lon_index @@ -237,28 +312,39 @@ def find_latitude_index(latitude, lat_list): ValueError If the latitude is not within the range covered by the list. """ - # Check if reversed or sorted - if lat_list[0] < lat_list[-1]: - # Deal with sorted lat_list - lat_index = bisect.bisect(lat_list, latitude) - else: - # Deal with reversed lat_list - lat_list.reverse() - lat_index = len(lat_list) - bisect.bisect_left(lat_list, latitude) - lat_list.reverse() - # Take care of latitude value equal to maximum longitude in the grid - if lat_index == len(lat_list) and lat_list[lat_index - 1] == latitude: - lat_index = lat_index - 1 + + def _coord_value(source, index): + return float(source[index]) + + lat_len = len(lat_list) + lat_start = _coord_value(lat_list, 0) + lat_end = _coord_value(lat_list, lat_len - 1) + is_ascending = lat_start < lat_end + + low = 0 + high = lat_len + while low < high: + mid = (low + high) // 2 + mid_value = _coord_value(lat_list, mid) + if (mid_value < latitude) if is_ascending else (mid_value > latitude): + low = mid + 1 + else: + high = mid + lat_index = low + + # Take care of latitude value equal to maximum latitude in the grid + if lat_index == lat_len and _coord_value(lat_list, lat_index - 1) == latitude: + lat_index -= 1 # Check if latitude value is inside the grid - if lat_index == 0 or lat_index == len(lat_list): + if lat_index in (0, lat_len): raise ValueError( f"Latitude {latitude} not inside region covered by file, " - f"which is from {lat_list[0]} to {lat_list[-1]}." + f"which is from {lat_start} to {lat_end}." ) return latitude, lat_index -def find_time_index(datetime_date, time_array): +def find_time_index(datetime_date, time_array): # pylint: disable=too-many-statements """Finds the index of the given datetime in a netCDF4 time array. Parameters @@ -280,26 +366,58 @@ def find_time_index(datetime_date, time_array): ValueError If the exact datetime is not available and the nearest datetime is used instead. """ - time_index = netCDF4.date2index( - datetime_date, time_array, calendar="gregorian", select="nearest" - ) - # Convert times do dates and numbers - input_time_num = netCDF4.date2num( - datetime_date, time_array.units, calendar="gregorian" - ) - file_time_num = time_array[time_index] - file_time_date = netCDF4.num2date( - time_array[time_index], time_array.units, calendar="gregorian" - ) + time_len = len(time_array) + time_units = time_array.units + input_time_num = netCDF4.date2num(datetime_date, time_units, calendar="gregorian") + + first_time_num = float(time_array[0]) + last_time_num = float(time_array[time_len - 1]) + is_ascending = first_time_num <= last_time_num + + # Binary search nearest index using scalar probing only. + low = 0 + high = time_len + while low < high: + mid = (low + high) // 2 + mid_time_num = float(time_array[mid]) + if ( + (mid_time_num < input_time_num) + if is_ascending + else (mid_time_num > input_time_num) + ): + low = mid + 1 + else: + high = mid + + right_index = min(max(low, 0), time_len - 1) + left_index = min(max(right_index - 1, 0), time_len - 1) + + right_time_num = float(time_array[right_index]) + left_time_num = float(time_array[left_index]) + if abs(input_time_num - left_time_num) <= abs(right_time_num - input_time_num): + time_index = left_index + file_time_num = left_time_num + else: + time_index = right_index + file_time_num = right_time_num + + file_time_date = netCDF4.num2date(file_time_num, time_units, calendar="gregorian") + # Check if time is inside range supplied by file - if time_index == 0 and input_time_num < file_time_num: + if time_index == 0 and ( + (is_ascending and input_time_num < file_time_num) + or (not is_ascending and input_time_num > file_time_num) + ): raise ValueError( f"The chosen launch time '{datetime_date.strftime('%Y-%m-%d-%H:')} UTC' is" " not available in the provided file. Please choose a time within the range" " of the file, which starts at " f"'{file_time_date.strftime('%Y-%m-%d-%H')} UTC'." ) - elif time_index == len(time_array) - 1 and input_time_num > file_time_num: + elif time_index == time_len - 1 and ( + (is_ascending and input_time_num > file_time_num) + or (not is_ascending and input_time_num < file_time_num) + ): raise ValueError( "Chosen launch time is not available in the provided file, " f"which ends at {file_time_date}." diff --git a/rocketpy/environment/weather_model_mapping.py b/rocketpy/environment/weather_model_mapping.py index 75089f577..c490fad9d 100644 --- a/rocketpy/environment/weather_model_mapping.py +++ b/rocketpy/environment/weather_model_mapping.py @@ -1,9 +1,42 @@ class WeatherModelMapping: - """Class to map the weather model variables to the variables used in the - Environment class. + """Map provider-specific variable names to RocketPy weather fields. + + RocketPy reads forecast/reanalysis/ensemble datasets using canonical keys + such as ``time``, ``latitude``, ``longitude``, ``level``, ``temperature``, + ``geopotential_height``, ``geopotential``, ``u_wind`` and ``v_wind``. + Each dictionary in this class maps those canonical keys to the actual + variable names in a specific data provider format. + + Mapping families + ---------------- + - Base names (for example ``GFS``, ``NAM``, ``RAP``) represent the current + default mappings used by the latest-model shortcuts and THREDDS-style + datasets. + - ``*_LEGACY`` names represent older NOMADS-style variable naming + conventions (for example ``lev``, ``tmpprs``, ``ugrdprs`` and + ``vgrdprs``) and are intended for archived or previously downloaded files. + + Notes + ----- + - Mappings can also include optional keys such as ``projection`` for + projected grids and ``ensemble`` for member dimensions. + - The :meth:`get` method is case-insensitive, so ``"gfs_legacy"`` and + ``"GFS_LEGACY"`` are equivalent. """ GFS = { + "time": "time", + "latitude": "lat", + "longitude": "lon", + "level": "isobaric", + "temperature": "Temperature_isobaric", + "surface_geopotential_height": "Geopotential_height_surface", + "geopotential_height": "Geopotential_height_isobaric", + "geopotential": None, + "u_wind": "u-component_of_wind_isobaric", + "v_wind": "v-component_of_wind_isobaric", + } + GFS_LEGACY = { "time": "time", "latitude": "lat", "longitude": "lon", @@ -16,6 +49,19 @@ class WeatherModelMapping: "v_wind": "vgrdprs", } NAM = { + "time": "time", + "latitude": "y", + "longitude": "x", + "projection": "LambertConformal_Projection", + "level": "isobaric", + "temperature": "Temperature_isobaric", + "surface_geopotential_height": None, + "geopotential_height": "Geopotential_height_isobaric", + "geopotential": None, + "u_wind": "u-component_of_wind_isobaric", + "v_wind": "v-component_of_wind_isobaric", + } + NAM_LEGACY = { "time": "time", "latitude": "lat", "longitude": "lon", @@ -54,6 +100,18 @@ class WeatherModelMapping: "v_wind": "v", } NOAA = { + "time": "time", + "latitude": "lat", + "longitude": "lon", + "level": "isobaric", + "temperature": "Temperature_isobaric", + "surface_geopotential_height": "Geopotential_height_surface", + "geopotential_height": "Geopotential_height_isobaric", + "geopotential": None, + "u_wind": "u-component_of_wind_isobaric", + "v_wind": "v-component_of_wind_isobaric", + } + NOAA_LEGACY = { "time": "time", "latitude": "lat", "longitude": "lon", @@ -66,6 +124,19 @@ class WeatherModelMapping: "v_wind": "vgrdprs", } RAP = { + "time": "time", + "latitude": "y", + "longitude": "x", + "projection": "LambertConformal_Projection", + "level": "isobaric", + "temperature": "Temperature_isobaric", + "surface_geopotential_height": None, + "geopotential_height": "Geopotential_height_isobaric", + "geopotential": None, + "u_wind": "u-component_of_wind_isobaric", + "v_wind": "v-component_of_wind_isobaric", + } + RAP_LEGACY = { "time": "time", "latitude": "lat", "longitude": "lon", @@ -90,6 +161,19 @@ class WeatherModelMapping: "u_wind": "ugrdprs", "v_wind": "vgrdprs", } + CMC_LEGACY = { + "time": "time", + "latitude": "lat", + "longitude": "lon", + "level": "lev", + "ensemble": "ens", + "temperature": "tmpprs", + "surface_geopotential_height": None, + "geopotential_height": "hgtprs", + "geopotential": None, + "u_wind": "ugrdprs", + "v_wind": "vgrdprs", + } GEFS = { "time": "time", "latitude": "lat", @@ -103,6 +187,19 @@ class WeatherModelMapping: "u_wind": "ugrdprs", "v_wind": "vgrdprs", } + GEFS_LEGACY = { + "time": "time", + "latitude": "lat", + "longitude": "lon", + "level": "lev", + "ensemble": "ens", + "temperature": "tmpprs", + "surface_geopotential_height": None, + "geopotential_height": "hgtprs", + "geopotential": None, + "u_wind": "ugrdprs", + "v_wind": "vgrdprs", + } HIRESW = { "time": "time", "latitude": "lat", @@ -114,6 +211,17 @@ class WeatherModelMapping: "u_wind": "ugrdprs", "v_wind": "vgrdprs", } + HIRESW_LEGACY = { + "time": "time", + "latitude": "lat", + "longitude": "lon", + "level": "lev", + "temperature": "tmpprs", + "surface_geopotential_height": "hgtsfc", + "geopotential_height": "hgtprs", + "u_wind": "ugrdprs", + "v_wind": "vgrdprs", + } MERRA2 = { "time": "time", "latitude": "lat", @@ -127,29 +235,78 @@ class WeatherModelMapping: "u_wind": "U", "v_wind": "V", } + MERRA2_LEGACY = { + "time": "time", + "latitude": "lat", + "longitude": "lon", + "level": "lev", + "temperature": "T", + "surface_geopotential_height": None, + "surface_geopotential": "PHIS", + "geopotential_height": "H", + "geopotential": None, + "u_wind": "U", + "v_wind": "V", + } def __init__(self): - """Initialize the class, creates a dictionary with all the weather models - available and their respective dictionaries with the variables.""" + """Build the lookup table with default and legacy mapping aliases.""" self.all_dictionaries = { "GFS": self.GFS, + "GFS_LEGACY": self.GFS_LEGACY, "NAM": self.NAM, + "NAM_LEGACY": self.NAM_LEGACY, "ECMWF_v0": self.ECMWF_v0, "ECMWF": self.ECMWF, "NOAA": self.NOAA, + "NOAA_LEGACY": self.NOAA_LEGACY, "RAP": self.RAP, + "RAP_LEGACY": self.RAP_LEGACY, "CMC": self.CMC, + "CMC_LEGACY": self.CMC_LEGACY, "GEFS": self.GEFS, + "GEFS_LEGACY": self.GEFS_LEGACY, "HIRESW": self.HIRESW, + "HIRESW_LEGACY": self.HIRESW_LEGACY, "MERRA2": self.MERRA2, + "MERRA2_LEGACY": self.MERRA2_LEGACY, } def get(self, model): + """Return a mapping dictionary by model alias (case-insensitive). + + Parameters + ---------- + model : str + Mapping alias name, such as ``"GFS"`` or ``"GFS_LEGACY"``. + + Returns + ------- + dict + Dictionary mapping RocketPy canonical weather keys to dataset + variable names. + + Raises + ------ + KeyError + If ``model`` is unknown or not a string. + """ + if not isinstance(model, str): + raise KeyError( + f"Model {model} not found in the WeatherModelMapping. " + f"The available models are: {self.all_dictionaries.keys()}" + ) + try: return self.all_dictionaries[model] - except KeyError as e: + except KeyError as exc: + model_casefold = model.casefold() + for key, value in self.all_dictionaries.items(): + if key.casefold() == model_casefold: + return value + raise KeyError( f"Model {model} not found in the WeatherModelMapping. " f"The available models are: {self.all_dictionaries.keys()}" - ) from e + ) from exc diff --git a/rocketpy/rocket/parachute.py b/rocketpy/rocket/parachute.py index 83b0ce0fd..4e0318d18 100644 --- a/rocketpy/rocket/parachute.py +++ b/rocketpy/rocket/parachute.py @@ -92,17 +92,25 @@ class Parachute: Function of noisy_pressure_signal. Parachute.clean_pressure_signal_function : Function Function of clean_pressure_signal. + Parachute.drag_coefficient : float + Drag coefficient of the inflated canopy shape, used only when + ``radius`` is not provided to estimate the parachute radius from + ``cd_s``: ``R = sqrt(cd_s / (drag_coefficient * pi))``. Typical + values: 1.4 for hemispherical canopies (default), 0.75 for flat + circular canopies, 1.5 for extended-skirt canopies. Parachute.radius : float Length of the non-unique semi-axis (radius) of the inflated hemispheroid - parachute in meters. - Parachute.height : float, None + parachute in meters. If not provided at construction time, it is + estimated from ``cd_s`` and ``drag_coefficient``. + Parachute.height : float Length of the unique semi-axis (height) of the inflated hemispheroid parachute in meters. Parachute.porosity : float - Geometric porosity of the canopy (ratio of open area to total canopy area), - in [0, 1]. Affects only the added-mass scaling during descent; it does - not change ``cd_s`` (drag). The default, 0.0432, yields an added-mass - of 1.0 (“neutral” behavior). + Geometric porosity of the canopy (ratio of open area to total canopy + area), in [0, 1]. Affects only the added-mass scaling during descent; + it does not change ``cd_s`` (drag). The default value of 0.0432 is + chosen so that the resulting ``added_mass_coefficient`` equals + approximately 1.0 ("neutral" added-mass behavior). Parachute.added_mass_coefficient : float Coefficient used to calculate the added-mass due to dragged air. It is calculated from the porosity of the parachute. @@ -116,9 +124,10 @@ def __init__( sampling_rate, lag=0, noise=(0, 0, 0), - radius=1.5, + radius=None, height=None, porosity=0.0432, + drag_coefficient=1.4, ): """Initializes Parachute class. @@ -172,25 +181,83 @@ def __init__( passed to the trigger function. Default value is ``(0, 0, 0)``. Units are in Pa. radius : float, optional - Length of the non-unique semi-axis (radius) of the inflated hemispheroid - parachute. Default value is 1.5. + Length of the non-unique semi-axis (radius) of the inflated + hemispheroid parachute. If not provided, it is estimated from + ``cd_s`` and ``drag_coefficient`` using: + ``radius = sqrt(cd_s / (drag_coefficient * pi))``. Units are in meters. height : float, optional Length of the unique semi-axis (height) of the inflated hemispheroid parachute. Default value is the radius of the parachute. Units are in meters. porosity : float, optional - Geometric porosity of the canopy (ratio of open area to total canopy area), - in [0, 1]. Affects only the added-mass scaling during descent; it does - not change ``cd_s`` (drag). The default, 0.0432, yields an added-mass - of 1.0 (“neutral” behavior). + Geometric porosity of the canopy (ratio of open area to total + canopy area), in [0, 1]. Affects only the added-mass scaling + during descent; it does not change ``cd_s`` (drag). The default + value of 0.0432 is chosen so that the resulting + ``added_mass_coefficient`` equals approximately 1.0 ("neutral" + added-mass behavior). + drag_coefficient : float, optional + Drag coefficient of the inflated canopy shape, used only when + ``radius`` is not provided. It relates the aerodynamic ``cd_s`` + to the physical canopy area via + ``cd_s = drag_coefficient * pi * radius**2``. Typical values: + + - **1.4** — hemispherical canopy (default, NASA SP-8066) + - **0.75** — flat circular canopy + - **1.5** — extended-skirt canopy + + Has no effect when ``radius`` is explicitly provided. """ + + # Save arguments as attributes self.name = name self.cd_s = cd_s self.trigger = trigger self.sampling_rate = sampling_rate self.lag = lag self.noise = noise + self.drag_coefficient = drag_coefficient + self.porosity = porosity + + # Initialize derived attributes + self.radius = self.__resolve_radius(radius, cd_s, drag_coefficient) + self.height = self.__resolve_height(height, self.radius) + self.added_mass_coefficient = self.__compute_added_mass_coefficient( + self.porosity + ) + self.__init_noise(noise) + self.__evaluate_trigger_function(trigger) + + # Prints and plots + self.prints = _ParachutePrints(self) + + def __resolve_radius(self, radius, cd_s, drag_coefficient): + """Resolves parachute radius from input or aerodynamic relation.""" + if radius is not None: + return radius + + # cd_s = Cd * S = Cd * pi * R^2 => R = sqrt(cd_s / (Cd * pi)) + return np.sqrt(cd_s / (drag_coefficient * np.pi)) + + def __resolve_height(self, height, radius): + """Resolves parachute height defaulting to radius when not provided.""" + return height or radius + + def __compute_added_mass_coefficient(self, porosity): + """Computes the added-mass coefficient from canopy porosity.""" + return 1.068 * ( + 1 - 1.465 * porosity - 0.25975 * porosity**2 + 1.2626 * porosity**3 + ) + + def __init_noise(self, noise): + """Initializes all noise-related attributes. + + Parameters + ---------- + noise : tuple, list + List in the format (mean, standard deviation, time-correlation). + """ self.noise_signal = [[-1e-6, np.random.normal(noise[0], noise[1])]] self.noisy_pressure_signal = [] self.clean_pressure_signal = [] @@ -200,32 +267,19 @@ def __init__( self.clean_pressure_signal_function = Function(0) self.noisy_pressure_signal_function = Function(0) self.noise_signal_function = Function(0) - self.radius = radius - self.height = height or radius - self.porosity = porosity - self.added_mass_coefficient = 1.068 * ( - 1 - - 1.465 * self.porosity - - 0.25975 * self.porosity**2 - + 1.2626 * self.porosity**3 - ) - alpha, beta = self.noise_corr self.noise_function = lambda: ( alpha * self.noise_signal[-1][1] + beta * np.random.normal(noise[0], noise[1]) ) - self.prints = _ParachutePrints(self) - - self.__evaluate_trigger_function(trigger) - def __evaluate_trigger_function(self, trigger): """This is used to set the triggerfunc attribute that will be used to interact with the Flight class. """ # pylint: disable=unused-argument, function-redefined - # The parachute is deployed by a custom function + + # Case 1: The parachute is deployed by a custom function if callable(trigger): # work around for having added sensors to parachute triggers # to avoid breaking changes @@ -238,9 +292,10 @@ def triggerfunc(p, h, y, sensors): self.triggerfunc = triggerfunc + # Case 2: The parachute is deployed at a given height elif isinstance(trigger, (int, float)): # The parachute is deployed at a given height - def triggerfunc(p, h, y, sensors): # pylint: disable=unused-argument + def triggerfunc(p, h, y, sensors): # p = pressure considering parachute noise signal # h = height above ground level considering parachute noise signal # y = [x, y, z, vx, vy, vz, e0, e1, e2, e3, w1, w2, w3] @@ -248,9 +303,10 @@ def triggerfunc(p, h, y, sensors): # pylint: disable=unused-argument self.triggerfunc = triggerfunc + # Case 3: The parachute is deployed at apogee elif trigger.lower() == "apogee": # The parachute is deployed at apogee - def triggerfunc(p, h, y, sensors): # pylint: disable=unused-argument + def triggerfunc(p, h, y, sensors): # p = pressure considering parachute noise signal # h = height above ground level considering parachute noise signal # y = [x, y, z, vx, vy, vz, e0, e1, e2, e3, w1, w2, w3] @@ -258,6 +314,7 @@ def triggerfunc(p, h, y, sensors): # pylint: disable=unused-argument self.triggerfunc = triggerfunc + # Case 4: Invalid trigger input else: raise ValueError( f"Unable to set the trigger function for parachute '{self.name}'. " @@ -289,7 +346,7 @@ def info(self): def all_info(self): """Prints all information about the Parachute class.""" self.info() - # self.plots.all() # Parachutes still doesn't have plots + # self.plots.all() # TODO: Parachutes still doesn't have plots def to_dict(self, **kwargs): allow_pickle = kwargs.get("allow_pickle", True) @@ -309,6 +366,7 @@ def to_dict(self, **kwargs): "lag": self.lag, "noise": self.noise, "radius": self.radius, + "drag_coefficient": self.drag_coefficient, "height": self.height, "porosity": self.porosity, } @@ -341,7 +399,8 @@ def from_dict(cls, data): sampling_rate=data["sampling_rate"], lag=data["lag"], noise=data["noise"], - radius=data.get("radius", 1.5), + radius=data.get("radius", None), + drag_coefficient=data.get("drag_coefficient", 1.4), height=data.get("height", None), porosity=data.get("porosity", 0.0432), ) diff --git a/rocketpy/rocket/rocket.py b/rocketpy/rocket/rocket.py index 86fa981a9..51719753d 100644 --- a/rocketpy/rocket/rocket.py +++ b/rocketpy/rocket/rocket.py @@ -1502,9 +1502,10 @@ def add_parachute( sampling_rate=100, lag=0, noise=(0, 0, 0), - radius=1.5, + radius=None, height=None, porosity=0.0432, + drag_coefficient=1.4, ): """Creates a new parachute, storing its parameters such as opening delay, drag coefficients and trigger function. @@ -1564,26 +1565,34 @@ def add_parachute( passed to the trigger function. Default value is (0, 0, 0). Units are in pascal. radius : float, optional - Length of the non-unique semi-axis (radius) of the inflated hemispheroid - parachute. Default value is 1.5. + Length of the non-unique semi-axis (radius) of the inflated + hemispheroid parachute. If not provided, it is estimated from + `cd_s` and `drag_coefficient` using: + `radius = sqrt(cd_s / (drag_coefficient * pi))`. Units are in meters. height : float, optional Length of the unique semi-axis (height) of the inflated hemispheroid parachute. Default value is the radius of the parachute. Units are in meters. porosity : float, optional - Geometric porosity of the canopy (ratio of open area to total canopy area), - in [0, 1]. Affects only the added-mass scaling during descent; it does - not change ``cd_s`` (drag). The default, 0.0432, yields an added-mass - of 1.0 (“neutral” behavior). + Geometric porosity of the canopy (ratio of open area to total + canopy area), in [0, 1]. Affects only the added-mass scaling + during descent; it does not change `cd_s` (drag). The default + value of 0.0432 yields an `added_mass_coefficient` of + approximately 1.0 ("neutral" added-mass behavior). + drag_coefficient : float, optional + Drag coefficient of the inflated canopy shape, used only when + `radius` is not provided. Typical values: 1.4 for hemispherical + canopies (default), 0.75 for flat circular canopies, 1.5 for + extended-skirt canopies. Has no effect when `radius` is given. Returns ------- parachute : Parachute - Parachute containing trigger, sampling_rate, lag, cd_s, noise, radius, - height, porosity and name. Furthermore, it stores clean_pressure_signal, - noise_signal and noisyPressureSignal which are filled in during - Flight simulation. + Parachute containing trigger, sampling_rate, lag, cd_s, noise, + radius, drag_coefficient, height, porosity and name. Furthermore, + it stores clean_pressure_signal, noise_signal and + noisyPressureSignal which are filled in during Flight simulation. """ parachute = Parachute( name, @@ -1595,6 +1604,7 @@ def add_parachute( radius, height, porosity, + drag_coefficient, ) self.parachutes.append(parachute) return self.parachutes[-1] diff --git a/rocketpy/stochastic/stochastic_parachute.py b/rocketpy/stochastic/stochastic_parachute.py index dea8a077d..038907187 100644 --- a/rocketpy/stochastic/stochastic_parachute.py +++ b/rocketpy/stochastic/stochastic_parachute.py @@ -31,6 +31,9 @@ class StochasticParachute(StochasticModel): List with the name of the parachute object. This cannot be randomized. radius : tuple, list, int, float Radius of the parachute in meters. + drag_coefficient : tuple, list, int, float + Drag coefficient of the inflated canopy shape, used only when + ``radius`` is not provided. height : tuple, list, int, float Height of the parachute in meters. porosity : tuple, list, int, float @@ -46,6 +49,7 @@ def __init__( lag=None, noise=None, radius=None, + drag_coefficient=None, height=None, porosity=None, ): @@ -74,6 +78,9 @@ def __init__( time-correlation). radius : tuple, list, int, float Radius of the parachute in meters. + drag_coefficient : tuple, list, int, float + Drag coefficient of the inflated canopy shape, used only when + ``radius`` is not provided. height : tuple, list, int, float Height of the parachute in meters. porosity : tuple, list, int, float @@ -86,6 +93,7 @@ def __init__( self.lag = lag self.noise = noise self.radius = radius + self.drag_coefficient = drag_coefficient self.height = height self.porosity = porosity @@ -100,6 +108,7 @@ def __init__( noise=noise, name=None, radius=radius, + drag_coefficient=drag_coefficient, height=height, porosity=porosity, ) diff --git a/tests/integration/environment/test_environment.py b/tests/integration/environment/test_environment.py index 3bdd5209a..37078b8fd 100644 --- a/tests/integration/environment/test_environment.py +++ b/tests/integration/environment/test_environment.py @@ -178,8 +178,11 @@ def test_gefs_atmosphere(mock_show, example_spaceport_env): # pylint: disable=u example_spaceport_env : rocketpy.Environment Example environment object to be tested. """ - example_spaceport_env.set_atmospheric_model(type="Ensemble", file="GEFS") - assert example_spaceport_env.all_info() is None + with pytest.raises( + ValueError, + match="GEFS latest-model shortcut is currently unavailable", + ): + example_spaceport_env.set_atmospheric_model(type="Ensemble", file="GEFS") @pytest.mark.slow @@ -234,13 +237,15 @@ def test_hiresw_ensemble_atmosphere(mock_show, example_spaceport_env): # pylint example_spaceport_env.set_date(date_info) - example_spaceport_env.set_atmospheric_model( - type="Forecast", - file="HIRESW", - dictionary="HIRESW", - ) - - assert example_spaceport_env.all_info() is None + with pytest.raises( + ValueError, + match="HIRESW latest-model shortcut is currently unavailable", + ): + example_spaceport_env.set_atmospheric_model( + type="Forecast", + file="HIRESW", + dictionary="HIRESW", + ) @pytest.mark.skip(reason="CMC model is currently not working") diff --git a/tests/integration/simulation/test_flight.py b/tests/integration/simulation/test_flight.py index 7e25a8927..66f0848a4 100644 --- a/tests/integration/simulation/test_flight.py +++ b/tests/integration/simulation/test_flight.py @@ -717,6 +717,48 @@ def invalid_controller_9_params( # pylint: disable=unused-argument ) +def make_controller_test_environment_access(methods_called): + def _call_env_methods(environment, altitude_asl): + _ = environment.elevation + methods_called["elevation"] = True + _ = environment.wind_velocity_x(altitude_asl) + methods_called["wind_velocity_x"] = True + _ = environment.wind_velocity_y(altitude_asl) + methods_called["wind_velocity_y"] = True + _ = environment.speed_of_sound(altitude_asl) + methods_called["speed_of_sound"] = True + _ = environment.pressure(altitude_asl) + methods_called["pressure"] = True + _ = environment.temperature(altitude_asl) + methods_called["temperature"] = True + + def controller( # pylint: disable=unused-argument + time, + sampling_rate, + state, + state_history, + observed_variables, + air_brakes, + sensors, + environment, + ): + """Controller that tests access to various environment methods.""" + altitude_asl = state[2] + + if time < 3.9: + return None + + try: + _call_env_methods(environment, altitude_asl) + air_brakes.deployment_level = 0.3 + except AttributeError as e: + raise AssertionError(f"Environment method not accessible: {e}") from e + + return (time, air_brakes.deployment_level) + + return controller + + def test_environment_methods_accessible_in_controller( calisto_robust, example_plain_env ): @@ -742,54 +784,13 @@ def test_environment_methods_accessible_in_controller( "temperature": False, } - def controller_test_environment_access( # pylint: disable=unused-argument - time, - sampling_rate, - state, - state_history, - observed_variables, - air_brakes, - sensors, - environment, - ): - """Controller that tests access to various environment methods.""" - altitude_asl = state[2] - - if time < 3.9: - return None - - # Test accessing various environment methods - try: - _ = environment.elevation - methods_called["elevation"] = True - - _ = environment.wind_velocity_x(altitude_asl) - methods_called["wind_velocity_x"] = True - - _ = environment.wind_velocity_y(altitude_asl) - methods_called["wind_velocity_y"] = True - - _ = environment.speed_of_sound(altitude_asl) - methods_called["speed_of_sound"] = True - - _ = environment.pressure(altitude_asl) - methods_called["pressure"] = True - - _ = environment.temperature(altitude_asl) - methods_called["temperature"] = True - - air_brakes.deployment_level = 0.3 - except AttributeError as e: - # If any method is not accessible, the test should fail - raise AssertionError(f"Environment method not accessible: {e}") from e - - return (time, air_brakes.deployment_level) + controller = make_controller_test_environment_access(methods_called) # Add air brakes with environment-testing controller calisto_robust.parachutes = [] calisto_robust.add_air_brakes( drag_coefficient_curve="data/rockets/calisto/air_brakes_cd.csv", - controller_function=controller_test_environment_access, + controller_function=controller, sampling_rate=10, clamp=True, ) diff --git a/tests/unit/environment/test_environment.py b/tests/unit/environment/test_environment.py index 6ad3e51db..6d04c089f 100644 --- a/tests/unit/environment/test_environment.py +++ b/tests/unit/environment/test_environment.py @@ -7,6 +7,7 @@ from rocketpy import Environment from rocketpy.environment.tools import geodesic_to_utm, utm_to_geodesic +from rocketpy.environment.weather_model_mapping import WeatherModelMapping @pytest.mark.parametrize( @@ -243,3 +244,70 @@ def test_environment_export_environment_exports_valid_environment_json( ) os.remove("environment.json") + + +class _DummyDataset: + """Small test double that mimics a netCDF dataset variables mapping.""" + + def __init__(self, variable_names): + self.variables = {name: object() for name in variable_names} + + +def test_resolve_dictionary_keeps_compatible_mapping(example_plain_env): + """Keep the user-selected mapping when it already matches dataset keys.""" + gfs_mapping = example_plain_env._Environment__weather_model_map.get("GFS") + dataset = _DummyDataset( + [ + "time", + "lat", + "lon", + "isobaric", + "Temperature_isobaric", + "Geopotential_height_isobaric", + "u-component_of_wind_isobaric", + "v-component_of_wind_isobaric", + ] + ) + + resolved = example_plain_env._Environment__resolve_dictionary_for_dataset( + gfs_mapping, dataset + ) + + assert resolved is gfs_mapping + + +def test_resolve_dictionary_falls_back_to_legacy_mapping(example_plain_env): + """Fallback to a compatible built-in mapping for legacy NOMADS-style files.""" + thredds_gfs_mapping = example_plain_env._Environment__weather_model_map.get("GFS") + dataset = _DummyDataset( + [ + "time", + "lat", + "lon", + "lev", + "tmpprs", + "hgtprs", + "ugrdprs", + "vgrdprs", + ] + ) + + resolved = example_plain_env._Environment__resolve_dictionary_for_dataset( + thredds_gfs_mapping, dataset + ) + + # Explicit legacy mappings should be preferred over unrelated model mappings. + assert resolved == example_plain_env._Environment__weather_model_map.get( + "GFS_LEGACY" + ) + assert resolved["level"] == "lev" + assert resolved["temperature"] == "tmpprs" + assert resolved["geopotential_height"] == "hgtprs" + + +def test_weather_model_mapping_exposes_legacy_aliases(): + """Legacy mapping names should be available and case-insensitive.""" + mapping = WeatherModelMapping() + + assert mapping.get("GFS_LEGACY")["temperature"] == "tmpprs" + assert mapping.get("gfs_legacy")["temperature"] == "tmpprs" diff --git a/tests/unit/rocket/test_parachute.py b/tests/unit/rocket/test_parachute.py new file mode 100644 index 000000000..e193b777b --- /dev/null +++ b/tests/unit/rocket/test_parachute.py @@ -0,0 +1,111 @@ +"""Unit tests for the Parachute class, focusing on the radius and +drag_coefficient parameters introduced in PR #889.""" + +import numpy as np +import pytest + +from rocketpy import Parachute + + +def _make_parachute(**kwargs): + defaults = { + "name": "test", + "cd_s": 10.0, + "trigger": "apogee", + "sampling_rate": 100, + } + defaults.update(kwargs) + return Parachute(**defaults) + + +class TestParachuteRadiusEstimation: + """Tests for auto-computed radius from cd_s and drag_coefficient.""" + + def test_radius_auto_computed_from_cd_s_default_drag_coefficient(self): + """When radius is not provided the radius is estimated using the + default drag_coefficient of 1.4 and the formula R = sqrt(cd_s / (Cd * pi)).""" + cd_s = 10.0 + parachute = _make_parachute(cd_s=cd_s) + expected_radius = np.sqrt(cd_s / (1.4 * np.pi)) + assert parachute.radius == pytest.approx(expected_radius, rel=1e-9) + + def test_radius_auto_computed_uses_custom_drag_coefficient(self): + """When drag_coefficient is provided and radius is not, the radius + must be estimated using the given drag_coefficient.""" + cd_s = 10.0 + custom_cd = 0.75 + parachute = _make_parachute(cd_s=cd_s, drag_coefficient=custom_cd) + expected_radius = np.sqrt(cd_s / (custom_cd * np.pi)) + assert parachute.radius == pytest.approx(expected_radius, rel=1e-9) + + def test_explicit_radius_overrides_estimation(self): + """When radius is explicitly provided, it must be used directly and + drag_coefficient must be ignored for the radius calculation.""" + explicit_radius = 2.5 + parachute = _make_parachute(radius=explicit_radius, drag_coefficient=0.5) + assert parachute.radius == explicit_radius + + def test_drag_coefficient_stored_on_instance(self): + """drag_coefficient must be stored as an attribute regardless of + whether radius is provided or not.""" + parachute = _make_parachute(drag_coefficient=0.75) + assert parachute.drag_coefficient == 0.75 + + def test_drag_coefficient_default_is_1_4(self): + """Default drag_coefficient must be 1.4 for backward compatibility.""" + parachute = _make_parachute() + assert parachute.drag_coefficient == pytest.approx(1.4) + + def test_drogue_radius_smaller_than_main(self): + """A drogue (cd_s=1.0) must have a smaller radius than a main (cd_s=10.0) + when using the same drag_coefficient.""" + main = _make_parachute(cd_s=10.0) + drogue = _make_parachute(cd_s=1.0) + assert drogue.radius < main.radius + + def test_drogue_radius_approximately_0_48(self): + """For cd_s=1.0 and drag_coefficient=1.4, the estimated radius + must be approximately 0.48 m (fixes the previous hard-coded 1.5 m).""" + drogue = _make_parachute(cd_s=1.0) + assert drogue.radius == pytest.approx(0.476, abs=1e-3) + + def test_main_radius_approximately_1_51(self): + """For cd_s=10.0 and drag_coefficient=1.4, the estimated radius + must be approximately 1.51 m, matching the old hard-coded value.""" + main = _make_parachute(cd_s=10.0) + assert main.radius == pytest.approx(1.508, abs=1e-3) + + +class TestParachuteSerialization: + """Tests for to_dict / from_dict round-trip including drag_coefficient.""" + + def test_to_dict_includes_drag_coefficient(self): + """to_dict must include the drag_coefficient key.""" + parachute = _make_parachute(drag_coefficient=0.75) + data = parachute.to_dict() + assert "drag_coefficient" in data + assert data["drag_coefficient"] == 0.75 + + def test_from_dict_round_trip_preserves_drag_coefficient(self): + """A Parachute serialized to dict and restored must have the same + drag_coefficient.""" + original = _make_parachute(cd_s=5.0, drag_coefficient=0.75) + data = original.to_dict() + restored = Parachute.from_dict(data) + assert restored.drag_coefficient == pytest.approx(0.75) + assert restored.radius == pytest.approx(original.radius, rel=1e-9) + + def test_from_dict_defaults_drag_coefficient_to_1_4_when_absent(self): + """Dicts serialized before drag_coefficient was added (no key) must + fall back to 1.4 for backward compatibility.""" + data = { + "name": "legacy", + "cd_s": 10.0, + "trigger": "apogee", + "sampling_rate": 100, + "lag": 0, + "noise": (0, 0, 0), + # no drag_coefficient key — simulates old serialized data + } + parachute = Parachute.from_dict(data) + assert parachute.drag_coefficient == pytest.approx(1.4)