Original Post: How to populate a 2-d numpy array with values from a third dimension?
New Post:
I'm trying to analyze interference between satellites using numpy and sgp4 python libraries, and want to start with a very basic model which simply shows geometric conjunctions between objects. In this case a conjunction is defined as the angle between two spacecraft's position vectors being equal to 0 degrees (or close to).
To keep it simple, my example shows my current code for one GPS satellite and 10 OneWeb ones, at two points during the day. The GP data needed to load these is provided at the bottom.
Question: I don't know how to do this efficiently - and using loops will get very slow very quickly as I start to introduce more timesteps and more satellites. How can I do this in a vectorized manner?
Example:
The following code gets the position vectors for each spacecraft and stores them in owR and gpsR:
# Imports
import numpy as np
from datetime import datetime as dt
from datetime import timedelta as td
from sgp4 import omm
from sgp4.api import jday, Satrec, SatrecArray
# Set up analysis period
startDate = dt(2025, 3, 17, 0)
endDate = startDate + td(days=1)
step = td(days=0.5)
dateArray = np.arange(startDate, endDate, step).astype(dt)
dateArray = np.array([np.array(jday(x.year, x.month, x.day, x.hour, x.minute, x.second)) for x in dateArray])
jdArray = np.ascontiguousarray(dateArray[:,0]) # Integer + 0.5 Julian day
frArray = np.ascontiguousarray(dateArray[:,1]) # Fraction past 1200UTC of Julian day
# Loading OMM data from CSV
owSats = {}
gpsSats = {}
with open('GP_DATA.csv') as f:
for fields in omm.parse_csv(f):
sat = Satrec()
omm.initialize(sat, fields)
if 'ONEWEB' in fields['OBJECT_NAME']:
owSats[fields['OBJECT_NAME']] = sat
elif 'NAVSTAR' in fields['OBJECT_NAME']:
gpsSats[fields['OBJECT_NAME']] = sat
# Arrays of Satrec objects for SGP4
owSatsArray = SatrecArray(list(owSats.values()))
gpsSatsArray = SatrecArray(list(gpsSats.values()))
# Calculate Position & Velocity for each satellite at each timestep, returns error codes too
owErr, owR, owV = owSatsArray.sgp4(jdArray, frArray)
gpsErr, gpsR, gpsV = gpsSatsArray.sgp4(jdArray, frArray)
The shape of the position vector arrays is (m, n, r), which is (10, 2, 3) and (1, 2, 3) for the OneWeb and GPS arrays respectively.
I then define a function, angle_of_separation(u, v), which takes two vectors and calculates the angle between them using np.dot:
def angle_of_separation(u, v):
cosTheta = np.dot(u, v) / (np.linalg.norm(u)*np.linalg.norm(v))
thetaRad = np.arccos(cosTheta)
thetaDeg = np.degrees(thetaRad)
return thetaRad, thetaDeg
Then, the bit I'm trying to improve, I loop through and find the angle between the GPS satellite and all the OneWeb satellites at all timesteps:
m, n, r = owR.shape
res = np.ones(shape=[m,n])
for i in range(m):
for j in range(n):
rad, deg = angle_of_separation(gpsR[j,:], owR[i,j,:])
res[i,j] = deg
The shape of deg is (i, j), so (10, 2) in this instance.
GP_DATA.csv
NAVSTAR 43 (USA 132),1997-035A,2025-03-18T23:55:30.358560,2.00562658,.0087553,55.7688,116.3810,53.8413,307.0378,0,U,24876,999,20282,0,.19E-6,0
ONEWEB-0012,2019-010A,2025-03-19T17:10:05.865888,13.16596004,.0002263,87.9081,319.3992,91.2044,268.9346,0,U,44057,999,29176,-.65808E-3,-.239E-5,0
ONEWEB-0010,2019-010B,2025-03-19T15:57:09.992736,13.16593707,.0001911,87.9085,319.3913,67.6783,292.4550,0,U,44058,999,29180,.17816E-3,.81E-6,0
ONEWEB-0008,2019-010C,2025-03-19T16:33:37.353888,13.16594232,.0001519,87.9084,319.3871,85.1200,275.0104,0,U,44059,999,29192,-.95577E-3,-.353E-5,0
ONEWEB-0007,2019-010D,2025-03-19T19:41:26.459232,13.15550961,.0002162,87.9140,349.7394,76.2670,283.8701,0,U,44060,999,29163,.44111E-3,.179E-5,0
ONEWEB-0006,2019-010E,2025-03-19T19:47:07.657152,13.15551703,.0002,87.9131,349.6962,72.9243,287.2106,0,U,44061,999,29168,.3057E-3,.128E-5,0
ONEWEB-0011,2019-010F,2025-03-19T19:03:48.012480,13.15548960,.0001684,87.9131,349.7085,71.9060,288.2254,0,U,44062,999,29167,-.43781E-4,-.4E-7,0
ONEWEB-0013,2020-008A,2025-03-19T16:06:50.603328,13.09419131,.0002857,87.8620,144.4872,114.4762,245.6663,0,U,45131,999,24658,-.20908E-2,-.707E-5,0
ONEWEB-0017,2020-008B,2025-03-19T16:08:33.371808,13.10369647,.0001894,87.8836,126.8843,89.6849,270.4495,0,U,45132,999,24842,.51902E-3,.194E-5,0
ONEWEB-0020,2020-008C,2025-03-19T15:25:31.641312,13.10377069,.0001741,87.8844,126.8744,100.3520,259.7803,0,U,45133,999,24881,.76903E-3,.281E-5,0
ONEWEB-0021,2020-008D,2025-03-19T15:47:03.423264,13.10373043,.0002105,87.8843,126.8999,96.8444,263.2922,0,U,45134,999,24912,.18627E-2,.663E-5,0
Here is the brute force solution:
import numpy as np
rng = np.random.default_rng(3289245982346534) # for generating random data
def angle(v1, v2, axis=-1):
# axis=-1, the default, will treat the last axis as the x, y, z
# coordinates.
dot = np.vecdot(v1, v2, axis=axis)
norm = np.linalg.norm(v1, axis=axis) * np.linalg.norm(v2, axis=axis)
return np.arccos(dot / norm)
v1 = rng.random((10, 2, 3))
v2 = rng.random((1, 2, 3))
angle(v1, v2)
This returns an array of shape (10, 2), the angle between the one GPS satellite and 10 OneWeb satellites.
It sounds like may have more than one GPS satellite. In this case:
v2 = rng.random((5, 2, 3))
v2 = v2[:, np.newaxis, ...]
angle(v1, v2)
We insert another dimension with np.newaxis because we want the axis that indexes the GPS satellites to be orthogonal to the axis that indexes the OneWeb satellites. This ensures that all pairs of calculations are performed.
This will return an array of shape (5, 10, 2), the pairwise angles between five GPS satellites and 10 OneWebs.