1   /* Copyright 2002-2021 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.displacement;
18  
19  import java.io.BufferedReader;
20  import java.io.IOException;
21  import java.io.InputStream;
22  import java.io.InputStreamReader;
23  import java.nio.charset.StandardCharsets;
24  import java.util.ArrayList;
25  import java.util.List;
26  import java.util.Optional;
27  import java.util.regex.Matcher;
28  import java.util.regex.Pattern;
29  import java.util.stream.Collectors;
30  
31  import org.hipparchus.util.FastMath;
32  import org.orekit.annotation.DefaultDataContext;
33  import org.orekit.bodies.GeodeticPoint;
34  import org.orekit.data.AbstractSelfFeedingLoader;
35  import org.orekit.data.DataContext;
36  import org.orekit.data.DataLoader;
37  import org.orekit.data.DataProvidersManager;
38  import org.orekit.errors.OrekitException;
39  import org.orekit.errors.OrekitMessages;
40  
41  /**
42   * Factory for ocean loading coefficients, using Onsala Space Observatory files in BLQ format.
43   * <p>
44   * Files in BLQ format can be generated using the form at the
45   * <a href="http://holt.oso.chalmers.se/loading/">Bos-Scherneck web site</a>,
46   * selecting BLQ as the output format.
47   * </p>
48   * <p>
49   * The sites names are extracted from the file content, not the file name, because the
50   * file can contain more than one station. As we expect existing files may have been
51   * stripped from headers and footers, we do not attempt to parse them. We only parse
52   * the series of 7 lines blocks starting with the lines with the station names and their
53   * coordinates and the 6 data lines that follows. Several such blocks may appear in the
54   * file. Copy-pasting the entire mail received from OSO after completing the web site
55   * form works, as intermediate lines between the 7 lines blocks are simply ignored.
56   * </p>
57   * @see OceanLoadingCoefficients
58   * @see OceanLoading
59   * @since 9.1
60   * @author Luc Maisonobe
61   */
62  public class OceanLoadingCoefficientsBLQFactory extends AbstractSelfFeedingLoader {
63  
64      /** Default supported files name pattern for Onsala Space Observatory files in BLQ format. */
65      public static final String DEFAULT_BLQ_SUPPORTED_NAMES = "^.+\\.blq$";
66  
67      /** Pattern for fields with real type. */
68      private static final String  REAL_TYPE_PATTERN = "[-+]?(?:(?:\\p{Digit}+(?:\\.\\p{Digit}*)?)|(?:\\.\\p{Digit}+))(?:[eE][-+]?\\p{Digit}+)?";
69  
70      /** Pattern for extracted real fields. */
71      private static final String  REAL_FIELD_PATTERN = "\\p{Space}*(" + REAL_TYPE_PATTERN + ")";
72  
73      /** Pattern for end of line. */
74      private static final String  END_OF_LINE_PATTERN = "\\p{Space}*$";
75  
76      /** Pattern for site name and coordinates lines. */
77      private static final String  SITE_LINE_PATTERN = "^\\$\\$ *([^,]*),\\p{Space}*(?:RADI TANG)?\\p{Space}*lon/lat:" +
78                                                       REAL_FIELD_PATTERN +
79                                                       REAL_FIELD_PATTERN +
80                                                       REAL_FIELD_PATTERN +
81                                                       END_OF_LINE_PATTERN;
82  
83      /** Pattern for coefficients lines. */
84      private static final String  DATA_LINE_PATTERN = "^" +
85                                                       REAL_FIELD_PATTERN + // M₂ tide
86                                                       REAL_FIELD_PATTERN + // S₂ tide
87                                                       REAL_FIELD_PATTERN + // N₂ tide
88                                                       REAL_FIELD_PATTERN + // K₂ tide
89                                                       REAL_FIELD_PATTERN + // K₁ tide
90                                                       REAL_FIELD_PATTERN + // O₁ tide
91                                                       REAL_FIELD_PATTERN + // P₁ tide
92                                                       REAL_FIELD_PATTERN + // Q₁ tide
93                                                       REAL_FIELD_PATTERN + // Mf tide
94                                                       REAL_FIELD_PATTERN + // Mm tide
95                                                       REAL_FIELD_PATTERN + // Ssa tide
96                                                       END_OF_LINE_PATTERN;
97  
98      /** Pattern for site name and coordinates lines. */
99      private static final Pattern SITE_PATTERN = Pattern.compile(SITE_LINE_PATTERN);
100 
101     /** Pattern for coefficients lines. */
102     private static final Pattern DATA_PATTERN = Pattern.compile(DATA_LINE_PATTERN);
103 
104     /** Main tides. */
105     private static final Tide[][] TIDES = {
106         {
107             Tide.SSA, Tide.MM, Tide.MF
108         }, {
109             Tide.Q1,  Tide.O1, Tide.P1, Tide.K1
110         }, {
111             Tide.N2,  Tide.M2, Tide.S2, Tide.K2
112         }
113     };
114 
115     /** Species index for each column. */
116     private static final int[] I = {
117         2, 2, 2, 2, 1, 1, 1, 1, 0, 0, 0
118     };
119 
120     /** Tides index for each column. */
121     private static final int[] J = {
122         1, 2, 0, 3, 3, 1, 2, 0, 2, 1, 0
123     };
124 
125     /** Parsed coefficients. */
126     private final List<OceanLoadingCoefficients> coefficients;
127 
128     /** Simple constructor. This constructor uses the {@link DataContext#getDefault()
129      * default data context}.
130      * <p>
131      * Files in BLQ format can be generated using the form at the
132      * <a href="http://holt.oso.chalmers.se/loading/">Bos-Scherneck web site</a>,
133      * selecting BLQ as the output format.
134      * </p>
135      * @param supportedNames regular expression for supported files names
136      * @see #DEFAULT_BLQ_SUPPORTED_NAMES
137      * @see #OceanLoadingCoefficientsBLQFactory(String, DataProvidersManager)
138      */
139     @DefaultDataContext
140     public OceanLoadingCoefficientsBLQFactory(final String supportedNames) {
141         this(supportedNames, DataContext.getDefault().getDataProvidersManager());
142     }
143 
144     /**
145      * This constructor allows specification of the source of the BLQ auxiliary data
146      * files.
147      *
148      * <p>
149      * Files in BLQ format can be generated using the form at the
150      * <a href="http://holt.oso.chalmers.se/loading/">Bos-Scherneck web site</a>,
151      * selecting BLQ as the output format.
152      * </p>
153      * @param supportedNames regular expression for supported files names
154      * @param dataProvidersManager provides access to auxiliary data files.
155      * @see #DEFAULT_BLQ_SUPPORTED_NAMES
156      * @since 10.1
157      */
158     public OceanLoadingCoefficientsBLQFactory(
159             final String supportedNames,
160             final DataProvidersManager dataProvidersManager) {
161         super(supportedNames, dataProvidersManager);
162 
163         this.coefficients   = new ArrayList<>();
164 
165     }
166 
167     /** Lazy loading of coefficients.
168      */
169     private void loadsIfNeeded() {
170         if (coefficients.isEmpty()) {
171             feed(new BLQParser());
172         }
173     }
174 
175     /** Get the list of sites for which we have found coefficients, in lexicographic order ignoring case.
176      * @return list of sites for which we have found coefficients, in lexicographic order ignoring case
177      */
178     public List<String> getSites() {
179 
180         loadsIfNeeded();
181 
182         // extract sites names from the map
183         final List<String> sites = coefficients.stream()
184                 .map(OceanLoadingCoefficients::getSiteName)
185                 // sort to ensure we have a reproducible order
186                 .sorted(String::compareToIgnoreCase)
187                 .collect(Collectors.toList());
188 
189         return sites;
190 
191     }
192 
193     /** Get the coefficients for a given site.
194      * @param site site name (as it appears in the Onsala Space Observatory files in BLQ format),
195      * ignoring case
196      * @return coefficients for the site
197      */
198     public OceanLoadingCoefficients getCoefficients(final String site) {
199 
200         loadsIfNeeded();
201 
202         final Optional<OceanLoadingCoefficients> optional =
203                         coefficients.stream().filter(c -> c.getSiteName().equalsIgnoreCase(site)).findFirst();
204         if (!optional.isPresent()) {
205             throw new OrekitException(OrekitMessages.STATION_NOT_FOUND,
206                                       site,
207                                       String.join(", ", getSites()));
208         }
209 
210         return optional.get();
211 
212     }
213 
214     /** Local parser for the Onsala Space Observatory BLQ files.
215      * <p>
216      * when completing the web site form, the email received as the following form:
217      * </p>
218      * <pre>{@literal
219      * $$ Ocean loading displacement
220      * $$
221      * $$ Calculated on holt using olfg/olmpp of H.-G. Scherneck
222      * $$
223      * $$ COLUMN ORDER:  M2  S2  N2  K2  K1  O1  P1  Q1  MF  MM SSA
224      * $$
225      * $$ ROW ORDER:
226      * $$ AMPLITUDES (m)
227      * $$   RADIAL
228      * $$   TANGENTL    EW
229      * $$   TANGENTL    NS
230      * $$ PHASES (degrees)
231      * $$   RADIAL
232      * $$   TANGENTL    EW
233      * $$   TANGENTL    NS
234      * $$
235      * $$ Displacement is defined positive in upwards, South and West direction.
236      * $$ The phase lag is relative to Greenwich and lags positive. The
237      * $$ Gutenberg-Bullen Greens function is used. In the ocean tide model the
238      * $$ deficit of tidal water mass has been corrected by subtracting a uniform
239      * $$ layer of water with a certain phase lag globally.
240      * $$
241      * $$ Complete <model name> : No interpolation of ocean model was necessary
242      * $$ <model name>_PP       : Ocean model has been interpolated near the station
243      * $$                         (PP = Post-Processing)
244      * $$
245      * $$ Ocean tide model: CSR4.0, long-period tides from FES99
246      * $$
247      * $$ END HEADER
248      * $$
249      *   Goldstone
250      * $$ Complete CSR4.0_f
251      * $$ Computed by OLFG, H.-G. Scherneck, Onsala Space Observatory 2017-Sep-28
252      * $$ Goldstone,                 RADI TANG  lon/lat:  243.1105   35.4259    0.000
253      *   .00130 .00155 .00064 .00052 .01031 .00661 .00339 .00119 .00005 .00002 .00003
254      *   .00136 .00020 .00024 .00004 .00322 .00202 .00106 .00036 .00007 .00003 .00001
255      *   .00372 .00165 .00082 .00045 .00175 .00113 .00057 .00022 .00004 .00002 .00003
256      *     -2.7 -106.3  -62.6 -106.8   41.6   27.3   40.4   24.0 -119.1 -123.2 -169.7
257      *   -145.3  -88.4  178.5  -66.3 -130.5 -145.3 -131.7 -148.7  124.3  139.6   23.3
258      *     90.7  111.1   74.1  111.3  176.9  165.3  175.8  164.0   48.9   25.3    4.5
259      * $$
260      *   ONSALA
261      * $$ CSR4.0_f_PP ID: 2017-09-28 15:01:14
262      * $$ Computed by OLMPP by H G Scherneck, Onsala Space Observatory, 2017
263      * $$ Onsala,                    RADI TANG  lon/lat:   11.9264   57.3958    0.000
264      *   .00344 .00121 .00078 .00031 .00189 .00116 .00064 .00004 .00090 .00048 .00041
265      *   .00143 .00035 .00035 .00008 .00053 .00051 .00018 .00009 .00013 .00006 .00007
266      *   .00086 .00023 .00023 .00006 .00029 .00025 .00010 .00008 .00003 .00001 .00000
267      *    -64.6  -50.3  -95.0  -53.1  -58.8 -152.4  -65.5 -133.8    9.8    5.8    2.1
268      *     85.4  115.2   56.7  114.7   99.5   15.9   94.2  -10.0 -166.3 -169.8 -177.7
269      *    110.7  147.1   93.9  148.6   49.4  -56.5   34.8 -169.9  -35.3   -3.7   10.1
270      * $$ END TABLE
271      * Errors:
272      * Warnings:
273      * }</pre>
274      * <p>
275      * We only parse blocks 7 lines blocks starting with the lines with the station names
276      * and their coordinates and the 6 data lines that follows. Several such blocks may
277      * appear in the file.
278      * </p>
279      */
280     private class BLQParser implements DataLoader {
281 
282         /** Simple constructor.
283          */
284         BLQParser() {
285             // empty constructor
286         }
287 
288         /** {@inheritDoc} */
289         @Override
290         public boolean stillAcceptsData() {
291             // as data for different stations may be in different files
292             // we always accept new data, even if we have already parsed
293             // some files
294             return true;
295         }
296 
297         /** {@inheritDoc} */
298         @Override
299         public void loadData(final InputStream input, final String name)
300             throws IOException, OrekitException {
301 
302             // temporary holders for parsed data
303             String         siteName     = null;
304             GeodeticPoint  siteLocation = null;
305             final double[][][] data     = new double[6][3][];
306             for (int i = 0; i < data.length; ++i) {
307                 data[i][0] = new double[3];
308                 data[i][1] = new double[4];
309                 data[i][2] = new double[4];
310             }
311 
312             try (BufferedReader reader = new BufferedReader(new InputStreamReader(input, StandardCharsets.UTF_8))) {
313 
314                 int     lineNumber = 0;
315                 int     dataLine   = -1;
316                 for (String line = reader.readLine(); line != null; line = reader.readLine()) {
317                     ++lineNumber;
318 
319                     if (dataLine < 0) {
320                         // we are looking for a site line
321                         final Matcher matcher = SITE_PATTERN.matcher(line);
322                         if (matcher.matches()) {
323                             // the current line is a site description line
324                             siteName = matcher.group(1);
325                             siteLocation = new GeodeticPoint(FastMath.toRadians(Double.parseDouble(matcher.group(3))),
326                                                              FastMath.toRadians(Double.parseDouble(matcher.group(2))),
327                                                              Double.parseDouble(matcher.group(4)));
328                             // next line must be line 0 of the data
329                             dataLine = 0;
330                         }
331                     } else {
332                         // we are looking for a data line
333                         final Matcher matcher = DATA_PATTERN.matcher(line);
334                         if (!matcher.matches()) {
335                             throw new OrekitException(OrekitMessages.UNABLE_TO_PARSE_LINE_IN_FILE,
336                                                       lineNumber, name, line);
337                         }
338                         for (int k = 0; k < I.length; ++k) {
339                             if (dataLine < 3) {
340                                 // amplitude data
341                                 data[dataLine][I[k]][J[k]] = Double.parseDouble(matcher.group(k + 1));
342                             } else {
343                                 // phase data (reversed to be negative for lags)
344                                 data[dataLine][I[k]][J[k]] = -FastMath.toRadians(Double.parseDouble(matcher.group(k + 1)));
345                             }
346                         }
347                         if (dataLine < data.length - 1) {
348                             // we need more data
349                             ++dataLine;
350                         } else {
351                             // it was the last data line
352                             coefficients.add(new OceanLoadingCoefficients(siteName, siteLocation, TIDES,
353                                                                           data[0], data[3],
354                                                                           data[1], data[4],
355                                                                           data[2], data[5]));
356                             dataLine = -1;
357                         }
358                     }
359                 }
360 
361                 if (dataLine >= 0) {
362                     // we were looking for a line that did not appear
363                     throw new OrekitException(OrekitMessages.UNEXPECTED_END_OF_FILE_AFTER_LINE,
364                                               name, lineNumber);
365                 }
366 
367             }
368         }
369 
370     }
371 
372 }