1 Million Star 3D Stellar Density Isosurface Data Visualization Using Python
by matt392 in Circuits > Computers
16 Views, 0 Favorites, 0 Comments
1 Million Star 3D Stellar Density Isosurface Data Visualization Using Python
Image of 1 million stars as a 3D stellar density isosurface sphere data visualization using Python and Claude.ai.
Data downloaded from GAIA DR3 stored on the Flatiron Institute's servers in NYC:
https://flathub.flatironinstitute.org/
Exact star count is 1,032,591. Python code attached and below.
import pandas as pd
import numpy as np
from scipy.ndimage import gaussian_filter
import pyvista as pv
# Load data
df = pd.read_csv('1-million-stars.csv')
# Convert parallax to distance (light years)
df['distance_ly'] = 3262 / df['parallax']
# Remove bad rows
df = df.dropna(subset=['ra', 'dec', 'distance_ly', 'parallax'])
df = df[np.isfinite(df['distance_ly'])]
# ── DISTANCE FILTER ──────────────────────────────────────────
MAX_DISTANCE_LY = 388 # Actual max from data: 388.33 ly
df = df[df['distance_ly'] <= MAX_DISTANCE_LY]
df = df.reset_index(drop=True)
print(f"Stars within {MAX_DISTANCE_LY} ly: {len(df):,}")
# Convert RA/Dec to Cartesian coordinates
ra_rad = np.radians(df['ra'])
dec_rad = np.radians(df['dec'])
x = df['distance_ly'] * np.cos(dec_rad) * np.cos(ra_rad)
y = df['distance_ly'] * np.cos(dec_rad) * np.sin(ra_rad)
z = df['distance_ly'] * np.sin(dec_rad)
# ── 3D GRID ──────────────────────────────────────────────────
grid_size = 150
lim = MAX_DISTANCE_LY + 5
x_range = np.linspace(-lim, lim, grid_size)
y_range = np.linspace(-lim, lim, grid_size)
z_range = np.linspace(-lim, lim, grid_size)
# Calculate density via 3D histogram
H, edges = np.histogramdd(
np.column_stack([x, y, z]),
bins=(x_range, y_range, z_range)
)
# Smooth the density field
H_smooth = gaussian_filter(H, sigma=2.5)
print(f"Density stats: min={H_smooth.min():.4f}, max={H_smooth.max():.4f}, mean={H_smooth.mean():.4f}")
print(f"Non-zero cells: {np.sum(H_smooth > 0):,}/{H_smooth.size:,}")
# ── DIAGNOSTIC: RA/Dec of densest regions ────────────────────
print("\n── Top 20 Densest Grid Cells ──")
x_centers = (x_range[:-1] + x_range[1:]) / 2
y_centers = (y_range[:-1] + y_range[1:]) / 2
z_centers = (z_range[:-1] + z_range[1:]) / 2
flat_indices = np.argsort(H_smooth.ravel())[-20:][::-1]
ix, iy, iz = np.unravel_index(flat_indices, H_smooth.shape)
print(f"{'Rank':<5} {'Density':<10} {'X (ly)':<10} {'Y (ly)':<10} {'Z (ly)':<10} {'RA (deg)':<10} {'Dec (deg)':<10} {'Dist (ly)':<10}")
for rank, (i, j, k) in enumerate(zip(ix, iy, iz), 1):
cx = x_centers[i]
cy = y_centers[j]
cz = z_centers[k]
dist = np.sqrt(cx**2 + cy**2 + cz**2)
ra_deg = np.degrees(np.arctan2(cy, cx)) % 360
dec_deg = np.degrees(np.arcsin(np.clip(cz / dist, -1, 1))) if dist > 0 else 0
density = H_smooth[i, j, k]
print(f"{rank:<5} {density:<10.2f} {cx:<10.1f} {cy:<10.1f} {cz:<10.1f} {ra_deg:<10.1f} {dec_deg:<10.1f} {dist:<10.1f}")
print("\nLook up these RA/Dec coordinates to identify known stellar structures.")
print("Hint: Scorpius-Centaurus ~ RA 240, Dec -40 | Orion OB1 ~ RA 84, Dec 0")
print(" Local Bubble boundary ~ 200-300 ly | Pleiades ~ RA 56, Dec +24\n")
# ── PYVISTA GRID ─────────────────────────────────────────────
grid = pv.ImageData()
grid.dimensions = H_smooth.shape
grid.spacing = (
(x_range.max() - x_range.min()) / (grid_size - 1),
(y_range.max() - y_range.min()) / (grid_size - 1),
(z_range.max() - z_range.min()) / (grid_size - 1)
)
grid.origin = (x_range.min(), y_range.min(), z_range.min())
grid.point_data['density'] = H_smooth.flatten(order='F')
# ── DENSITY THRESHOLDS ───────────────────────────────────────
d_max = H_smooth.max()
thresholds = [
(d_max * 0.005, 'darkblue', 0.15, 'Ultra sparse (0.5%)'),
(d_max * 0.01, 'blue', 0.20, 'Very sparse (1%)'),
(d_max * 0.02, 'cyan', 0.30, 'Sparse (2%)'),
(d_max * 0.05, 'green', 0.45, 'Medium (5%)'),
(d_max * 0.10, 'yellow', 0.55, 'Dense (10%)'),
(d_max * 0.20, 'orange', 0.70, 'Very dense (20%)'),
(d_max * 0.40, 'red', 0.85, 'Maximum (40%)'),
]
print(f"Auto-scaled thresholds (max density = {d_max:.2f}):")
for t in thresholds:
print(f" {t[3]}: {t[0]:.4f}")
# ── PLOTTER ──────────────────────────────────────────────────
plotter = pv.Plotter()
plotter.set_background('black')
surfaces_added = 0
for threshold, color, opacity, label in thresholds:
try:
contour = grid.contour([threshold], scalars='density')
if contour.n_points > 0:
plotter.add_mesh(
contour,
color=color,
opacity=opacity,
show_edges=False,
smooth_shading=True,
label=label
)
surfaces_added += 1
print(f"Added surface at {threshold:.4f}: {contour.n_points:,} points")
else:
print(f"No surface at threshold {threshold:.4f}")
except Exception as e:
print(f"Failed at threshold {threshold:.4f}: {e}")
print(f"\nTotal surfaces added: {surfaces_added}")
# ── ANNOTATIONS ──────────────────────────────────────────────
plotter.add_text(
f'3D Stellar Density Isosurfaces\n1,032,591 Stars within {MAX_DISTANCE_LY} Light Years (Actual Data Extent)',
position='upper_edge',
font_size=12,
color='white'
)
plotter.add_text(
'Thin shells near Earth = Stellar Desert\nThick shells at distance = Dense regions',
position='lower_left',
font_size=9,
color='white'
)
plotter.show_axes()
# Earth marker
origin_sphere = pv.Sphere(radius=8, center=(0, 0, 0))
plotter.add_mesh(origin_sphere, color='white', opacity=0.9)
plotter.add_point_labels(
[[0, 0, 0]],
['Earth'],
font_size=12,
text_color='white',
point_color='white',
point_size=0,
always_visible=True,
shape_opacity=0
)
# ── DISTANCE SCALE RINGS ──────────────────────────────────────
# Flat equatorial circles at 100, 200, 300, 400, 500 ly
ring_distances = [100, 200, 300, 388] # 388 = actual max distance in dataset
for dist in ring_distances:
# Build a circle in the XY plane (galactic equator)
theta = np.linspace(0, 2 * np.pi, 180)
ring_x = dist * np.cos(theta)
ring_y = dist * np.sin(theta)
ring_z = np.zeros_like(theta)
points = np.column_stack([ring_x, ring_y, ring_z])
# Create polyline
n = len(points)
lines = np.zeros((n, 3), dtype=int)
lines[:, 0] = 2
lines[:, 1] = np.arange(n)
lines[:, 2] = np.roll(np.arange(n), -1)
ring = pv.PolyData(points)
ring.lines = lines.flatten()
plotter.add_mesh(ring, color='gray', opacity=0.35, line_width=1)
# Label at positive X edge of each ring
plotter.add_point_labels(
[[dist, 0, 0]],
[f'{dist} ly'],
font_size=9,
text_color='lightgray',
point_color='lightgray',
point_size=0,
always_visible=True,
shape_opacity=0,
)
# ── SCO-CEN LABEL ─────────────────────────────────────────────
# Sco-Cen centroid from diagnostic: RA~269, Dec~-35, dist~350 ly
# Convert to Cartesian for label placement
sco_cen_ra = np.radians(269.0)
sco_cen_dec = np.radians(-35.0)
sco_cen_dist = 410.0 # Push label just outside the 388 ly boundary
sco_cen_x = sco_cen_dist * np.cos(sco_cen_dec) * np.cos(sco_cen_ra)
sco_cen_y = sco_cen_dist * np.cos(sco_cen_dec) * np.sin(sco_cen_ra)
sco_cen_z = sco_cen_dist * np.sin(sco_cen_dec)
plotter.add_point_labels(
[[sco_cen_x, sco_cen_y, sco_cen_z]],
['Scorpius-Centaurus\nOB Association\n(~350-420 ly)'],
font_size=11,
text_color='yellow',
point_color='yellow',
point_size=12,
always_visible=True,
shape_opacity=0
)
# ── COLOR LEGEND ─────────────────────────────────────────────
legend_entries = [
["Ultra sparse (0.5%)", "darkblue"],
["Very sparse (1%)", "blue"],
["Sparse (2%)", "cyan"],
["Medium density (5%)", "green"],
["Dense (10%)", "yellow"],
["Very dense (20%)", "orange"],
["Maximum (40%)", "red"],
]
plotter.add_legend(
legend_entries,
bcolor=(0.1, 0.1, 0.1),
border=True,
loc="center left",
size=(0.18, 0.22),
)
# ── INTERACTIVE WINDOW ───────────────────────────────────────
# Mouse controls:
# Left-click + drag → rotate
# Right-click + drag → zoom
# Middle-click + drag → pan
# Scroll wheel → zoom in/out
# R → reset camera
# Q or close window → quit
plotter.camera_position = 'iso'
plotter.camera.zoom(1.2)
print("\nOpening interactive window — rotate and zoom freely...")
plotter.show()
print(f"\nSummary:")
print(f" Total stars loaded: {len(df):,}")
print(f" Distance cap: {MAX_DISTANCE_LY} ly")
print(f" Grid resolution: {grid_size}^3 = {grid_size**3:,} cells")
print(f" Density range: {H_smooth.min():.2f} to {H_smooth.max():.2f}")
print(f" Surfaces generated: {surfaces_added}")
Downloads
Supplies
- Python
- Pandas
- Numpy
- Ccipy
- Pyvista
- 1 million star dataset from GAIA DR3
1 Million Stars 3D Stellar Density Isosurface Data Visualization Using Python
import pandas as pd
import numpy as np
from scipy.ndimage import gaussian_filter
import pyvista as pv
# Load data
df = pd.read_csv('1-million-stars.csv')
# Convert parallax to distance (light years)
df['distance_ly'] = 3262 / df['parallax']
# Remove bad rows
df = df.dropna(subset=['ra', 'dec', 'distance_ly', 'parallax'])
df = df[np.isfinite(df['distance_ly'])]
# ── DISTANCE FILTER ──────────────────────────────────────────
MAX_DISTANCE_LY = 388 # Actual max from data: 388.33 ly
df = df[df['distance_ly'] <= MAX_DISTANCE_LY]
df = df.reset_index(drop=True)
print(f"Stars within {MAX_DISTANCE_LY} ly: {len(df):,}")
# Convert RA/Dec to Cartesian coordinates
ra_rad = np.radians(df['ra'])
dec_rad = np.radians(df['dec'])
x = df['distance_ly'] * np.cos(dec_rad) * np.cos(ra_rad)
y = df['distance_ly'] * np.cos(dec_rad) * np.sin(ra_rad)
z = df['distance_ly'] * np.sin(dec_rad)
# ── 3D GRID ──────────────────────────────────────────────────
grid_size = 150
lim = MAX_DISTANCE_LY + 5
x_range = np.linspace(-lim, lim, grid_size)
y_range = np.linspace(-lim, lim, grid_size)
z_range = np.linspace(-lim, lim, grid_size)
# Calculate density via 3D histogram
H, edges = np.histogramdd(
np.column_stack([x, y, z]),
bins=(x_range, y_range, z_range)
)
# Smooth the density field
H_smooth = gaussian_filter(H, sigma=2.5)
print(f"Density stats: min={H_smooth.min():.4f}, max={H_smooth.max():.4f}, mean={H_smooth.mean():.4f}")
print(f"Non-zero cells: {np.sum(H_smooth > 0):,}/{H_smooth.size:,}")
# ── DIAGNOSTIC: RA/Dec of densest regions ────────────────────
print("\n── Top 20 Densest Grid Cells ──")
x_centers = (x_range[:-1] + x_range[1:]) / 2
y_centers = (y_range[:-1] + y_range[1:]) / 2
z_centers = (z_range[:-1] + z_range[1:]) / 2
flat_indices = np.argsort(H_smooth.ravel())[-20:][::-1]
ix, iy, iz = np.unravel_index(flat_indices, H_smooth.shape)
print(f"{'Rank':<5} {'Density':<10} {'X (ly)':<10} {'Y (ly)':<10} {'Z (ly)':<10} {'RA (deg)':<10} {'Dec (deg)':<10} {'Dist (ly)':<10}")
for rank, (i, j, k) in enumerate(zip(ix, iy, iz), 1):
cx = x_centers[i]
cy = y_centers[j]
cz = z_centers[k]
dist = np.sqrt(cx**2 + cy**2 + cz**2)
ra_deg = np.degrees(np.arctan2(cy, cx)) % 360
dec_deg = np.degrees(np.arcsin(np.clip(cz / dist, -1, 1))) if dist > 0 else 0
density = H_smooth[i, j, k]
print(f"{rank:<5} {density:<10.2f} {cx:<10.1f} {cy:<10.1f} {cz:<10.1f} {ra_deg:<10.1f} {dec_deg:<10.1f} {dist:<10.1f}")
print("\nLook up these RA/Dec coordinates to identify known stellar structures.")
print("Hint: Scorpius-Centaurus ~ RA 240, Dec -40 | Orion OB1 ~ RA 84, Dec 0")
print(" Local Bubble boundary ~ 200-300 ly | Pleiades ~ RA 56, Dec +24\n")
# ── PYVISTA GRID ─────────────────────────────────────────────
grid = pv.ImageData()
grid.dimensions = H_smooth.shape
grid.spacing = (
(x_range.max() - x_range.min()) / (grid_size - 1),
(y_range.max() - y_range.min()) / (grid_size - 1),
(z_range.max() - z_range.min()) / (grid_size - 1)
)
grid.origin = (x_range.min(), y_range.min(), z_range.min())
grid.point_data['density'] = H_smooth.flatten(order='F')
# ── DENSITY THRESHOLDS ───────────────────────────────────────
d_max = H_smooth.max()
thresholds = [
(d_max * 0.005, 'darkblue', 0.15, 'Ultra sparse (0.5%)'),
(d_max * 0.01, 'blue', 0.20, 'Very sparse (1%)'),
(d_max * 0.02, 'cyan', 0.30, 'Sparse (2%)'),
(d_max * 0.05, 'green', 0.45, 'Medium (5%)'),
(d_max * 0.10, 'yellow', 0.55, 'Dense (10%)'),
(d_max * 0.20, 'orange', 0.70, 'Very dense (20%)'),
(d_max * 0.40, 'red', 0.85, 'Maximum (40%)'),
]
print(f"Auto-scaled thresholds (max density = {d_max:.2f}):")
for t in thresholds:
print(f" {t[3]}: {t[0]:.4f}")
# ── PLOTTER ──────────────────────────────────────────────────
plotter = pv.Plotter()
plotter.set_background('black')
surfaces_added = 0
for threshold, color, opacity, label in thresholds:
try:
contour = grid.contour([threshold], scalars='density')
if contour.n_points > 0:
plotter.add_mesh(
contour,
color=color,
opacity=opacity,
show_edges=False,
smooth_shading=True,
label=label
)
surfaces_added += 1
print(f"Added surface at {threshold:.4f}: {contour.n_points:,} points")
else:
print(f"No surface at threshold {threshold:.4f}")
except Exception as e:
print(f"Failed at threshold {threshold:.4f}: {e}")
print(f"\nTotal surfaces added: {surfaces_added}")
# ── ANNOTATIONS ──────────────────────────────────────────────
plotter.add_text(
f'3D Stellar Density Isosurfaces\n1,032,591 Stars within {MAX_DISTANCE_LY} Light Years (Actual Data Extent)',
position='upper_edge',
font_size=12,
color='white'
)
plotter.add_text(
'Thin shells near Earth = Stellar Desert\nThick shells at distance = Dense regions',
position='lower_left',
font_size=9,
color='white'
)
plotter.show_axes()
# Earth marker
origin_sphere = pv.Sphere(radius=8, center=(0, 0, 0))
plotter.add_mesh(origin_sphere, color='white', opacity=0.9)
plotter.add_point_labels(
[[0, 0, 0]],
['Earth'],
font_size=12,
text_color='white',
point_color='white',
point_size=0,
always_visible=True,
shape_opacity=0
)
# ── DISTANCE SCALE RINGS ──────────────────────────────────────
# Flat equatorial circles at 100, 200, 300, 400, 500 ly
ring_distances = [100, 200, 300, 388] # 388 = actual max distance in dataset
for dist in ring_distances:
# Build a circle in the XY plane (galactic equator)
theta = np.linspace(0, 2 * np.pi, 180)
ring_x = dist * np.cos(theta)
ring_y = dist * np.sin(theta)
ring_z = np.zeros_like(theta)
points = np.column_stack([ring_x, ring_y, ring_z])
# Create polyline
n = len(points)
lines = np.zeros((n, 3), dtype=int)
lines[:, 0] = 2
lines[:, 1] = np.arange(n)
lines[:, 2] = np.roll(np.arange(n), -1)
ring = pv.PolyData(points)
ring.lines = lines.flatten()
plotter.add_mesh(ring, color='gray', opacity=0.35, line_width=1)
# Label at positive X edge of each ring
plotter.add_point_labels(
[[dist, 0, 0]],
[f'{dist} ly'],
font_size=9,
text_color='lightgray',
point_color='lightgray',
point_size=0,
always_visible=True,
shape_opacity=0,
)
# ── SCO-CEN LABEL ─────────────────────────────────────────────
# Sco-Cen centroid from diagnostic: RA~269, Dec~-35, dist~350 ly
# Convert to Cartesian for label placement
sco_cen_ra = np.radians(269.0)
sco_cen_dec = np.radians(-35.0)
sco_cen_dist = 410.0 # Push label just outside the 388 ly boundary
sco_cen_x = sco_cen_dist * np.cos(sco_cen_dec) * np.cos(sco_cen_ra)
sco_cen_y = sco_cen_dist * np.cos(sco_cen_dec) * np.sin(sco_cen_ra)
sco_cen_z = sco_cen_dist * np.sin(sco_cen_dec)
plotter.add_point_labels(
[[sco_cen_x, sco_cen_y, sco_cen_z]],
['Scorpius-Centaurus\nOB Association\n(~350-420 ly)'],
font_size=11,
text_color='yellow',
point_color='yellow',
point_size=12,
always_visible=True,
shape_opacity=0
)
# ── COLOR LEGEND ─────────────────────────────────────────────
legend_entries = [
["Ultra sparse (0.5%)", "darkblue"],
["Very sparse (1%)", "blue"],
["Sparse (2%)", "cyan"],
["Medium density (5%)", "green"],
["Dense (10%)", "yellow"],
["Very dense (20%)", "orange"],
["Maximum (40%)", "red"],
]
plotter.add_legend(
legend_entries,
bcolor=(0.1, 0.1, 0.1),
border=True,
loc="center left",
size=(0.18, 0.22),
)
# ── INTERACTIVE WINDOW ───────────────────────────────────────
# Mouse controls:
# Left-click + drag → rotate
# Right-click + drag → zoom
# Middle-click + drag → pan
# Scroll wheel → zoom in/out
# R → reset camera
# Q or close window → quit
plotter.camera_position = 'iso'
plotter.camera.zoom(1.2)
print("\nOpening interactive window — rotate and zoom freely...")
plotter.show()
print(f"\nSummary:")
print(f" Total stars loaded: {len(df):,}")
print(f" Distance cap: {MAX_DISTANCE_LY} ly")
print(f" Grid resolution: {grid_size}^3 = {grid_size**3:,} cells")
print(f" Density range: {H_smooth.min():.2f} to {H_smooth.max():.2f}")
print(f" Surfaces generated: {surfaces_added}")