tolerate multiple blank header lines starting file

made sp3 reading more robust
geospace-code · Feb 10, 2020 · a4a0f5e · a4a0f5e
1 parent 30be8e5
commit a4a0f5e
Show file tree

Hide file tree

Showing 10 changed files with 3,196 additions and 35 deletions.
diff --git a/georinex/rio.py b/georinex/rio.py
@@ -36,7 +36,7 @@ def opener(fn: typing.Union[TextIO, Path],
 
         if fn.suffix == '.gz':
             with gzip.open(fn, 'rt') as f:
-                version, is_crinex = rinex_version(f.readline(80))
+                version, is_crinex = rinex_version(first_nonblank_line(f))
                 f.seek(0)
 
                 if is_crinex and not header:
@@ -57,7 +57,7 @@ def opener(fn: typing.Union[TextIO, Path],
                     yield f
         else:  # assume not compressed (or Hatanaka)
             with fn.open('r', encoding='ascii', errors='ignore') as f:
-                version, is_crinex = rinex_version(f.readline(80))
+                version, is_crinex = rinex_version(first_nonblank_line(f))
                 f.seek(0)
 
                 if is_crinex and not header:
@@ -67,6 +67,28 @@ def opener(fn: typing.Union[TextIO, Path],
         raise OSError(f'Unsure what to do with input of type: {type(fn)}')
 
 
+def first_nonblank_line(f: TextIO, max_lines: int = 10) -> str:
+    """ return first non-blank 80 character line in file
+
+    Parameters
+    ----------
+
+    max_lines: int
+        maximum number of blank lines
+    """
+
+    line = ""
+    for _i in range(max_lines):
+        line = f.readline(80)
+        if line.strip():
+            break
+
+    if _i == max_lines - 1 or not line:
+        raise ValueError(f"could not find first valid header line in {f.name}")
+
+    return line
+
+
 def rinexinfo(f: typing.Union[Path, TextIO]) -> typing.Dict[str, typing.Any]:
     """verify RINEX version"""
 
@@ -90,7 +112,7 @@ def rinexinfo(f: typing.Union[Path, TextIO]) -> typing.Dict[str, typing.Any]:
     f.seek(0)
 
     try:
-        line = f.readline(80)  # don't choke on binary files
+        line = first_nonblank_line(f)  # don't choke on binary files
 
         if line.startswith('#c'):
             return {'version': 'c',
@@ -150,7 +172,7 @@ def rinex_version(s: str) -> typing.Tuple[typing.Union[float, str], bool]:
     if not isinstance(s, str):
         raise TypeError('need first line of RINEX file as string')
     if len(s) < 2:
-        raise ValueError(f'first line of file is corrupted {s}')
+        raise ValueError(f'cannot decode RINEX version from line:\n{s}')
 
     if len(s) >= 80:
         if s[60:80] not in ('RINEX VERSION / TYPE', 'CRINEX VERS   / TYPE'):

diff --git a/georinex/sp3.py b/georinex/sp3.py
@@ -5,6 +5,7 @@
 import xarray
 import numpy as np
 import logging
+from .rio import first_nonblank_line
 from pathlib import Path
 from datetime import datetime
 import typing
@@ -16,10 +17,11 @@
 def load_sp3(fn: Path, outfn: Path) -> xarray.Dataset:
     dat: typing.Dict[str, typing.Any] = {}
     with fn.open("r") as f:
-        ln = f.readline()
+        ln = first_nonblank_line(f)
         assert ln[0] == "#", f"failed to read {fn} line 1"
         dat["t0"] = sp3dt(ln)
-        Nepoch = int(ln[32:39])
+        # Nepoch != number of time steps, at least for some files
+        dat["Nepoch"] = int(ln[32:39])
         dat["coord_sys"] = ln[46:51]
         dat["orbit_type"] = ln[52:55]
         dat["agency"] = ln[56:60]
@@ -45,30 +47,35 @@ def load_sp3(fn: Path, outfn: Path) -> xarray.Dataset:
         if not ln.startswith("*"):  # EOF
             raise ValueError(f"{fn} appears to be badly malformed")
         # the rest of the file is data, punctuated by epoch lines
-        ecef = np.empty((Nepoch, Nsv, 3))
-        ecef.fill(np.nan)
-        clock = np.empty((Nepoch, Nsv, 2))
-        clock.fill(np.nan)
-        vel = np.empty((Nepoch, Nsv, 3))
-        vel.fill(np.nan)
-        times = [sp3dt(ln)]
+        ecefs = []
+        clocks = []
+        vels = []
+        ecef = np.empty((Nsv, 3))
+        clock = np.empty((Nsv, 2))
+        vel = np.empty((Nsv, 3))
         i = 0
-        j = 0
+
+        times = [sp3dt(ln)]
 
         for ln in f:
             if ln[0] == "*":
                 times.append(sp3dt(ln))
-                i += 1
-                j = 0
+                ecefs.append(ecef)
+                clocks.append(clock)
+                vels.append(vel)
+                ecef = np.empty((Nsv, 3))
+                clock = np.empty((Nsv, 2))
+                vel = np.empty((Nsv, 3))
+                i = 0
                 continue
 
             if ln[0] == "P":
-                ecef[i, j, :] = (float(ln[4:18]), float(ln[18:32]), float(ln[32:46]))
-                clock[i, j, 0] = float(ln[46:60])
-                j += 1
+                ecef[i, :] = (float(ln[4:18]), float(ln[18:32]), float(ln[32:46]))
+                clock[i, 0] = float(ln[46:60])
+                i += 1
             elif ln[0] == "V":
-                vel[i, j - 1, :] = (float(ln[4:18]), float(ln[18:32]), float(ln[32:46]))
-                clock[i, j - 1, 1] = float(ln[46:60])
+                vel[i - 1, :] = (float(ln[4:18]), float(ln[18:32]), float(ln[32:46]))
+                clock[i - 1, 1] = float(ln[46:60])
             elif ln[:2] in ("EP", "EV"):
                 # let us know if you want these data types
                 pass
@@ -79,13 +86,19 @@ def load_sp3(fn: Path, outfn: Path) -> xarray.Dataset:
             else:
                 logging.info(f"unknown data {ln}")
 
+    # assemble the last time step
+    ecefs.append(ecef)
+    clocks.append(clock)
+    vels.append(vel)
+    aclock = np.asarray(clocks)
+
     # assemble into final xarray.Dataset
     ds = xarray.Dataset(coords={"time": times, "sv": svs, "ECEF": ["x", "y", "z"]})
-    ds["position"] = (("time", "sv", "ECEF"), ecef)
-    ds["clock"] = (("time", "sv"), clock[:, :, 0])
+    ds["position"] = (("time", "sv", "ECEF"), ecefs)
+    ds["clock"] = (("time", "sv"), aclock[:, :, 0])
     if not np.isnan(vel).all():
-        ds["velocity"] = (("time", "sv", "ECEF"), vel)
-        ds["dclock"] = (("time", "sv"), clock[:, :, 1])
+        ds["velocity"] = (("time", "sv", "ECEF"), vels)
+        ds["dclock"] = (("time", "sv"), aclock[:, :, 1])
 
     ds.attrs = dat
 

diff --git a/tests/data/blank1st.10n b/tests/data/blank1st.10n
@@ -0,0 +1,3 @@
+
+     2.11           N: GPS NAV. MESSAGE                     RINEX VERSION / TYPE
+                                                            END OF HEADER
diff --git a/tests/data/blank1st.10o b/tests/data/blank1st.10o
@@ -0,0 +1,5 @@
+
+     2.11           OBSERVATION DATA    M (MIXED)           RINEX VERSION / TYPE
+     7    L1    L2    P1    P2    C1    S1    S2            # / TYPES OF OBSERV
+  2010     3     5     0     0     0.0000000     GPS        TIME OF FIRST OBS
+                                                            END OF HEADER
diff --git a/tests/data/blank3_1st.10n b/tests/data/blank3_1st.10n
@@ -0,0 +1,3 @@
+
+     3.01           N: GNSS NAV DATA    S: SBAS             RINEX VERSION / TYPE
+                                                            END OF HEADER
diff --git a/tests/data/blank3_1st.10o b/tests/data/blank3_1st.10o
@@ -0,0 +1,7 @@
+
+     3.01           OBSERVATION DATA    M (MIXED)           RINEX VERSION / TYPE
+G    7 L1C L2P C1P C2P C1C S1P S2P                          SYS / # / OBS TYPES
+R    3 L1C C1C S1C                                          SYS / # / OBS TYPES
+S    3 L1C C1C S1C                                          SYS / # / OBS TYPES
+  2010     3     5     0     0     0.0000000     GPS        TIME OF FIRST OBS
+                                                            END OF HEADER
Original file line number	Diff line number	Diff line change
		@@ -0,0 +1,3 @@

		2.11 N: GPS NAV. MESSAGE RINEX VERSION / TYPE
		END OF HEADER
Original file line number	Diff line number	Diff line change
		@@ -0,0 +1,3 @@

		3.01 N: GNSS NAV DATA S: SBAS RINEX VERSION / TYPE
		END OF HEADER