ViennaModelCoefficientsLoader.java

  1. /* Copyright 2002-2020 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.troposphere;

  18. import java.io.BufferedReader;
  19. import java.io.IOException;
  20. import java.io.InputStream;
  21. import java.io.InputStreamReader;
  22. import java.nio.charset.StandardCharsets;
  23. import java.text.ParseException;
  24. import java.util.ArrayList;
  25. import java.util.regex.Pattern;

  26. import org.hipparchus.analysis.interpolation.BilinearInterpolatingFunction;
  27. import org.hipparchus.util.FastMath;
  28. import org.hipparchus.util.MathUtils;
  29. import org.orekit.annotation.DefaultDataContext;
  30. import org.orekit.data.AbstractSelfFeedingLoader;
  31. import org.orekit.data.DataContext;
  32. import org.orekit.data.DataLoader;
  33. import org.orekit.data.DataProvidersManager;
  34. import org.orekit.errors.OrekitException;
  35. import org.orekit.errors.OrekitMessages;
  36. import org.orekit.time.DateTimeComponents;

  37. /** Loads Vienna tropospheric coefficients a given input stream.
  38.  * A stream contains, for a given day and a given hour, the hydrostatic and wet zenith delays
  39.  * and the ah and aw coefficients used for the computation of the mapping function.
  40.  * The coefficients are given with a time interval of 6 hours.
  41.  * <p>
  42.  * A bilinear interpolation is performed the case of the user initialize the latitude and the
  43.  * longitude with values that are not contained in the stream.
  44.  * </p>
  45.  * <p>
  46.  * The coefficients are obtained from <a href="http://vmf.geo.tuwien.ac.at/trop_products/GRID/">Vienna Mapping Functions Open Access Data</a>.
  47.  * Find more on the files at the <a href="http://vmf.geo.tuwien.ac.at/readme.txt">VMF Model Documentation</a>.
  48.  * <p>
  49.  * The files have to be extracted to UTF-8 text files before being read by this loader.
  50.  * <p>
  51.  * After extraction, it is assumed they are named VMFG_YYYYMMDD.Hhh for {@link ViennaOneModel} and VMF3_YYYYMMDD.Hhh {@link ViennaThreeModel}.
  52.  * Where YYYY is the 4-digits year, MM the month, DD the day and hh the 2-digits hour.
  53.  *
  54.  * <p>
  55.  * The format is always the same, with and example shown below for VMF1 model.
  56.  * <p>
  57.  * Example:
  58.  * </p>
  59.  * <pre>
  60.  * ! Version:            1.0
  61.  * ! Source:             J. Boehm, TU Vienna (created: 2018-11-20)
  62.  * ! Data_types:         VMF1 (lat lon ah aw zhd zwd)
  63.  * ! Epoch:              2018 11 19 18 00  0.0
  64.  * ! Scale_factor:       1.e+00
  65.  * ! Range/resolution:   -90 90 0 360 2 2.5
  66.  * ! Comment:            http://vmf.geo.tuwien.ac.at/trop_products/GRID/2.5x2/VMF1/VMF1_OP/
  67.  *  90.0   0.0 0.00116059  0.00055318  2.3043  0.0096
  68.  *  90.0   2.5 0.00116059  0.00055318  2.3043  0.0096
  69.  *  90.0   5.0 0.00116059  0.00055318  2.3043  0.0096
  70.  *  90.0   7.5 0.00116059  0.00055318  2.3043  0.0096
  71.  *  90.0  10.0 0.00116059  0.00055318  2.3043  0.0096
  72.  *  90.0  12.5 0.00116059  0.00055318  2.3043  0.0096
  73.  *  90.0  15.0 0.00116059  0.00055318  2.3043  0.0096
  74.  *  90.0  17.5 0.00116059  0.00055318  2.3043  0.0096
  75.  *  90.0  20.0 0.00116059  0.00055318  2.3043  0.0096
  76.  *  90.0  22.5 0.00116059  0.00055318  2.3043  0.0096
  77.  *  90.0  25.0 0.00116059  0.00055318  2.3043  0.0096
  78.  *  90.0  27.5 0.00116059  0.00055318  2.3043  0.0096
  79.  * </pre>
  80.  *
  81.  * <p>It is not safe for multiple threads to share a single instance of this class.
  82.  *
  83.  * @author Bryan Cazabonne
  84.  */
  85. public class ViennaModelCoefficientsLoader extends AbstractSelfFeedingLoader
  86.         implements DataLoader {

  87.     /** Default supported files name pattern. */
  88.     public static final String DEFAULT_SUPPORTED_NAMES = "VMF*_\\\\*\\*\\.*H$";

  89.     /** Pattern for delimiting regular expressions. */
  90.     private static final Pattern SEPARATOR = Pattern.compile("\\s+");

  91.     /** The hydrostatic and wet a coefficients loaded. */
  92.     private double[] coefficientsA;

  93.     /** The hydrostatic and wet zenith delays loaded. */
  94.     private double[] zenithDelay;

  95.     /** Geodetic site latitude, radians.*/
  96.     private double latitude;

  97.     /** Geodetic site longitude, radians.*/
  98.     private double longitude;

  99.     /** Vienna tropospheric model type.*/
  100.     private ViennaModelType type;

  101.     /** Constructor with supported names given by user. This constructor uses the
  102.      * {@link DataContext#getDefault() default data context}.
  103.      *
  104.      * @param supportedNames Supported names
  105.      * @param latitude geodetic latitude of the station, in radians
  106.      * @param longitude geodetic latitude of the station, in radians
  107.      * @param type the type of Vienna tropospheric model (one or three)
  108.      * @see #ViennaModelCoefficientsLoader(String, double, double, ViennaModelType, DataProvidersManager)
  109.      */
  110.     @DefaultDataContext
  111.     public ViennaModelCoefficientsLoader(final String supportedNames, final double latitude,
  112.                                          final double longitude, final ViennaModelType type) {
  113.         this(supportedNames, latitude, longitude, type, DataContext.getDefault().getDataProvidersManager());
  114.     }

  115.     /**
  116.      * Constructor with supported names and source of mapping function files given by the
  117.      * user.
  118.      *
  119.      * @param supportedNames Supported names
  120.      * @param latitude       geodetic latitude of the station, in radians
  121.      * @param longitude      geodetic latitude of the station, in radians
  122.      * @param type           the type of Vienna tropospheric model (one or three)
  123.      * @param dataProvidersManager provides access to auxiliary files.
  124.      * @since 10.1
  125.      */
  126.     public ViennaModelCoefficientsLoader(final String supportedNames,
  127.                                          final double latitude,
  128.                                          final double longitude,
  129.                                          final ViennaModelType type,
  130.                                          final DataProvidersManager dataProvidersManager) {
  131.         super(supportedNames, dataProvidersManager);
  132.         this.coefficientsA  = null;
  133.         this.zenithDelay    = null;
  134.         this.type           = type;
  135.         this.latitude       = latitude;

  136.         // Normalize longitude between 0 and 2π
  137.         this.longitude = MathUtils.normalizeAngle(longitude, FastMath.PI);
  138.     }

  139.     /** Constructor with default supported names. This constructor uses the
  140.      * {@link DataContext#getDefault() default data context}.
  141.      *
  142.      * @param latitude geodetic latitude of the station, in radians
  143.      * @param longitude geodetic latitude of the station, in radians
  144.      * @param type the type of Vienna tropospheric model (one or three)
  145.      * @see #ViennaModelCoefficientsLoader(String, double, double, ViennaModelType, DataProvidersManager)
  146.      */
  147.     @DefaultDataContext
  148.     public ViennaModelCoefficientsLoader(final double latitude, final double longitude,
  149.                                          final ViennaModelType type) {
  150.         this(DEFAULT_SUPPORTED_NAMES, latitude, longitude, type);
  151.     }

  152.     /** Returns the a coefficients array.
  153.      * <ul>
  154.      * <li>double[0] = a<sub>h</sub>
  155.      * <li>double[1] = a<sub>w</sub>
  156.      * </ul>
  157.      * @return the a coefficients array
  158.      */
  159.     public double[] getA() {
  160.         return coefficientsA.clone();
  161.     }

  162.     /** Returns the zenith delay array.
  163.      * <ul>
  164.      * <li>double[0] = D<sub>hz</sub> → zenith hydrostatic delay
  165.      * <li>double[1] = D<sub>wz</sub> → zenith wet delay
  166.      * </ul>
  167.      * @return the zenith delay array
  168.      */
  169.     public double[] getZenithDelay() {
  170.         return zenithDelay.clone();
  171.     }

  172.     @Override
  173.     public String getSupportedNames() {
  174.         return super.getSupportedNames();
  175.     }

  176.     /** Load the data using supported names .
  177.      */
  178.     public void loadViennaCoefficients() {
  179.         feed(this);

  180.         // Throw an exception if ah, ah, zh or zw were not loaded properly
  181.         if (coefficientsA == null || zenithDelay == null) {
  182.             throw new OrekitException(OrekitMessages.VIENNA_ACOEF_OR_ZENITH_DELAY_NOT_LOADED,
  183.                     getSupportedNames());
  184.         }
  185.     }

  186.     /** Load the data for a given day.
  187.      * @param dateTimeComponents date and time component.
  188.      */
  189.     public void loadViennaCoefficients(final DateTimeComponents dateTimeComponents) {

  190.         // The files are named VMFG_YYYYMMDD.Hhh for Vienna-1 model and VMF3_YYYYMMDD.Hhh for Vienna-3 model.
  191.         // Where YYYY is the 4-digits year, MM the month, DD the day of the month and hh the 2-digits hour.
  192.         // Coefficients are only available for hh = 00 or 06 or 12 or 18.
  193.         final int    year        = dateTimeComponents.getDate().getYear();
  194.         final int    month       = dateTimeComponents.getDate().getMonth();
  195.         final int    day         = dateTimeComponents.getDate().getDay();
  196.         final int    hour        = dateTimeComponents.getTime().getHour();

  197.         // Correct month format is with 2-digits.
  198.         final String monthString;
  199.         if (month < 10) {
  200.             monthString = "0" + month;
  201.         } else {
  202.             monthString = String.valueOf(month);
  203.         }

  204.         // Correct day format is with 2-digits.
  205.         final String dayString;
  206.         if (day < 10) {
  207.             dayString = "0" + day;
  208.         } else {
  209.             dayString = String.valueOf(day);
  210.         }

  211.         // Correct hour format is with 2-digits.
  212.         final String hourString;
  213.         if (hour < 10) {
  214.             hourString = "0" + hour;
  215.         } else {
  216.             hourString = String.valueOf(hour);
  217.         }

  218.         // Name of the file is different between VMF1 and VMF3.
  219.         // For VMF1 it starts with "VMFG" whereas with VMF3 it starts with "VMF3"
  220.         switch (type) {
  221.             case VIENNA_ONE:
  222.                 setSupportedNames(String.format("VMFG_%04d%2s%2s.H%2s", year, monthString, dayString, hourString));
  223.                 break;
  224.             case VIENNA_THREE:
  225.                 setSupportedNames(String.format("VMF3_%04d%2s%2s.H%2s", year, monthString, dayString, hourString));
  226.                 break;
  227.             default:
  228.                 break;
  229.         }

  230.         try {
  231.             this.loadViennaCoefficients();
  232.         } catch (OrekitException oe) {
  233.             throw new OrekitException(oe,
  234.                                       OrekitMessages.VIENNA_ACOEF_OR_ZENITH_DELAY_NOT_AVAILABLE_FOR_DATE,
  235.                                       dateTimeComponents.toString());
  236.         }
  237.     }

  238.     @Override
  239.     public boolean stillAcceptsData() {
  240.         return true;
  241.     }

  242.     @Override
  243.     public void loadData(final InputStream input, final String name)
  244.         throws IOException, ParseException {

  245.         int lineNumber = 0;
  246.         String line = null;

  247.         // Initialize Lists
  248.         final ArrayList<Double> latitudes  = new ArrayList<>();
  249.         final ArrayList<Double> longitudes = new ArrayList<>();
  250.         final ArrayList<Double> ah         = new ArrayList<>();
  251.         final ArrayList<Double> aw         = new ArrayList<>();
  252.         final ArrayList<Double> zhd        = new ArrayList<>();
  253.         final ArrayList<Double> zwd        = new ArrayList<>();

  254.         // Open stream and parse data
  255.         try (BufferedReader br = new BufferedReader(new InputStreamReader(input, StandardCharsets.UTF_8))) {

  256.             for (line = br.readLine(); line != null; line = br.readLine()) {
  257.                 ++lineNumber;
  258.                 line = line.trim();

  259.                 // Fill latitudes and longitudes lists
  260.                 if (line.length() > 0 && line.startsWith("! Range/resolution:")) {
  261.                     final String[] range_line = SEPARATOR.split(line);

  262.                     // Latitudes list
  263.                     for (double lat = Double.parseDouble(range_line[2]); lat <= Double.parseDouble(range_line[3]); lat = lat + Double.parseDouble(range_line[6])) {
  264.                         latitudes.add(FastMath.toRadians(lat));
  265.                     }

  266.                     // Longitude list
  267.                     for (double lon = Double.parseDouble(range_line[4]); lon <= Double.parseDouble(range_line[5]); lon = lon + Double.parseDouble(range_line[7])) {
  268.                         longitudes.add(FastMath.toRadians(lon));
  269.                         // For VFM1 files, header specify that longitudes end at 360°
  270.                         // In reality they end at 357.5°. That is why we stop the loop when the longitude
  271.                         // reaches 357.5°.
  272.                         if (type == ViennaModelType.VIENNA_ONE && lon >= 357.5) {
  273.                             break;
  274.                         }
  275.                     }
  276.                 }

  277.                 // Fill ah, aw, zhd and zwd lists
  278.                 if (line.length() > 0 && !line.startsWith("!")) {
  279.                     final String[] values_line = SEPARATOR.split(line);
  280.                     ah.add(Double.valueOf(values_line[2]));
  281.                     aw.add(Double.valueOf(values_line[3]));
  282.                     zhd.add(Double.valueOf(values_line[4]));
  283.                     zwd.add(Double.valueOf(values_line[5]));
  284.                 }
  285.             }

  286.         } catch (NumberFormatException nfe) {
  287.             throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
  288.                                       lineNumber, name, line);
  289.         }

  290.         final int dimLat = latitudes.size();
  291.         final int dimLon = longitudes.size();

  292.         // Change Lists to Arrays
  293.         final double[] xVal = new double[dimLat];
  294.         for (int i = 0; i < dimLat; i++) {
  295.             xVal[i] = latitudes.get(i);
  296.         }

  297.         final double[] yVal = new double[dimLon];
  298.         for (int j = 0; j < dimLon; j++) {
  299.             yVal[j] = longitudes.get(j);
  300.         }

  301.         final double[][] fvalAH = new double[dimLat][dimLon];
  302.         final double[][] fvalAW = new double[dimLat][dimLon];
  303.         final double[][] fvalZH = new double[dimLat][dimLon];
  304.         final double[][] fvalZW = new double[dimLat][dimLon];

  305.         int index = dimLon * dimLat;
  306.         for (int x = 0; x < dimLat; x++) {
  307.             for (int y = dimLon - 1; y >= 0; y--) {
  308.                 index = index - 1;
  309.                 fvalAH[x][y] = ah.get(index);
  310.                 fvalAW[x][y] = aw.get(index);
  311.                 fvalZH[x][y] = zhd.get(index);
  312.                 fvalZW[x][y] = zwd.get(index);
  313.             }
  314.         }

  315.         // Build Bilinear Interpolation Functions
  316.         final BilinearInterpolatingFunction functionAH = new BilinearInterpolatingFunction(xVal, yVal, fvalAH);
  317.         final BilinearInterpolatingFunction functionAW = new BilinearInterpolatingFunction(xVal, yVal, fvalAW);
  318.         final BilinearInterpolatingFunction functionZH = new BilinearInterpolatingFunction(xVal, yVal, fvalZH);
  319.         final BilinearInterpolatingFunction functionZW = new BilinearInterpolatingFunction(xVal, yVal, fvalZW);

  320.         coefficientsA = new double[2];
  321.         zenithDelay   = new double[2];

  322.         // Get the values for the given latitude and longitude
  323.         coefficientsA[0] = functionAH.value(latitude, longitude);
  324.         coefficientsA[1] = functionAW.value(latitude, longitude);
  325.         zenithDelay[0]   = functionZH.value(latitude, longitude);
  326.         zenithDelay[1]   = functionZW.value(latitude, longitude);

  327.         // Check that ah, aw, zh and zw were found
  328.         if (coefficientsA == null || zenithDelay == null) {
  329.             throw new OrekitException(OrekitMessages.NO_VIENNA_ACOEF_OR_ZENITH_DELAY_IN_FILE, name);
  330.         }

  331.     }

  332. }