复现arxiv:1812.04739论文Fig.3时NVSS与SUMSS天图匹配异常的代码修正咨询

Let's walk through the key issues in your code and fix them to get the correct NVSS/SUMSS comparison map matching Fig.3 of the paper:

Key Issues Identified #

1. Count array size mismatch: You initialized count1 and count2 with the length of your source lists, but they should be the size of the healpix map (npix) to store counts per pixel. Accessing indices beyond the source list length caused incorrect counts or silent errors.
2. Mismatched sorted indices and fluxes: When you sorted ind1 and ind2, you didn't sort the corresponding flux arrays (fs1, fs2), so your flux filter (>15) was applied to the wrong sources entirely.
3. Inefficient and error-prone counting: Your loop to count overlapping sources used an arbitrary 1000-window check, which isn't reliable for grouping sources by healpix pixel.
4. Incorrect map combination: You overwrote NVSS counts with SUMSS counts instead of computing the ratio/difference that Fig.3 likely shows (a direct comparison between the two surveys).

Corrected Code #

Here's the revised code with all fixes applied, plus optimized vectorized operations:

import numpy as np
import healpy as hp
import matplotlib.pyplot as plt
from matplotlib import cm
import pandas as pd  # Added missing pandas import

nside = 64
npix = hp.nside2npix(nside)

# Process NVSS data
mp1 = pd.read_table('NVSStest.text', delim_whitespace=True, error_bad_lines=False, usecols=range(7))
mp1 = mp1.apply(pd.to_numeric, errors='coerce').dropna()

# Vectorized coordinate conversion (no loops needed!)
RA1 = 15 * (mp1['1'] + mp1['2']/60 + mp1['3']/3600)
DEC1 = mp1['4'] + mp1['5']/60 + mp1['6']/3600
phi1 = np.deg2rad(RA1)
theta1 = np.deg2rad(90 - DEC1)  # Healpy uses colatitude (theta) instead of declination
indices1 = hp.ang2pix(nside, theta1, phi1)

# Count sources with flux >15 per pixel using numpy's efficient bincount
fs1 = mp1['7'].values
mask1 = fs1 > 15
count1 = np.bincount(indices1[mask1], minlength=npix)

# Process SUMSS data
mp2 = pd.read_table('sumsstest.txt', delim_whitespace=True, error_bad_lines=False, usecols=range(11))
mp2 = mp2.apply(pd.to_numeric, errors='coerce').dropna()

# Vectorized coordinate conversion
RA2 = 15 * (mp2['1'] + mp2['2']/60 + mp2['3']/3600)
DEC2 = mp2['4'] + mp2['5']/60 + mp2['6']/3600
phi2 = np.deg2rad(RA2)
theta2 = np.deg2rad(90 - DEC2)
indices2 = hp.ang2pix(nside, theta2, phi2)

# Convert SUMSS flux to NVSS frequency (1400 MHz) using spectral index -0.81
fs2 = mp2['9'].values * (843/1400)**0.81
mask2 = fs2 > 15
count2 = np.bincount(indices2[mask2], minlength=npix)

# Create the comparison map (log ratio to highlight differences, matching Fig.3's likely format)
# Adding 1 to avoid division by zero for pixels with no sources
hpxmap = np.log10((count1 + 1) / (count2 + 1))

# Plot the map
hp.mollview(hpxmap, title="NVSS/SUMSS Source Count Ratio (log10)", norm="hist", cmap=cm.jet)
plt.show()

Additional Notes #

• Vectorization: Replacing loops with numpy vector operations makes the code faster and reduces the chance of off-by-one errors.
• Counting: np.bincount is the standard, reliable way to count occurrences of indices in numpy, ensuring all sources in the same healpix pixel are grouped correctly.
• Map Combination: The log ratio emphasizes regions where one survey detects more sources than the other. If Fig.3 uses a different metric (like total flux density instead of source count), adjust the map calculation to sum fluxes per pixel instead of counting sources.
• Data Validation: Double-check that your input files have the correct column indices (you used mp1['7'] for NVSS flux and mp2['9'] for SUMSS flux—confirm these match the actual data formats from the survey sites).

内容的提问来源于stack exchange，提问作者Krishna