initial commit

96813390_header.24O

@@ -0,0 +1,25 @@
+     3.02           OBSERVATION DATA    M (MIXED)           RINEX VERSION / TYPE
+NetR9 5.01          Receiver Operator   04-DEC-24 00:00:00  PGM / RUN BY / DATE
+TQBL                                                        MARKER NAME
+TQBL                                                        MARKER NUMBER
+GEODETIC                                                    MARKER TYPE         
+TQBL                TQBL                                    OBSERVER / AGENCY   
+5034K69681          Trimble NetR9       5.01                REC # / TYPE / VERS
+                    TRM59900.00     NONE                    ANT # / TYPE
+ -1626343.5660  5730606.2916  2271881.4271                  APPROX POSITION XYZ
+        0.0000        0.0000        0.0000                  ANTENNA: DELTA H/E/N
+G   12 C1C L1C S1C C2W L2W S2W C2X L2X S2X C5X L5X S5X      SYS / # / OBS TYPES
+S    3 C1C L1C S1C                                          SYS / # / OBS TYPES
+R   12 C1C L1C S1C C1P L1P S1P C2C L2C S2C C2P L2P S2P      SYS / # / OBS TYPES
+E   12 C1X L1X S1X C5X L5X S5X C7X L7X S7X C8X L8X S8X      SYS / # / OBS TYPES
+C    6 C1I L1I S1I C7I L7I S7I                              SYS / # / OBS TYPES
+    10.000                                                  INTERVAL
+  2024    12     4     0     0    0.0000000     GPS         TIME OF FIRST OBS
+G L2X -0.25000                                              SYS / PHASE SHIFT   
+R L1P  0.25000                                              SYS / PHASE SHIFT   
+R L2C -0.25000                                              SYS / PHASE SHIFT   
+J L2X -0.25000                                              SYS / PHASE SHIFT   
+GIOVE-A if present is mapped to satellite ID 51             COMMENT             
+GIOVE-B if present is mapped to satellite ID 52             COMMENT             
+DBHZ                                                        SIGNAL STRENGTH UNIT
+                                                            END OF HEADER       

ephemeris.py

@@ -0,0 +1,89 @@
+import georinex as gr
+import numpy as np
+from datetime import datetime
+import math
+
+def ecef_to_lla(x, y, z):
+    """Convert ECEF coordinates to LLA (Latitude, Longitude, Altitude)"""
+    a = 6378137.0  # WGS84 semi-major axis
+    f = 1/298.257223563  # WGS84 flattening
+    e2 = 2*f - f**2  # eccentricity squared
+    
+    # Calculate longitude
+    lon = math.atan2(y, x)
+    
+    # Initial latitude guess
+    p = math.sqrt(x**2 + y**2)
+    lat = math.atan2(z, p * (1-e2))
+    
+    # Iterative latitude calculation
+    for _ in range(5):
+        N = a / math.sqrt(1 - e2 * math.sin(lat)**2)
+        h = p / math.cos(lat) - N
+        lat = math.atan2(z, p * (1 - e2 * N/(N + h)))
+    
+    return np.degrees(lat), np.degrees(lon), h
+
+def calculate_azimuth_elevation(receiver_xyz, satellite_xyz):
+    """Calculate azimuth and elevation angles between receiver and satellite"""
+    # Convert receiver ECEF to LLA
+    rec_lat, rec_lon, rec_h = ecef_to_lla(receiver_xyz[0], receiver_xyz[1], receiver_xyz[2])
+    rec_lat_rad = math.radians(rec_lat)
+    rec_lon_rad = math.radians(rec_lon)
+    
+    # Calculate local vector between receiver and satellite
+    dx = satellite_xyz[0] - receiver_xyz[0]
+    dy = satellite_xyz[1] - receiver_xyz[1]
+    dz = satellite_xyz[2] - receiver_xyz[2]
+    
+    # Transform to ENU coordinates
+    E = -math.sin(rec_lon_rad)*dx + math.cos(rec_lon_rad)*dy
+    N = -math.sin(rec_lat_rad)*math.cos(rec_lon_rad)*dx - \
+         math.sin(rec_lat_rad)*math.sin(rec_lon_rad)*dy + \
+         math.cos(rec_lat_rad)*dz
+    U = math.cos(rec_lat_rad)*math.cos(rec_lon_rad)*dx + \
+        math.cos(rec_lat_rad)*math.sin(rec_lon_rad)*dy + \
+        math.sin(rec_lat_rad)*dz
+    
+    # Calculate azimuth and elevation
+    azimuth = math.degrees(math.atan2(E, N)) % 360
+    elevation = math.degrees(math.asin(U/math.sqrt(E**2 + N**2 + U**2)))
+    
+    return azimuth, elevation
+
+def main():
+    # Read navigation and observation files
+    nav_file = "96813390.24N"
+    obs_file = "96813390_small.24O"
+    
+    nav = gr.load(nav_file)
+    obs = gr.load(obs_file)
+    
+    # Get receiver position from observation file header
+    receiver_pos = obs.position
+    
+    # Process each epoch
+    for time in obs.time:
+        print(f"\nTime: {time}")
+        
+        # Get visible satellites at this epoch
+        visible_sats = obs.sel(time=time).sv.values
+        
+        for sat in visible_sats:
+            if sat[0] == 'G':  # Process only GPS satellites
+                # Get satellite position from navigation data
+                try:
+                    sat_pos = nav.sel(sv=sat, time=time)
+                    sat_xyz = np.array([sat_pos.position[0], 
+                                      sat_pos.position[1], 
+                                      sat_pos.position[2]])
+                    
+                    # Calculate azimuth and elevation
+                    az, el = calculate_azimuth_elevation(receiver_pos, sat_xyz)
+                    
+                    print(f"Satellite {sat}: Azimuth = {az:.2f}°, Elevation = {el:.2f}°")
+                except:
+                    continue
+
+if __name__ == "__main__":
+    main()

pseudorange.py

@@ -0,0 +1,61 @@
+import georinex as gr
+import matplotlib.pyplot as plt
+import numpy as np
+from datetime import datetime
+import pandas as pd
+
+# Read the RINEX file
+rinex_file = '96813390_small.24O'
+data = gr.load(rinex_file)
+
+# Get available GPS satellites, start with 'G'
+available_sats = [sv for sv in data.sv.values if sv.startswith('G')]
+
+# Extract GPS pseudorange data
+gps_data = data.sel(sv=available_sats)
+# For now, get only the C1C pseudorange data
+c1c_data = gps_data['C1C'].dropna(dim='sv', how='all')
+
+# Convert time index to datetime
+times = c1c_data.time.values
+satellites = c1c_data.sv.values
+
+# Create the plot
+plt.figure(figsize=(15, 10))
+
+# Plot each satellite's pseudorange
+for sat in satellites:
+    sat_data = c1c_data.sel(sv=sat)
+    valid_data = sat_data.dropna(dim='time')
+    
+    if len(valid_data) > 0:  # Only plot if we have valid data
+        plt.plot(valid_data.time, valid_data.values, 
+                 label=sat, marker='.',
+                 linestyle='-', markersize=2)
+
+# Customize the plot
+plt.title('GPS Pseudorange Measurements (C1C)', fontsize=14)
+plt.xlabel('Time (UTC)', fontsize=12)
+plt.ylabel('Pseudorange', fontsize=12)
+plt.grid(True, linestyle='--', alpha=0.7)
+plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
+plt.tight_layout()
+
+# Rotate x-axis labels for better readability
+plt.xticks(rotation=45)
+
+# Show the plot
+plt.show()
+
+# Print some basic statistics
+print("\nBasic Statistics for GPS Pseudorange Measurements:")
+print("Number of GPS satellites observed:", len(satellites))
+print("Observation period:", times[0], "to", times[-1])
+print("\nPseudorange ranges for each satellite:")
+for sat in satellites:
+    sat_data = c1c_data.sel(sv=sat).dropna(dim='time')
+    if len(sat_data) > 0:
+        print(f"{sat}: {sat_data.min().values:.2f} - {sat_data.max().values:.2f} meters")
+
+# Optionally, print available satellites
+print("\nAvailable GPS satellites:", available_sats)

requirements.txt

@@ -0,0 +1,5 @@
+georinex
+matplotlib
+numpy
+pandas
+xarray