ams: scripts to compte satellite tracks for each day. Now focusing on…

CHANGELOG.md

@@ -4,6 +4,7 @@
 
 ### Added
 
+- scripts to plot satellite tracks for each day - focused on ASIA-AQ regions, but can be extended.
 - ASIA-AQ mission directory initialized with FLUID instance from FIREX-AQ
 - pyabc module added in GMAO_aeropybs. the snket class is needed in pyabc for loading ffnet .net files
 - VIIRS and MODIS NNR angstrom exponent predictions are now saved in gridded output files (Level 3)

src/Components/missions/ASIA-AQ/TLE/aqua.tle

@@ -0,0 +1,2 @@
+1 27424U 02022A   24035.85570116  .00001446  00000-0  31909-3 0  9992
+2 27424  98.3154 344.9007 0001509 108.7279 310.6966 14.58828120157279

src/Components/missions/ASIA-AQ/TLE/noaa20.tle

@@ -0,0 +1,2 @@
+1 43013U 17073A   24035.86217735  .00000126  00000-0  80368-4 0  9990
+2 43013  98.7367 335.8403 0001157 115.1763 244.9533 14.19541299321985

src/Components/missions/ASIA-AQ/TLE/s5p.tle

@@ -0,0 +1,2 @@
+1 42969U 17064A   24035.40636003  .00000083  00000-0  59952-4 0  9999
+2 42969  98.7287 336.6467 0001638  78.2149 281.9211 14.19557998327012

src/Components/missions/ASIA-AQ/TLE/sat_tracks.py

@@ -0,0 +1,72 @@
+import matplotlib.pyplot as plt
+
+import numpy as np
+import cartopy.crs as ccrs
+import cartopy.feature as cfeature
+
+from pyobs.tle import TLE
+from datetime  import datetime, timedelta
+
+TLEs = dict( Aqua = 'aqua.tle',
+             NOAA20 = 'noaa20.tle',
+             Sentinel5P = 's5p.tle',
+             SNPP = 'snpp.tle',
+             Terra = 'terra.tle'
+            )
+
+
+def trim(t, x, y, bbox,t_offset):
+
+    I = (x>=bbox[0])&(x<=bbox[1])&(y>=bbox[2])&(y<=bbox[3])
+    t[~I] = np.nan
+    x[~I] = np.nan
+    y[~I] = np.nan
+    ts, xs, ys = [], [], []
+    i = 0
+    t_, x_, y_ = t[I], x[I], y[I]
+    for i in range(len(t_)):
+        ts.append((t_[i]+t_offset).isoformat().split('T')[1][0:5]+'L')
+        xs.append(x_[i])
+        ys.append(y_[i])
+
+
+    return t, x, y, ts, xs, ys
+
+def plot_traj(year,month,day,t_offset):
+
+    fig = plt.figure(figsize=[15, 15])
+    ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree())
+    bbox = [90, 130, -6, 25] # Southeast Asia
+    ax.set_extent(bbox, crs=ccrs.PlateCarree())
+
+    # Put a background image on for nice sea rendering.
+    ax.stock_img()
+
+    ax.add_feature(cfeature.LAND)
+    ax.add_feature(cfeature.COASTLINE)
+
+    dp = 0.2 # delta pixel for text
+    for sat in TLEs:
+        tle = TLE(TLEs[sat])
+        t1, t2, dt = datetime(year,month,day,0), datetime(year,month,day,12), timedelta(seconds=60)
+        t, x, y = tle.getSubpoint(t1,t2,dt)
+        t, x, y, ts, xs, ys = trim(t,x,y,bbox,t_offset)
+        ax.plot(x,y,'-o',label=sat)
+        for i in range(len(ts)):
+           ax.text(xs[i]+dp,ys[i]-dp,ts[i])
+
+
+    ax.legend(loc='lower left')
+    ax.title.set_text('Satellite Tracks for %0d-%02d-%02d'%(year,month,day))            
+    #plt.show()
+
+    fname = 'sattrack_%d-%02d-%02d'%(year,month,day)+'.png'
+    print('Saving',fname)
+    plt.savefig(fname,bbox_inches='tight')
+
+    return
+
+if __name__ == '__main__':
+    t_offset = timedelta(hours = 8) # Manila offset
+    for day in range(5,16):
+        plot_traj(2024,2,day,t_offset)

src/Components/missions/ASIA-AQ/TLE/snpp.tle

@@ -0,0 +1,2 @@
+1 37849U 11061A   24035.82717653  .00000130  00000-0  82355-4 0  9994
+2 37849  98.7288 335.8341 0000823  62.3175 297.8085 14.19551710635949

src/Components/missions/ASIA-AQ/TLE/terra.tle

@@ -0,0 +1,2 @@
+1 25994U 99068A   24035.84026503  .00000599  00000-0  13539-3 0  9990
+2 25994  98.0667 103.1655 0002986  73.9052  51.0127 14.59570635283711