ViennaModelCoefficientsLoader.java
/* Copyright 2002-2019 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (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;
import java.io.BufferedReader;
import java.io.IOException;
import java.io.InputStream;
import java.io.InputStreamReader;
import java.text.ParseException;
import java.util.ArrayList;
import org.hipparchus.analysis.interpolation.BilinearInterpolatingFunction;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathUtils;
import org.orekit.data.DataLoader;
import org.orekit.data.DataProvidersManager;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitMessages;
import org.orekit.time.DateTimeComponents;
/** Loads Vienna tropospheric coefficients a given input stream.
* A stream contains, for a given day and a given hour, the hydrostatic and wet zenith delays
* and the ah and aw coefficients used for the computation of the mapping function.
* The coefficients are given with a time interval of 6 hours.
* <p>
* A bilinear interpolation is performed the case of the user initialize the latitude and the
* longitude with values that are not contained in the stream.
* </p>
* <p>
* The coefficients are obtained from <a href="http://vmf.geo.tuwien.ac.at/trop_products/GRID/">Vienna Mapping Functions Open Access Data</a>.
* Find more on the files at the <a href="http://vmf.geo.tuwien.ac.at/readme.txt">VMF Model Documentation</a>.
* <p>
* The files have to be extracted to UTF-8 text files before being read by this loader.
* <p>
* After extraction, it is assumed they are named VMFG_YYYYMMDD.Hhh for {@link ViennaOneModel} and VMF3_YYYYMMDD.Hhh {@link ViennaThreeModel}.
* Where YYYY is the 4-digits year, MM the month, DD the day and hh the 2-digits hour.
*
* <p>
* The format is always the same, with and example shown below for VMF1 model.
* <p>
* Example:
* </p>
* <pre>
* ! Version: 1.0
* ! Source: J. Boehm, TU Vienna (created: 2018-11-20)
* ! Data_types: VMF1 (lat lon ah aw zhd zwd)
* ! Epoch: 2018 11 19 18 00 0.0
* ! Scale_factor: 1.e+00
* ! Range/resolution: -90 90 0 360 2 2.5
* ! Comment: http://vmf.geo.tuwien.ac.at/trop_products/GRID/2.5x2/VMF1/VMF1_OP/
* 90.0 0.0 0.00116059 0.00055318 2.3043 0.0096
* 90.0 2.5 0.00116059 0.00055318 2.3043 0.0096
* 90.0 5.0 0.00116059 0.00055318 2.3043 0.0096
* 90.0 7.5 0.00116059 0.00055318 2.3043 0.0096
* 90.0 10.0 0.00116059 0.00055318 2.3043 0.0096
* 90.0 12.5 0.00116059 0.00055318 2.3043 0.0096
* 90.0 15.0 0.00116059 0.00055318 2.3043 0.0096
* 90.0 17.5 0.00116059 0.00055318 2.3043 0.0096
* 90.0 20.0 0.00116059 0.00055318 2.3043 0.0096
* 90.0 22.5 0.00116059 0.00055318 2.3043 0.0096
* 90.0 25.0 0.00116059 0.00055318 2.3043 0.0096
* 90.0 27.5 0.00116059 0.00055318 2.3043 0.0096
* </pre>
* @author Bryan Cazabonne
*/
public class ViennaModelCoefficientsLoader implements DataLoader {
/** Default supported files name pattern. */
public static final String DEFAULT_SUPPORTED_NAMES = "VMF*_\\\\*\\*\\.*H$";
/** Regular expression for supported file name. */
private String supportedNames;
/** The hydrostatic and wet a coefficients loaded. */
private double[] coefficientsA;
/** The hydrostatic and wet zenith delays loaded. */
private double[] zenithDelay;
/** Geodetic site latitude, radians.*/
private double latitude;
/** Geodetic site longitude, radians.*/
private double longitude;
/** Vienna tropospheric model type.*/
private ViennaModelType type;
/** Constructor with supported names given by user.
* @param supportedNames Supported names
* @param latitude geodetic latitude of the station, in radians
* @param longitude geodetic latitude of the station, in radians
* @param type the type of Vienna tropospheric model (one or three)
*/
public ViennaModelCoefficientsLoader(final String supportedNames, final double latitude,
final double longitude, final ViennaModelType type) {
this.coefficientsA = null;
this.zenithDelay = null;
this.supportedNames = supportedNames;
this.type = type;
this.latitude = latitude;
// Normalize longitude between 0 and 2π
this.longitude = MathUtils.normalizeAngle(longitude, FastMath.PI);
}
/** Constructor with default supported names.
* @param latitude geodetic latitude of the station, in radians
* @param longitude geodetic latitude of the station, in radians
* @param type the type of Vienna tropospheric model (one or three)
*/
public ViennaModelCoefficientsLoader(final double latitude, final double longitude,
final ViennaModelType type) {
this(DEFAULT_SUPPORTED_NAMES, latitude, longitude, type);
}
/** Returns the a coefficients array.
* <ul>
* <li>double[0] = a<sub>h</sub>
* <li>double[1] = a<sub>w</sub>
* </ul>
* @return the a coefficients array
*/
public double[] getA() {
return coefficientsA.clone();
}
/** Returns the zenith delay array.
* <ul>
* <li>double[0] = D<sub>hz</sub> → zenith hydrostatic delay
* <li>double[1] = D<sub>wz</sub> → zenith wet delay
* </ul>
* @return the zenith delay array
*/
public double[] getZenithDelay() {
return zenithDelay.clone();
}
/** Returns the supported names of the loader.
* @return the supported names
*/
public String getSupportedNames() {
return supportedNames;
}
/** Load the data using supported names .
*/
public void loadViennaCoefficients() {
DataProvidersManager.getInstance().feed(supportedNames, this);
// Throw an exception if ah, ah, zh or zw were not loaded properly
if (coefficientsA == null || zenithDelay == null) {
throw new OrekitException(OrekitMessages.VIENNA_ACOEF_OR_ZENITH_DELAY_NOT_LOADED, supportedNames);
}
}
/** Load the data for a given day.
* @param dateTimeComponents date and time component.
* @throws OrekitException if the coefficients could not be loaded
*/
public void loadViennaCoefficients(final DateTimeComponents dateTimeComponents) {
// The files are named VMFG_YYYYMMDD.Hhh for Vienna-1 model and VMF3_YYYYMMDD.Hhh for Vienna-3 model.
// Where YYYY is the 4-digits year, MM the month, DD the day of the month and hh the 2-digits hour.
// Coefficients are only available for hh = 00 or 06 or 12 or 18.
final int year = dateTimeComponents.getDate().getYear();
final int month = dateTimeComponents.getDate().getMonth();
final int day = dateTimeComponents.getDate().getDay();
final int hour = dateTimeComponents.getTime().getHour();
// Correct month format is with 2-digits.
final String monthString;
if (month < 10) {
monthString = "0" + String.valueOf(month);
} else {
monthString = String.valueOf(month);
}
// Correct day format is with 2-digits.
final String dayString;
if (day < 10) {
dayString = "0" + String.valueOf(day);
} else {
dayString = String.valueOf(day);
}
// Correct hour format is with 2-digits.
final String hourString;
if (hour < 10) {
hourString = "0" + String.valueOf(hour);
} else {
hourString = String.valueOf(hour);
}
// Name of the file is different between VMF1 and VMF3.
// For VMF1 it starts with "VMFG" whereas with VMF3 it starts with "VMF3"
switch (type) {
case VIENNA_ONE:
this.supportedNames = String.format("VMFG_%04d%2s%2s.H%2s", year, monthString, dayString, hourString);
break;
case VIENNA_THREE:
this.supportedNames = String.format("VMF3_%04d%2s%2s.H%2s", year, monthString, dayString, hourString);
break;
default:
break;
}
try {
this.loadViennaCoefficients();
} catch (OrekitException oe) {
throw new OrekitException(oe,
OrekitMessages.VIENNA_ACOEF_OR_ZENITH_DELAY_NOT_AVAILABLE_FOR_DATE,
dateTimeComponents.toString());
}
}
@Override
public boolean stillAcceptsData() {
return true;
}
@Override
public void loadData(final InputStream input, final String name)
throws IOException, ParseException {
// Open stream and parse data
final BufferedReader br = new BufferedReader(new InputStreamReader(input, "UTF-8"));
int lineNumber = 0;
final String splitter = "\\s+";
// Initialize Lists
final ArrayList<Double> latitudes = new ArrayList<>();
final ArrayList<Double> longitudes = new ArrayList<>();
final ArrayList<Double> ah = new ArrayList<>();
final ArrayList<Double> aw = new ArrayList<>();
final ArrayList<Double> zhd = new ArrayList<>();
final ArrayList<Double> zwd = new ArrayList<>();
for (String line = br.readLine(); line != null; line = br.readLine()) {
++lineNumber;
line = line.trim();
try {
// Fill latitudes and longitudes lists
if (line.length() > 0 && line.startsWith("! Range/resolution:")) {
final String[] range_line = line.split(splitter);
// Latitudes list
for (double lat = Double.valueOf(range_line[2]); lat <= Double.valueOf(range_line[3]); lat = lat + Double.valueOf(range_line[6])) {
latitudes.add(FastMath.toRadians(lat));
}
// Longitude list
for (double lon = Double.valueOf(range_line[4]); lon <= Double.valueOf(range_line[5]); lon = lon + Double.valueOf(range_line[7])) {
longitudes.add(FastMath.toRadians(lon));
// For VFM1 files, header specify that longitudes end at 360°
// In reality they end at 357.5°. That is why we stop the loop when the longitude
// reaches 357.5°.
if (type == ViennaModelType.VIENNA_ONE && lon >= 357.5) {
break;
}
}
}
// Fill ah, aw, zhd and zwd lists
if (line.length() > 0 && !line.startsWith("!")) {
final String[] values_line = line.split(splitter);
ah.add(Double.valueOf(values_line[2]));
aw.add(Double.valueOf(values_line[3]));
zhd.add(Double.valueOf(values_line[4]));
zwd.add(Double.valueOf(values_line[5]));
}
} catch (NumberFormatException nfe) {
throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
lineNumber, name, line);
}
}
// Close the stream after reading
input.close();
final int dimLat = latitudes.size();
final int dimLon = longitudes.size();
// Change Lists to Arrays
final double[] xVal = new double[dimLat];
for (int i = 0; i < dimLat; i++) {
xVal[i] = latitudes.get(i);
}
final double[] yVal = new double[dimLon];
for (int j = 0; j < dimLon; j++) {
yVal[j] = longitudes.get(j);
}
final double[][] fvalAH = new double[dimLat][dimLon];
final double[][] fvalAW = new double[dimLat][dimLon];
final double[][] fvalZH = new double[dimLat][dimLon];
final double[][] fvalZW = new double[dimLat][dimLon];
int index = dimLon * dimLat;
for (int x = 0; x < dimLat; x++) {
for (int y = dimLon - 1; y >= 0; y--) {
index = index - 1;
fvalAH[x][y] = ah.get(index);
fvalAW[x][y] = aw.get(index);
fvalZH[x][y] = zhd.get(index);
fvalZW[x][y] = zwd.get(index);
}
}
// Build Bilinear Interpolation Functions
final BilinearInterpolatingFunction functionAH = new BilinearInterpolatingFunction(xVal, yVal, fvalAH);
final BilinearInterpolatingFunction functionAW = new BilinearInterpolatingFunction(xVal, yVal, fvalAW);
final BilinearInterpolatingFunction functionZH = new BilinearInterpolatingFunction(xVal, yVal, fvalZH);
final BilinearInterpolatingFunction functionZW = new BilinearInterpolatingFunction(xVal, yVal, fvalZW);
coefficientsA = new double[2];
zenithDelay = new double[2];
// Get the values for the given latitude and longitude
coefficientsA[0] = functionAH.value(latitude, longitude);
coefficientsA[1] = functionAW.value(latitude, longitude);
zenithDelay[0] = functionZH.value(latitude, longitude);
zenithDelay[1] = functionZW.value(latitude, longitude);
// Check that ah, aw, zh and zw were found
if (coefficientsA == null || zenithDelay == null) {
throw new OrekitException(OrekitMessages.NO_VIENNA_ACOEF_OR_ZENITH_DELAY_IN_FILE, name);
}
}
}