AbstractGlobalPressureTemperature.java
/* Copyright 2002-2024 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package org.orekit.models.earth.weather;
import java.io.BufferedInputStream;
import java.io.IOException;
import java.io.InputStream;
import org.hipparchus.CalculusFieldElement;
import org.orekit.bodies.FieldGeodeticPoint;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.data.DataSource;
import org.orekit.models.earth.troposphere.AzimuthalGradientCoefficients;
import org.orekit.models.earth.troposphere.AzimuthalGradientProvider;
import org.orekit.models.earth.troposphere.FieldAzimuthalGradientCoefficients;
import org.orekit.models.earth.troposphere.FieldViennaACoefficients;
import org.orekit.models.earth.troposphere.TroposphericModelUtils;
import org.orekit.models.earth.troposphere.ViennaACoefficients;
import org.orekit.models.earth.troposphere.ViennaAProvider;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.FieldAbsoluteDate;
/** Base class for Global Pressure and Temperature 2, 2w and 3 models.
* These models are empirical models that provide the temperature, the pressure and the water vapor pressure
* of a site depending its latitude and longitude. These models also {@link ViennaACoefficients provide}
* the a<sub>h</sub> and a<sub>w</sub> coefficients for Vienna models.
* <p>
* The requisite coefficients for the computation of the weather parameters are provided by the
* Department of Geodesy and Geoinformation of the Vienna University. They are based on an
* external grid file like "gpt2_1.grd" (1° x 1°), "gpt2_5.grd" (5° x 5°), "gpt2_1w.grd" (1° x 1°),
* "gpt2_5w.grd" (5° x 5°), "gpt3_1.grd" (1° x 1°), or "gpt3_5.grd" (5° x 5°) available at:
* <a href="https://vmf.geo.tuwien.ac.at/codes/"> link</a>
* </p>
* <p>
* A bilinear interpolation is performed in order to obtained the correct values of the weather parameters.
* </p>
* <p>
* The format is always the same, with and example shown below for the pressure and the temperature.
* The "GPT2w" model (w stands for wet) also provide humidity parameters and the "GPT3" model also
* provides horizontal gradient, so the number of columns vary depending on the model.
* <p>
* Example:
* </p>
* <pre>
* % lat lon p:a0 A1 B1 A2 B2 T:a0 A1 B1 A2 B2
* 87.5 2.5 101421 21 409 -217 -122 259.2 -13.2 -6.1 2.6 0.3
* 87.5 7.5 101416 21 411 -213 -120 259.3 -13.1 -6.1 2.6 0.3
* 87.5 12.5 101411 22 413 -209 -118 259.3 -13.1 -6.1 2.6 0.3
* 87.5 17.5 101407 23 415 -205 -116 259.4 -13.0 -6.1 2.6 0.3
* ...
* </pre>
*
* @see "K. Lagler, M. Schindelegger, J. Böhm, H. Krasna, T. Nilsson (2013),
* GPT2: empirical slant delay model for radio space geodetic techniques. Geophys
* Res Lett 40(6):1069–1073. doi:10.1002/grl.50288"
*
* @author Bryan Cazabonne
* @author Luc Maisonobe
* @since 12.1
*/
public abstract class AbstractGlobalPressureTemperature
implements ViennaAProvider, AzimuthalGradientProvider, PressureTemperatureHumidityProvider {
/** Loaded grid. */
private final Grid grid;
/**
* Constructor with source of GPTn auxiliary data given by user.
*
* @param source grid data source
* @param expected expected seasonal models types
* @exception IOException if grid data cannot be read
*/
protected AbstractGlobalPressureTemperature(final DataSource source, final SeasonalModelType... expected)
throws IOException {
// load the grid data
try (InputStream is = source.getOpener().openStreamOnce();
BufferedInputStream bis = new BufferedInputStream(is)) {
final GptNParser parser = new GptNParser(expected);
parser.loadData(bis, source.getName());
grid = parser.getGrid();
}
}
/** Get duration since reference date.
* @param date date
* @return duration since reference date
* @since 13.0
*/
protected abstract double deltaRef(AbsoluteDate date);
/** Get duration since reference date.
* @param <T> type of the field elements
* @param date date
* @return duration since reference date
* @since 13.0
*/
protected abstract <T extends CalculusFieldElement<T>> T deltaRef(FieldAbsoluteDate<T> date);
/** {@inheritDoc} */
@Override
public ViennaACoefficients getA(final GeodeticPoint location, final AbsoluteDate date) {
// set up interpolation parameters
final CellInterpolator interpolator =
grid.getInterpolator(location.getLatitude(), location.getLongitude(), location.getAltitude(),
deltaRef(date));
// ah and aw coefficients
return new ViennaACoefficients(interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.AH)) * 0.001,
interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.AW)) * 0.001);
}
/** {@inheritDoc} */
@Override
public PressureTemperatureHumidity getWeatherParameters(final GeodeticPoint location, final AbsoluteDate date) {
// set up interpolation parameters
final CellInterpolator interpolator =
grid.getInterpolator(location.getLatitude(), location.getLongitude(), location.getAltitude(),
deltaRef(date));
// Temperature [K]
final double temperature = interpolator.interpolate(EvaluatedGridEntry::getTemperature);
// Pressure [hPa]
final double pressure = interpolator.interpolate(EvaluatedGridEntry::getPressure);
// water vapor decrease factor
final double lambda = grid.hasModels(SeasonalModelType.LAMBDA) ?
interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.LAMBDA)) :
Double.NaN;
// Water vapor pressure [hPa]
final double el;
if (Double.isNaN(lambda)) {
// GPT2
final double qv = interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.QV)) * 0.001;
el = qv * pressure / (0.622 + 0.378 * qv);
} else {
// GPT2w and GPT3
el = interpolator.interpolate(EvaluatedGridEntry::getWaterVaporPressure);
}
// mean temperature weighted with water vapor pressure
final double tm = grid.hasModels(SeasonalModelType.TM) ?
interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.TM)) :
Double.NaN;
return new PressureTemperatureHumidity(location.getAltitude(),
TroposphericModelUtils.HECTO_PASCAL.toSI(pressure),
temperature,
TroposphericModelUtils.HECTO_PASCAL.toSI(el),
tm,
lambda);
}
/** {@inheritDoc} */
@Override
public AzimuthalGradientCoefficients getGradientCoefficients(final GeodeticPoint location, final AbsoluteDate date) {
if (grid.hasModels(SeasonalModelType.GN_H, SeasonalModelType.GE_H, SeasonalModelType.GN_W, SeasonalModelType.GE_W)) {
// set up interpolation parameters
final CellInterpolator interpolator = grid.getInterpolator(location.getLatitude(), location.getLongitude(),
location.getAltitude(), deltaRef(date));
return new AzimuthalGradientCoefficients(interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.GN_H) * 0.00001),
interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.GE_H) * 0.00001),
interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.GN_W) * 0.00001),
interpolator.interpolate(e -> e.getEvaluatedModel(SeasonalModelType.GE_W) * 0.00001));
} else {
return null;
}
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> FieldViennaACoefficients<T> getA(final FieldGeodeticPoint<T> location,
final FieldAbsoluteDate<T> date) {
// set up interpolation parameters
final FieldCellInterpolator<T> interpolator =
grid.getInterpolator(location.getLatitude(), location.getLongitude(), location.getAltitude(),
deltaRef(date));
// ah and aw coefficients
return new FieldViennaACoefficients<>(interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.AH)).multiply(0.001),
interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.AW)).multiply(0.001));
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> FieldPressureTemperatureHumidity<T> getWeatherParameters(final FieldGeodeticPoint<T> location,
final FieldAbsoluteDate<T> date) {
// set up interpolation parameters
final FieldCellInterpolator<T> interpolator =
grid.getInterpolator(location.getLatitude(), location.getLongitude(), location.getAltitude(),
deltaRef(date));
// Temperature [K]
final T temperature = interpolator.fieldInterpolate(FieldEvaluatedGridEntry::getTemperature);
// Pressure [hPa]
final T pressure = interpolator.fieldInterpolate(FieldEvaluatedGridEntry::getPressure);
// water vapor decrease factor
final T lambda = grid.hasModels(SeasonalModelType.LAMBDA) ?
interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.LAMBDA)) :
date.getField().getZero().newInstance(Double.NaN);
// Water vapor pressure [hPa]
final T el;
if (lambda.isNaN()) {
// GPT2
final T qv = interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.QV)).multiply(0.001);
el = qv.multiply(pressure).divide(qv.multiply(0.378).add(0.622));
} else {
// GPT3
el = interpolator.fieldInterpolate(FieldEvaluatedGridEntry::getWaterVaporPressure);
}
// mean temperature weighted with water vapor pressure
final T tm = grid.hasModels(SeasonalModelType.TM) ?
interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.TM)) :
date.getField().getZero().newInstance(Double.NaN);
return new FieldPressureTemperatureHumidity<>(location.getAltitude(),
TroposphericModelUtils.HECTO_PASCAL.toSI(pressure),
temperature,
TroposphericModelUtils.HECTO_PASCAL.toSI(el),
tm,
lambda);
}
/** {@inheritDoc} */
@Override
public <T extends CalculusFieldElement<T>> FieldAzimuthalGradientCoefficients<T> getGradientCoefficients(final FieldGeodeticPoint<T> location,
final FieldAbsoluteDate<T> date) {
if (grid.hasModels(SeasonalModelType.GN_H, SeasonalModelType.GE_H, SeasonalModelType.GN_W, SeasonalModelType.GE_W)) {
// set up interpolation parameters
final FieldCellInterpolator<T> interpolator =
grid.getInterpolator(location.getLatitude(), location.getLongitude(), location.getAltitude(),
deltaRef(date));
return new FieldAzimuthalGradientCoefficients<>(interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.GN_H)).multiply(0.00001),
interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.GE_H)).multiply(0.00001),
interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.GN_W)).multiply(0.00001),
interpolator.fieldInterpolate(e -> e.getEvaluatedModel(SeasonalModelType.GE_W)).multiply(0.00001));
} else {
return null;
}
}
}