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
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
- I can imagine solutions with different computational complexity, e g. brute force vs use of specialized data structures, but the best answer will depend on the scale of your problem. How many satellites in each set are there going to be? Are you working in a small region of the sky, or all around the globe (I.e. periodic boundary conditions become important)? Although it can probably be inferred from your code, please describe the format of the data. – Matt Haberland Commented Mar 20 at 14:57
- @MattHaberland to begin with, a specific region in the sky would suffice (I suppose for example if I took a GEO satellite and wanted interference with the OW constellation, I just need to look at a segment of the sky. Ideally though, I want to be able to build a bigger model of RF interference, which would mean looking on a global scale. That global scale might have to be a summation of segments however... – SPZHunter Commented Mar 20 at 16:49
2 Answers
Reset to default 2If you have a lot of satellites in a small enough region of the sky to ignore some of the challenges of working with spherical data, I'll outline a more efficient solution: convert your coordinates to a spherical coordinate system (ignoring the radius), then use a KD-tree (e.g. scipy.spatial.KDTree) to efficiently find the pairs of satellites that are close together (in an angular sense). Then you can use the techniques in my other answer to find the actual angles between those pairs of satellites.
For instance, at a single time point, suppose you have 10 OneWeb satellites and 5 GPS satellites. After converting to a spherical coordinate system and ignoring the radius:
import numpy as np
from scipy.spatial import KDTree
rng = np.random.default_rng(328924982346534)
x = rng.random(size=(10, 2)) # theta and phi for OneWeb satellites
y = rng.random(size=(5, 2)) # theta and phi for GPS satellites
kdtree = KDTree(x)
kdtree.query(y) # returns "distances" between each GPS satellite and closest OneWeb
kdtree.query_ball_point(y, r=0.3) # returns ONLY the information about pairs with "distance" less than `r`
The "distances" above are not valid angles because they treat the spherical coordinates as if they were regular vectors and take the Euclidean norm of the difference between the coordinates. That's OK. You can use this to get information about satellites that are "close" in an angular sense, and then calculate the actual angles between them.
The advantage is that instead of calculating all pairwise angle, you will only need to calculate angles between pairs of satellites that are sufficiently close in the sky to have a chance of being considered a "conjunction".
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.