AbstractGlobalPressureTemperature.java

  1. /* Copyright 2002-2025 CS GROUP
  2.  * Licensed to CS GROUP (CS) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * CS licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *   http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.orekit.models.earth.weather;

  19. import java.io.BufferedInputStream;
  20. import java.io.IOException;
  21. import java.io.InputStream;

  23. import org.hipparchus.CalculusFieldElement;
  24. import org.orekit.bodies.FieldGeodeticPoint;
  25. import org.orekit.bodies.GeodeticPoint;
  26. import org.orekit.data.DataSource;
  27. import org.orekit.models.earth.troposphere.AzimuthalGradientCoefficients;
  28. import org.orekit.models.earth.troposphere.AzimuthalGradientProvider;
  29. import org.orekit.models.earth.troposphere.FieldAzimuthalGradientCoefficients;
  30. import org.orekit.models.earth.troposphere.FieldViennaACoefficients;
  31. import org.orekit.models.earth.troposphere.TroposphericModelUtils;
  32. import org.orekit.models.earth.troposphere.ViennaACoefficients;
  33. import org.orekit.models.earth.troposphere.ViennaAProvider;
  34. import org.orekit.time.AbsoluteDate;
  35. import org.orekit.time.FieldAbsoluteDate;

  37. /** Base class for Global Pressure and Temperature 2, 2w and 3 models.
  38.  * These models are empirical models that provide the temperature, the pressure and the water vapor pressure
  39.  * of a site depending its latitude and  longitude. These models also {@link ViennaACoefficients provide}
  40.  * the a<sub>h</sub> and a<sub>w</sub> coefficients for Vienna models.
  41.  * <p>
  42.  * The requisite coefficients for the computation of the weather parameters are provided by the
  43.  * Department of Geodesy and Geoinformation of the Vienna University. They are based on an
  44.  * external grid file like "gpt2_1.grd" (1° x 1°), "gpt2_5.grd" (5° x 5°), "gpt2_1w.grd" (1° x 1°),
  45.  * "gpt2_5w.grd" (5° x 5°), "gpt3_1.grd" (1° x 1°), or "gpt3_5.grd" (5° x 5°) available at:
  46.  * <a href="https://vmf.geo.tuwien.ac.at/codes/"> link</a>
  47.  * </p>
  48.  * <p>
  49.  * A bilinear interpolation is performed in order to obtained the correct values of the weather parameters.
  50.  * </p>
  51.  * <p>
  52.  * The format is always the same, with and example shown below for the pressure and the temperature.
  53.  * The "GPT2w" model (w stands for wet) also provide humidity parameters and the "GPT3" model also
  54.  * provides horizontal gradient, so the number of columns vary depending on the model.
  55.  * <p>
  56.  * Example:
  57.  * </p>
  58.  * <pre>
  59.  * %  lat    lon   p:a0    A1   B1   A2   B2  T:a0    A1   B1   A2   B2
  60.  *   87.5    2.5 101421    21  409 -217 -122 259.2 -13.2 -6.1  2.6  0.3
  61.  *   87.5    7.5 101416    21  411 -213 -120 259.3 -13.1 -6.1  2.6  0.3
  62.  *   87.5   12.5 101411    22  413 -209 -118 259.3 -13.1 -6.1  2.6  0.3
  63.  *   87.5   17.5 101407    23  415 -205 -116 259.4 -13.0 -6.1  2.6  0.3
  64.  *   ...
  65.  * </pre>
  66.  *
  67.  * @see "K. Lagler, M. Schindelegger, J. Böhm, H. Krasna, T. Nilsson (2013),
  68.  * GPT2: empirical slant delay model for radio space geodetic techniques. Geophys
  69.  * Res Lett 40(6):1069–1073. doi:10.1002/grl.50288"
  70.  *
  71.  * @author Bryan Cazabonne
  72.  * @author Luc Maisonobe
  73.  * @since 12.1
  74.  */
  75. public abstract class AbstractGlobalPressureTemperature
  76.     implements ViennaAProvider, AzimuthalGradientProvider, PressureTemperatureHumidityProvider {

  78.     /** Loaded grid. */
  79.     private final Grid grid;

  81.     /**
  82.      * Constructor with source of GPTn auxiliary data given by user.
  83.      *
  84.      * @param source grid data source
  85.      * @param expected expected seasonal models types
  86.      * @exception IOException if grid data cannot be read
  87.      */
  88.     protected AbstractGlobalPressureTemperature(final DataSource source, final SeasonalModelType... expected)
  89.         throws IOException {

  91.         // load the grid data
  92.         try (InputStream         is     = source.getOpener().openStreamOnce();
  93.              BufferedInputStream bis    = new BufferedInputStream(is)) {
  94.             final GptNParser     parser = new GptNParser(expected);
  95.             parser.loadData(bis, source.getName());
  96.             grid = parser.getGrid();
  97.         }

  99.     }

  101.     /** Get duration since reference date.
  102.      * @param date date
  103.      * @return duration since reference date
  104.      * @since 13.0
  105.      */
  106.     protected abstract double deltaRef(AbsoluteDate date);

  108.     /** Get duration since reference date.
  109.      * @param <T> type of the field elements
  110.      * @param date date
  111.      * @return duration since reference date
  112.      * @since 13.0
  113.      */
  114.     protected abstract <T extends CalculusFieldElement<T>> T deltaRef(FieldAbsoluteDate<T> date);

  116.     /** {@inheritDoc} */
  117.     @Override
  118.     public ViennaACoefficients getA(final GeodeticPoint location, final AbsoluteDate date) {

  120.         // set up interpolation parameters
  121.         final CellInterpolator interpolator =
  122.             grid.getInterpolator(location.getLatitude(), location.getLongitude(), location.getAltitude(),
  123.                                  deltaRef(date));

  125.         // ah and aw coefficients
  126.         return new ViennaACoefficients(interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.AH)) * 0.001,
  127.                                        interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.AW)) * 0.001);

  129.     }

  131.     /** {@inheritDoc} */
  132.     @Override
  133.     public PressureTemperatureHumidity getWeatherParameters(final GeodeticPoint location, final AbsoluteDate date) {

  135.         // set up interpolation parameters
  136.         final CellInterpolator interpolator =
  137.             grid.getInterpolator(location.getLatitude(), location.getLongitude(), location.getAltitude(),
  138.                                  deltaRef(date));

  140.         // Temperature [K]
  141.         final double temperature = interpolator.interpolate(EvaluatedGridEntry::getTemperature);

  143.         // Pressure [hPa]
  144.         final double pressure = interpolator.interpolate(EvaluatedGridEntry::getPressure);

  146.         // water vapor decrease factor
  147.         final double lambda = grid.hasModels(SeasonalModelType.LAMBDA) ?
  148.                               interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.LAMBDA)) :
  149.                               Double.NaN;

  151.         // Water vapor pressure [hPa]
  152.         final double el;
  153.         if (Double.isNaN(lambda)) {
  154.             // GPT2
  155.             final double qv = interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.QV)) * 0.001;
  156.             el = qv * pressure / (0.622 + 0.378 * qv);
  157.         } else {
  158.             // GPT2w and GPT3
  159.             el = interpolator.interpolate(EvaluatedGridEntry::getWaterVaporPressure);
  160.         }

  162.         // mean temperature weighted with water vapor pressure
  163.         final double tm = grid.hasModels(SeasonalModelType.TM) ?
  164.                           interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.TM)) :
  165.                           Double.NaN;

  167.         return new PressureTemperatureHumidity(location.getAltitude(),
  168.                                                TroposphericModelUtils.HECTO_PASCAL.toSI(pressure),
  169.                                                temperature,
  170.                                                TroposphericModelUtils.HECTO_PASCAL.toSI(el),
  171.                                                tm,
  172.                                                lambda);

  174.     }

  176.     /** {@inheritDoc} */
  177.     @Override
  178.     public AzimuthalGradientCoefficients getGradientCoefficients(final GeodeticPoint location, final AbsoluteDate date) {

  180.         if (grid.hasModels(SeasonalModelType.GN_H, SeasonalModelType.GE_H, SeasonalModelType.GN_W, SeasonalModelType.GE_W)) {
  181.             // set up interpolation parameters
  182.             final CellInterpolator interpolator = grid.getInterpolator(location.getLatitude(), location.getLongitude(),
  183.                                                                        location.getAltitude(), deltaRef(date));

  185.             return new AzimuthalGradientCoefficients(interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.GN_H) * 0.00001),
  186.                                                      interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.GE_H) * 0.00001),
  187.                                                      interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.GN_W) * 0.00001),
  188.                                                      interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.GE_W) * 0.00001));
  189.         } else {
  190.             return null;
  191.         }

  193.     }

  195.     /** {@inheritDoc} */
  196.     @Override
  197.     public <T extends CalculusFieldElement<T>> FieldViennaACoefficients<T> getA(final FieldGeodeticPoint<T> location,
  198.                                                                                 final FieldAbsoluteDate<T> date) {

  200.         // set up interpolation parameters
  201.         final FieldCellInterpolator<T> interpolator =
  202.             grid.getInterpolator(location.getLatitude(), location.getLongitude(), location.getAltitude(),
  203.                                  deltaRef(date));

  205.         // ah and aw coefficients
  206.         return new FieldViennaACoefficients<>(interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.AH)).multiply(0.001),
  207.                                               interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.AW)).multiply(0.001));

  209.     }

  211.     /** {@inheritDoc} */
  212.     @Override
  213.     public <T extends CalculusFieldElement<T>> FieldPressureTemperatureHumidity<T> getWeatherParameters(final FieldGeodeticPoint<T> location,
  214.                                                                                                         final FieldAbsoluteDate<T> date) {

  216.         // set up interpolation parameters
  217.         final FieldCellInterpolator<T> interpolator =
  218.             grid.getInterpolator(location.getLatitude(), location.getLongitude(), location.getAltitude(),
  219.                                  deltaRef(date));

  221.         // Temperature [K]
  222.         final T temperature = interpolator.fieldInterpolate(FieldEvaluatedGridEntry::getTemperature);

  224.         // Pressure [hPa]
  225.         final T pressure = interpolator.fieldInterpolate(FieldEvaluatedGridEntry::getPressure);

  227.         // water vapor decrease factor
  228.         final T lambda = grid.hasModels(SeasonalModelType.LAMBDA) ?
  229.                          interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.LAMBDA)) :
  230.                          date.getField().getZero().newInstance(Double.NaN);

  232.         // Water vapor pressure [hPa]
  233.         final T el;
  234.         if (lambda.isNaN()) {
  235.             // GPT2
  236.             final T qv = interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.QV)).multiply(0.001);
  237.             el = qv.multiply(pressure).divide(qv.multiply(0.378).add(0.622));
  238.         } else {
  239.             // GPT3
  240.             el = interpolator.fieldInterpolate(FieldEvaluatedGridEntry::getWaterVaporPressure);
  241.         }

  243.         // mean temperature weighted with water vapor pressure
  244.         final T tm = grid.hasModels(SeasonalModelType.TM) ?
  245.                      interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.TM)) :
  246.                      date.getField().getZero().newInstance(Double.NaN);

  248.         return new FieldPressureTemperatureHumidity<>(location.getAltitude(),
  249.                                                       TroposphericModelUtils.HECTO_PASCAL.toSI(pressure),
  250.                                                       temperature,
  251.                                                       TroposphericModelUtils.HECTO_PASCAL.toSI(el),
  252.                                                       tm,
  253.                                                       lambda);

  255.     }

  257.     /** {@inheritDoc} */
  258.     @Override
  259.     public <T extends CalculusFieldElement<T>> FieldAzimuthalGradientCoefficients<T> getGradientCoefficients(final FieldGeodeticPoint<T> location,
  260.                                                                                                              final FieldAbsoluteDate<T> date) {

  262.         if (grid.hasModels(SeasonalModelType.GN_H, SeasonalModelType.GE_H, SeasonalModelType.GN_W, SeasonalModelType.GE_W)) {
  263.             // set up interpolation parameters
  264.             final FieldCellInterpolator<T> interpolator =
  265.                 grid.getInterpolator(location.getLatitude(), location.getLongitude(), location.getAltitude(),
  266.                                      deltaRef(date));

  268.             return new FieldAzimuthalGradientCoefficients<>(interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.GN_H)).multiply(0.00001),
  269.                                                             interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.GE_H)).multiply(0.00001),
  270.                                                             interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.GN_W)).multiply(0.00001),
  271.                                                             interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.GE_W)).multiply(0.00001));
  272.         } else {
  273.             return null;
  274.         }

  276.     }

  278. }
Created with JaCoCo 0.8.12.202403310830