OceanLoadingCoefficientsBlqParser.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.displacement;

import org.hipparchus.exception.LocalizedCoreFormats;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.data.DataSource;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitMessages;

import java.io.BufferedReader;
import java.io.IOException;
import java.util.ArrayList;
import java.util.List;
import java.util.regex.Matcher;
import java.util.regex.Pattern;

/**
 * Parser for ocean loading coefficients, using Onsala Space Observatory files in BLQ format.
 * <p>
 * Files in BLQ format can be generated using the form at the
 * <a href="http://holt.oso.chalmers.se/loading/">Bos-Scherneck web site</a>,
 * selecting BLQ as the output format.
 * </p>
 * <p>
 * The sites names are extracted from the file content, not the file name, because the
 * file can contain more than one station. As we expect existing files may have been
 * stripped from headers and footers, we do not attempt to parse them. We only parse
 * the series of 7 lines blocks starting with the lines with the station names and their
 * coordinates and the 6 data lines that follows. Several such blocks may appear in the
 * file. Copy-pasting the entire mail received from OSO after completing the web site
 * form works, as intermediate lines between the 7 lines blocks are simply ignored.
 * </p>
 * @see OceanLoadingCoefficients
 * @see OceanLoading
 * @since 9.1
 * @author Luc Maisonobe
 */
public class OceanLoadingCoefficientsBlqParser {

    /** Pattern for fields with real type. */
    private static final String  REAL_TYPE_PATTERN = "[-+]?(?:\\p{Digit}+(?:\\.\\p{Digit}*)?|\\.\\p{Digit}+)(?:[eE][-+]?\\p{Digit}+)?";

    /** Pattern for extracted real fields. */
    private static final String  REAL_FIELD_PATTERN = "\\p{Space}*(" + REAL_TYPE_PATTERN + ")";

    /** Pattern for end of line. */
    private static final String  END_OF_LINE_PATTERN = "\\p{Space}*$";

    /** Pattern for site name and coordinates lines. */
    private static final String  SITE_LINE_PATTERN = "^\\$\\$ *([^,]*),\\p{Space}*(?:RADI TANG)?\\p{Space}*lon/lat:" +
                                                     REAL_FIELD_PATTERN +
                                                     REAL_FIELD_PATTERN +
                                                     REAL_FIELD_PATTERN +
                                                     END_OF_LINE_PATTERN;

    /** Pattern for coefficients lines. */
    private static final String  DATA_LINE_PATTERN = "^" +
                                                     REAL_FIELD_PATTERN + // M₂ tide
                                                     REAL_FIELD_PATTERN + // S₂ tide
                                                     REAL_FIELD_PATTERN + // N₂ tide
                                                     REAL_FIELD_PATTERN + // K₂ tide
                                                     REAL_FIELD_PATTERN + // K₁ tide
                                                     REAL_FIELD_PATTERN + // O₁ tide
                                                     REAL_FIELD_PATTERN + // P₁ tide
                                                     REAL_FIELD_PATTERN + // Q₁ tide
                                                     REAL_FIELD_PATTERN + // Mf tide
                                                     REAL_FIELD_PATTERN + // Mm tide
                                                     REAL_FIELD_PATTERN + // Ssa tide
                                                     END_OF_LINE_PATTERN;

    /** Pattern for site name and coordinates lines. */
    private static final Pattern SITE_PATTERN = Pattern.compile(SITE_LINE_PATTERN);

    /** Pattern for coefficients lines. */
    private static final Pattern DATA_PATTERN = Pattern.compile(DATA_LINE_PATTERN);

    /** Main tides. */
    private static final Tide[][] TIDES = {
        {
            Tide.SSA, Tide.MM, Tide.MF
        }, {
            Tide.Q1,  Tide.O1, Tide.P1, Tide.K1
        }, {
            Tide.N2,  Tide.M2, Tide.S2, Tide.K2
        }
    };

    /** Species index for each column. */
    private static final int[] I = {
        2, 2, 2, 2, 1, 1, 1, 1, 0, 0, 0
    };

    /** Tides index for each column. */
    private static final int[] J = {
        1, 2, 0, 3, 3, 1, 2, 0, 2, 1, 0
    };

    /** Parse a BLQ file.
     * <p>
     * Files in BLQ format can be generated using the form at the
     * <a href="https://holt.oso.chalmers.se/loading/">Bos-Scherneck web site</a>,
     * selecting BLQ as the output format.
     * </p>
     * <p>
     * when completing the web site form, the email received as the following form:
     * </p>
     * <pre>{@literal
     * $$ Ocean loading displacement
     * $$
     * $$ Calculated on holt using olfg/olmpp of H.-G. Scherneck
     * $$
     * $$ COLUMN ORDER:  M2  S2  N2  K2  K1  O1  P1  Q1  MF  MM SSA
     * $$
     * $$ ROW ORDER:
     * $$ AMPLITUDES (m)
     * $$   RADIAL
     * $$   TANGENTL    EW
     * $$   TANGENTL    NS
     * $$ PHASES (degrees)
     * $$   RADIAL
     * $$   TANGENTL    EW
     * $$   TANGENTL    NS
     * $$
     * $$ Displacement is defined positive in upwards, South and West direction.
     * $$ The phase lag is relative to Greenwich and lags positive. The
     * $$ Gutenberg-Bullen Greens function is used. In the ocean tide model the
     * $$ deficit of tidal water mass has been corrected by subtracting a uniform
     * $$ layer of water with a certain phase lag globally.
     * $$
     * $$ Complete <model name> : No interpolation of ocean model was necessary
     * $$ <model name>_PP       : Ocean model has been interpolated near the station
     * $$                         (PP = Post-Processing)
     * $$
     * $$ Ocean tide model: CSR4.0, long-period tides from FES99
     * $$
     * $$ END HEADER
     * $$
     *   Goldstone
     * $$ Complete CSR4.0_f
     * $$ Computed by OLFG, H.-G. Scherneck, Onsala Space Observatory 2017-Sep-28
     * $$ Goldstone,                 RADI TANG  lon/lat:  243.1105   35.4259    0.000
     *   .00130 .00155 .00064 .00052 .01031 .00661 .00339 .00119 .00005 .00002 .00003
     *   .00136 .00020 .00024 .00004 .00322 .00202 .00106 .00036 .00007 .00003 .00001
     *   .00372 .00165 .00082 .00045 .00175 .00113 .00057 .00022 .00004 .00002 .00003
     *     -2.7 -106.3  -62.6 -106.8   41.6   27.3   40.4   24.0 -119.1 -123.2 -169.7
     *   -145.3  -88.4  178.5  -66.3 -130.5 -145.3 -131.7 -148.7  124.3  139.6   23.3
     *     90.7  111.1   74.1  111.3  176.9  165.3  175.8  164.0   48.9   25.3    4.5
     * $$
     *   ONSALA
     * $$ CSR4.0_f_PP ID: 2017-09-28 15:01:14
     * $$ Computed by OLMPP by H G Scherneck, Onsala Space Observatory, 2017
     * $$ Onsala,                    RADI TANG  lon/lat:   11.9264   57.3958    0.000
     *   .00344 .00121 .00078 .00031 .00189 .00116 .00064 .00004 .00090 .00048 .00041
     *   .00143 .00035 .00035 .00008 .00053 .00051 .00018 .00009 .00013 .00006 .00007
     *   .00086 .00023 .00023 .00006 .00029 .00025 .00010 .00008 .00003 .00001 .00000
     *    -64.6  -50.3  -95.0  -53.1  -58.8 -152.4  -65.5 -133.8    9.8    5.8    2.1
     *     85.4  115.2   56.7  114.7   99.5   15.9   94.2  -10.0 -166.3 -169.8 -177.7
     *    110.7  147.1   93.9  148.6   49.4  -56.5   34.8 -169.9  -35.3   -3.7   10.1
     * $$ END TABLE
     * Errors:
     * Warnings:
     * }</pre>
     * <p>
     * We only parse blocks 7 lines blocks starting with the lines with the station names
     * and their coordinates and the 6 data lines that follows. Several such blocks may
     * appear in the file.
     * </p>
     * @param source source for BLQ data
     * @return parsed coefficients
     */
    public List<OceanLoadingCoefficients> parse(final DataSource source) {

        final List<OceanLoadingCoefficients> coefficients = new ArrayList<>();

        // temporary holders for parsed data
        String siteName = null;
        GeodeticPoint siteLocation = null;
        final double[][][] data = new double[6][3][];
        for (int i = 0; i < data.length; ++i) {
            data[i][0] = new double[3];
            data[i][1] = new double[4];
            data[i][2] = new double[4];
        }

        try (BufferedReader reader = new BufferedReader(
            source.getOpener().openReaderOnce())) {

            int lineNumber = 0;
            int dataLine = -1;
            for (String line = reader.readLine(); line != null; line = reader.readLine()) {
                ++lineNumber;

                if (dataLine < 0) {
                    // we are looking for a site line
                    final Matcher matcher = SITE_PATTERN.matcher(line);
                    if (matcher.matches()) {
                        // the current line is a site description line
                        siteName = matcher.group(1);
                        siteLocation =
                            new GeodeticPoint(FastMath.toRadians(Double.parseDouble(matcher.group(3))),
                                              FastMath.toRadians(Double.parseDouble(matcher.group(2))),
                                              Double.parseDouble(matcher.group(4)));
                        // next line must be line 0 of the data
                        dataLine = 0;
                    }
                } else {
                    // we are looking for a data line
                    final Matcher matcher = DATA_PATTERN.matcher(line);
                    if (!matcher.matches()) {
                        throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
                                                  lineNumber, source.getName(), line);
                    }
                    for (int k = 0; k < I.length; ++k) {
                        if (dataLine < 3) {
                            // amplitude data
                            data[dataLine][I[k]][J[k]] = Double.parseDouble(matcher.group(k + 1));
                        } else {
                            // phase data (reversed to be negative for lags)
                            data[dataLine][I[k]][J[k]] = -FastMath.toRadians(Double.parseDouble(matcher.group(k + 1)));
                        }
                    }
                    if (dataLine < data.length - 1) {
                        // we need more data
                        ++dataLine;
                    } else {
                        // it was the last data line
                        coefficients.add(new OceanLoadingCoefficients(siteName, siteLocation,
                                                                      TIDES, data[0],
                                                                      data[3], data[1],
                                                                      data[4], data[2],
                                                                      data[5]));
                        dataLine = -1;
                    }
                }
            }

            if (dataLine >= 0) {
                // we were looking for a line that did not appear
                throw new OrekitException(OrekitMessages.UNEXPECTED_END_OF_FILE_AFTER_LINE,
                                          source.getName(), lineNumber);
            }

        } catch (IOException ioe) {
            throw new OrekitException(ioe, LocalizedCoreFormats.SIMPLE_MESSAGE,
                                      ioe.getLocalizedMessage());
        }

        return coefficients;

    }

}