基于半球上点的线密度图
Density map of lines based on points on a half sphere
我是运行一个模拟，输出每个时间步长的各种对象的alt和az位置。我可以将它们绘制成线条，但我 运行 遇到了两个问题，试图从中绘制密度图。
似乎密度图是基于点而不是基于线
半球上的面积不等于面积
也许我必须将数据从 alt/az 转换为 x,y？
可以找到数据here
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import math
df = pd.read_csv('file.csv')
print(df)
starlink_name = df.loc[:,'Name']
starlink_alt = df.loc[:,'starlink_alt']
starlink_az = df.loc[:,'starlink_az']
name = starlink_name.values
alt = starlink_alt.values
az = starlink_az.values
print(name)
print(df['Name'].nunique())
df['Date'] = pd.to_datetime(df['Date'])
for name, df_name in df.groupby('Name'):
print(name)
df_grouped = df.groupby('Name')
list_of_names = list(df_grouped.groups)
#########################################################################################
#LinePlot
#########################################################################################
fig = plt.figure()
ax = fig.add_axes([0.1,0.1,0.8,0.8], polar=True)
# ax.invert_yaxis()
ax.set_theta_zero_location('N')
ax.set_rlim(90, 60, 1)
ax.set_yticks(np.arange(0, 91, 15))
ax.set_rlim(bottom=90, top=0)
def legend_without_duplicate_labels(ax):
handles, labels = ax.get_legend_handles_labels()
unique = [(h, l) for i, (h, l) in enumerate(zip(handles, labels)) if l not in labels[:i]]
ax.legend(*zip(*unique),loc= 'upper right')
ax.plot(np.deg2rad(df2['starlink_az']), df2['starlink_alt'], linestyle='solid', marker='.',linewidth=0.5, markersize=0.1 )
legend_without_duplicate_labels(ax)
plt.show()
#########################################################################################
#DensityPlot
#########################################################################################
# define binning
rbins = np.linspace(0,90, 45)
abins = np.linspace(0,2*np.pi, 360)
#calculate histogram
hist, _, _ = np.histogram2d(np.deg2rad(az), alt, bins=(abins, rbins))
A, R = np.meshgrid(abins, rbins)
# plot
fig, ax = plt.subplots(subplot_kw=dict(projection="polar"))
# ax.invert_yaxis()
ax.set_theta_zero_location('N')
ax.set_rlim(90, 0, 1)
# Note: you must set the end of arange to be slightly larger than 90 or it won't include 90
ax.set_yticks(np.arange(0, 91, 15))
ax.set_rlim(bottom=90, top=0)
pc = ax.pcolormesh(A, R, hist.T, cmap="jet")
fig.colorbar(pc)
plt.show()```
好的，我已经完成了大部分。
如评论中所述，我们想要计算线与极坐标中的 bin 的交点。我们可以尝试解决与圆的交点，即对于 2 个点 B 解决 (1 - t) * A + t * B = r² 所有需要的 r。此解决方案适用于另一种方式，计算线与定义 bin 边缘的“辐条”的交点。
我只给你一个完整注释的代码，它应该是不言自明的：
# params: bin sizes
rmax = 90
rbin = rmax / 45
abin = 2 * np.pi / 360
# load data
df = pd.read_csv('9dbd6e7a.csv', parse_dates=['Date'])
df['Satellite'] = df['Name'].str.split('*', 1, expand=True)[0]
# Compute both coordinate systems of all points
points = pd.concat(axis='columns', objs={
# Radians all the way
'angle': df['starlink_az'].transform(np.deg2rad),
# NB. we flip coordinates to have 0 at the center
'dist': rmax - df['starlink_alt'],
})
points['x'] = points['dist'] * points['angle'].transform(np.cos)
points['y'] = points['dist'] * points['angle'].transform(np.sin)
# Compute segments here as pairs of points
segments = pd.concat(axis='columns', objs=[
points.add_suffix('_a'),
# .groupby() is optional, allows to skip segments between consecutive points of different origins
points.add_suffix('_b').groupby(df['Satellite']).shift(-1),
df['Satellite'],
]).dropna()
# Compute segment length and orthogonal vector
segments['dx'] = segments['x_b'] - segments['x_a']
segments['dy'] = segments['y_b'] - segments['y_a']
segments['length'] = (segments['dx'] ** 2 + segments['dy'] ** 2).transform(np.sqrt)
segments[['x_u', 'y_u']] = segments[['dy', 'dx']].mul({'dy': -1, 'dx': 1}).div(segments['length'], axis='index')
# From vector and one point, we get the angle and distance for polar formulation of line
segments['dist_p'] = (segments[['x_u', 'y_u']] * segments[['x_a', 'y_a']].values).sum(axis='columns')
segments[['x_u', 'y_u', 'dist_p']] = segments[['x_u', 'y_u', 'dist_p']].mul(segments['dist_p'].lt(0).map({True: -1, False: 1}), axis='index')
segments['angle_p'] = np.arctan2(segments['y_u'], segments['x_u'])
# True iff the intercept is in the segment, i.e. between A and B on the line
segments['in'] = segments['dist_p'].lt(segments['dist_a']) & segments['dist_p'].lt(segments['dist_b'])
segments['d_angle'] = (segments['angle_b'] - segments['angle_a']).transform(lambda s: s + (s.lt(-np.pi) * 2. - s.gt(np.pi) * 2.) * np.pi)
angle_rem = segments['angle_a'] % abin
# Generate angles of all line intersections with spokes
steps = pd.concat(axis='columns', objs={
'angle_a': segments['angle_a'],
'angle_b': segments['angle_a'] + segments['d_angle'], # Ensure the bound is in the right direction even if we had to wrap
'start': segments['angle_a'] - segments['angle_a'] % abin + segments['d_angle'].gt(0) * abin,
'end': segments['angle_a'] + segments['d_angle'] - segments['angle_b'] % abin + segments['d_angle'].lt(0) * abin, # no rounding needed at the end
'step': segments['d_angle'].transform(np.sign) * abin,
}).agg(lambda s: [s['angle_a'], *np.arange(s['start'], s['end'], s['step']), s['angle_b']], axis='columns')
steps = steps.explode().infer_objects().to_frame('angle')
# Use polar line formulation to compute the distances at each bin edge
steps['dist'] = segments['dist_p'].reindex_like(steps) / (steps['angle'] - segments['angle_p'].reindex_like(steps)).transform(np.cos)
# This is our histogram, and the part that could be greatly improved, taking into account corner cases.
# We have interpolation at every bin edge « laterally » (along rotation) but we may skip some vertically (along distance)
hist = pd.crosstab(steps['dist'] // rbin, steps['angle'] // abin).reindex(index=np.arange(45), columns=np.arange(360)).fillna(0, downcast='infer')
现在可以完善了，但是给出的结果是这样的：
用以下代码作图：
# plot lines and heatmap side by side
fig, (ax, bx) = plt.subplots(figsize=(20, 8), subplot_kw=dict(projection='polar'), ncols=2)
ticks = np.linspace(0, rmax, 7) # 6 ticks + last edge
# Plot unmodified data with 90 close to center and 0 outside, for verification
ax.plot(df['starlink_az'].transform(np.deg2rad), df['starlink_alt'], linestyle='solid', linewidth=.1)
ax.set_rlim(rmax, 0)
ax.set_yticks(ticks)
ax.set_theta_zero_location('N')
# Plot histogram data with 0 close to center and 90 at outside
grid = np.meshgrid(np.arange(0, 2 * np.pi + abin / 2, abin), np.arange(0, 90 + rbin / 2, rbin))
pc = bx.pcolormesh(*grid, hist, cmap="jet")
bx.set_rlim(0, rmax)
bx.set_yticks(ticks)
# Invert tick display as we inverted values: 0 outside, 90 in the center
bx.set_yticklabels(ticks[::-1])
bx.set_theta_zero_location('N')
fig.colorbar(pc)
plt.show()
我是运行一个模拟，输出每个时间步长的各种对象的alt和az位置。我可以将它们绘制成线条，但我 运行 遇到了两个问题，试图从中绘制密度图。
似乎密度图是基于点而不是基于线
半球上的面积不等于面积
也许我必须将数据从 alt/az 转换为 x,y？
可以找到数据here
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import math
df = pd.read_csv('file.csv')
print(df)
starlink_name = df.loc[:,'Name']
starlink_alt = df.loc[:,'starlink_alt']
starlink_az = df.loc[:,'starlink_az']
name = starlink_name.values
alt = starlink_alt.values
az = starlink_az.values
print(name)
print(df['Name'].nunique())
df['Date'] = pd.to_datetime(df['Date'])
for name, df_name in df.groupby('Name'):
print(name)
df_grouped = df.groupby('Name')
list_of_names = list(df_grouped.groups)
#########################################################################################
#LinePlot
#########################################################################################
fig = plt.figure()
ax = fig.add_axes([0.1,0.1,0.8,0.8], polar=True)
# ax.invert_yaxis()
ax.set_theta_zero_location('N')
ax.set_rlim(90, 60, 1)
ax.set_yticks(np.arange(0, 91, 15))
ax.set_rlim(bottom=90, top=0)
def legend_without_duplicate_labels(ax):
handles, labels = ax.get_legend_handles_labels()
unique = [(h, l) for i, (h, l) in enumerate(zip(handles, labels)) if l not in labels[:i]]
ax.legend(*zip(*unique),loc= 'upper right')
ax.plot(np.deg2rad(df2['starlink_az']), df2['starlink_alt'], linestyle='solid', marker='.',linewidth=0.5, markersize=0.1 )
legend_without_duplicate_labels(ax)
plt.show()
#########################################################################################
#DensityPlot
#########################################################################################
# define binning
rbins = np.linspace(0,90, 45)
abins = np.linspace(0,2*np.pi, 360)
#calculate histogram
hist, _, _ = np.histogram2d(np.deg2rad(az), alt, bins=(abins, rbins))
A, R = np.meshgrid(abins, rbins)
# plot
fig, ax = plt.subplots(subplot_kw=dict(projection="polar"))
# ax.invert_yaxis()
ax.set_theta_zero_location('N')
ax.set_rlim(90, 0, 1)
# Note: you must set the end of arange to be slightly larger than 90 or it won't include 90
ax.set_yticks(np.arange(0, 91, 15))
ax.set_rlim(bottom=90, top=0)
pc = ax.pcolormesh(A, R, hist.T, cmap="jet")
fig.colorbar(pc)
plt.show()```
好的，我已经完成了大部分。
如评论中所述，我们想要计算线与极坐标中的 bin 的交点。我们可以尝试解决与圆的交点，即对于 2 个点 B 解决 (1 - t) * A + t * B = r² 所有需要的 r。此解决方案适用于另一种方式，计算线与定义 bin 边缘的“辐条”的交点。
我只给你一个完整注释的代码，它应该是不言自明的：
# params: bin sizes
rmax = 90
rbin = rmax / 45
abin = 2 * np.pi / 360
# load data
df = pd.read_csv('9dbd6e7a.csv', parse_dates=['Date'])
df['Satellite'] = df['Name'].str.split('*', 1, expand=True)[0]
# Compute both coordinate systems of all points
points = pd.concat(axis='columns', objs={
# Radians all the way
'angle': df['starlink_az'].transform(np.deg2rad),
# NB. we flip coordinates to have 0 at the center
'dist': rmax - df['starlink_alt'],
})
points['x'] = points['dist'] * points['angle'].transform(np.cos)
points['y'] = points['dist'] * points['angle'].transform(np.sin)
# Compute segments here as pairs of points
segments = pd.concat(axis='columns', objs=[
points.add_suffix('_a'),
# .groupby() is optional, allows to skip segments between consecutive points of different origins
points.add_suffix('_b').groupby(df['Satellite']).shift(-1),
df['Satellite'],
]).dropna()
# Compute segment length and orthogonal vector
segments['dx'] = segments['x_b'] - segments['x_a']
segments['dy'] = segments['y_b'] - segments['y_a']
segments['length'] = (segments['dx'] ** 2 + segments['dy'] ** 2).transform(np.sqrt)
segments[['x_u', 'y_u']] = segments[['dy', 'dx']].mul({'dy': -1, 'dx': 1}).div(segments['length'], axis='index')
# From vector and one point, we get the angle and distance for polar formulation of line
segments['dist_p'] = (segments[['x_u', 'y_u']] * segments[['x_a', 'y_a']].values).sum(axis='columns')
segments[['x_u', 'y_u', 'dist_p']] = segments[['x_u', 'y_u', 'dist_p']].mul(segments['dist_p'].lt(0).map({True: -1, False: 1}), axis='index')
segments['angle_p'] = np.arctan2(segments['y_u'], segments['x_u'])
# True iff the intercept is in the segment, i.e. between A and B on the line
segments['in'] = segments['dist_p'].lt(segments['dist_a']) & segments['dist_p'].lt(segments['dist_b'])
segments['d_angle'] = (segments['angle_b'] - segments['angle_a']).transform(lambda s: s + (s.lt(-np.pi) * 2. - s.gt(np.pi) * 2.) * np.pi)
angle_rem = segments['angle_a'] % abin
# Generate angles of all line intersections with spokes
steps = pd.concat(axis='columns', objs={
'angle_a': segments['angle_a'],
'angle_b': segments['angle_a'] + segments['d_angle'], # Ensure the bound is in the right direction even if we had to wrap
'start': segments['angle_a'] - segments['angle_a'] % abin + segments['d_angle'].gt(0) * abin,
'end': segments['angle_a'] + segments['d_angle'] - segments['angle_b'] % abin + segments['d_angle'].lt(0) * abin, # no rounding needed at the end
'step': segments['d_angle'].transform(np.sign) * abin,
}).agg(lambda s: [s['angle_a'], *np.arange(s['start'], s['end'], s['step']), s['angle_b']], axis='columns')
steps = steps.explode().infer_objects().to_frame('angle')
# Use polar line formulation to compute the distances at each bin edge
steps['dist'] = segments['dist_p'].reindex_like(steps) / (steps['angle'] - segments['angle_p'].reindex_like(steps)).transform(np.cos)
# This is our histogram, and the part that could be greatly improved, taking into account corner cases.
# We have interpolation at every bin edge « laterally » (along rotation) but we may skip some vertically (along distance)
hist = pd.crosstab(steps['dist'] // rbin, steps['angle'] // abin).reindex(index=np.arange(45), columns=np.arange(360)).fillna(0, downcast='infer')
现在可以完善了，但是给出的结果是这样的：
用以下代码作图：
# plot lines and heatmap side by side
fig, (ax, bx) = plt.subplots(figsize=(20, 8), subplot_kw=dict(projection='polar'), ncols=2)
ticks = np.linspace(0, rmax, 7) # 6 ticks + last edge
# Plot unmodified data with 90 close to center and 0 outside, for verification
ax.plot(df['starlink_az'].transform(np.deg2rad), df['starlink_alt'], linestyle='solid', linewidth=.1)
ax.set_rlim(rmax, 0)
ax.set_yticks(ticks)
ax.set_theta_zero_location('N')
# Plot histogram data with 0 close to center and 90 at outside
grid = np.meshgrid(np.arange(0, 2 * np.pi + abin / 2, abin), np.arange(0, 90 + rbin / 2, rbin))
pc = bx.pcolormesh(*grid, hist, cmap="jet")
bx.set_rlim(0, rmax)
bx.set_yticks(ticks)
# Invert tick display as we inverted values: 0 outside, 90 in the center
bx.set_yticklabels(ticks[::-1])
bx.set_theta_zero_location('N')
fig.colorbar(pc)
plt.show()