import numpy as np
import healpy as hp
import libpysal
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap, BoundaryNorm
from scipy import sparse
from scipy.sparse.csgraph import connected_components
from scipy.sparse.csgraph import shortest_path
from esda.moran import Moran_Local_BV
class AstroPySal:
"""
Spatial analysis tools for astrophysica HEALPix data using PySal.
"""
def __init__(
self,
nside,
keepMask=None,
dropIslands=True,
alpha=0.05,
permutations=999
):
self.nside = nside
self.dropIslands = dropIslands
self.alpha = alpha
self.permutations = permutations
npix = hp.nside2npix(nside)
if keepMask is None:
self.keepMask = np.ones(npix, dtype=bool)
else:
self.keepMask = keepMask.copy()
self.keepIds = None
self.w = None
self.lisaCmap = ListedColormap(
[
"#e0e0e0", # Not significant (Neutral light gray)
"#E50000", # HH: Bright Red
"#FFA500", # LH: Soft Pink/Light Coral
"#0343DF", # LL: Dark Blue
"#029386", # HL: Deep Purple
]
)
self.lisaNorm = BoundaryNorm(
[-0.5, 0.5, 1.5, 2.5, 3.5, 4.5],
self.lisaCmap.N,
)
self.lisaLabels = [
"Not Significant",
"HH (High-High)",
"LH (Low-High)",
"LL (Low-Low)",
"HL (High-Low)",
]
def build_healpix_weights(self):
"""
Build row-standardized spatial weights for the retained HEALPix pixels.
Returns
--------
libpysal.weights.W
Spatial weights object.
"""
npix = hp.nside2npix(self.nside)
keepMask = self.keepMask.copy()
while True:
keepIds = np.where(keepMask)[0]
nKeep = len(keepIds)
idMap = np.full(npix, -1, dtype=np.int64)
idMap[keepIds] = np.arange(nKeep)
allNeighbors = hp.get_all_neighbours(self.nside, keepIds)
rowIndices = np.repeat(np.arange(nKeep), 8)
colIndicesOld = allNeighbors.T.flatten()
validMask = colIndicesOld != -1
rowIndices = rowIndices[validMask]
colIndicesOld = colIndicesOld[validMask]
colIndices = idMap[colIndicesOld]
keepEdge = colIndices != -1
rowIndices = rowIndices[keepEdge]
colIndices = colIndices[keepEdge]
degree = np.zeros(nKeep, dtype=int)
np.add.at(degree, rowIndices, 1)
islandMask = degree == 0
if not self.dropIslands or not np.any(islandMask):
break
keepMask[keepIds[islandMask]]=False
data = np.ones(len(rowIndices), dtype=float)
adjMatrix = sparse.coo_matrix(
(data, (rowIndices, colIndices)),
shape=(nKeep, nKeep),
)
w = libpysal.weights.W.from_sparse(adjMatrix)
w.transform="R"
assert w.n == nKeep, (
f"internal error: w.n={w.n} != nKeep={nKeep}"
)
self.keepMask = keepMask
self.keepIds = keepIds
self.w = w
return w
def compute_local_means(
self,
dataMap,
radiusDeg=5.0,
minValidFraction=0.5,
):
"""
Compute a circular local mean for each HEALPix pixel.
Parameters
----------
dataMap: array-like
full-sky HEALPix map at self.nside. Invalid pixels represented as hp.UNSEEN or nan values.
radiusDeg: float
angular radius of the local window in degrees.
minValidFraction : float
Minimum fraction of pixels in a local window that must contain valid data for a local mean to
be calculated.
Returns:
_________
LocalMean: numpy.ndarray
Full-sky HEALPix array containing the local mean for each pixel.
pixels without sufficient valid data are set to hp.UNSEEN
"""
dataMap = np.asarray(dataMap, dtype=np.float64)
npix = hp.nside2npix(self.nside)
if len(dataMap) != npix:
raise ValueError(
f"dataMap has {len(dataMap)} pixels, "
f"but nside ={self.nside} require {npix} pixels."
)
if radiusDeg <= 0:
raise ValueError("radiusDeg must be greater than zero.")
if not 0 < minValidFraction <= 1:
raise ValueError(
"minValidFraction must be greater than 0 and less than or equal to 1."
)
validMask = (
np.isfinite(dataMap)
& (dataMap != hp.UNSEEN)
)
localMean = np.full(npix, hp.UNSEEN, dtype=np.float64)
radiusRad = np.radians(radiusDeg)
validIds = np.where(validMask)[0]
for pixId in validIds:
centerVec = hp.pix2vec(
self.nside,
pixId,
)
windowPixIds = hp.query_disc(
self.nside,
centerVec,
radiusRad,
inclusive=True
)
windowValidMask = validMask[windowPixIds]
nValid = np.sum(windowValidMask)
nTotal = len(windowPixIds)
validFraction = nValid / nTotal
if validFraction < minValidFraction:
continue
localMean[pixId] = np.mean(
dataMap[windowPixIds[windowValidMask]]
)
return localMean
def local_lisa_classification(
self,
xArray,
yArray,
radiusDeg=5.0,
minValidFraction=0.5
):
"""Classify bivariate HEALPix pixels using locally normalized means.
The central X pixel is classified relative to its local X mean.
The neighboring Y field is classified relative to the local Y mean
at each neighboring pixel.
Quadrant codes:_
1 = HH
2 = LH
3 = LL_
4 = HL
Pixels without sufficient local data are set to hp.UNSEEN.
Returns
-------
dict
Local LISA classification results and intermediate maps.
"""
if self.w is None or self.keepIds is None:_
self.build_healpix_weights()
xArray = np.asarray(xArray, dtype=np.float64)
yArray = np.asarray(yArray, dtype=np.float64)
npix = hp.nside2npix(self.nside)
if len(xArray) != npix:
raise ValueError(
f"xArray has {len(xArray)} pixels, "
f"but nside={self.nside} requires {npix} pixels."
)
if len(yArray) != npix:
raise ValueError(
f"yArray has {len(yArray)} pixels, "
)
localMeanX = self.compute_local_means(
xArray,
radiusDeg=radiusDeg,
minValidFraction=minValidFraction,
)
localMeanY = self.compute_local_means(
yArray,
radiusDeg = radiusDeg,
minValidFraction=minValidFraction,
)
validX = (
np.isfinite(xArray)
& (xArray != hp.UNSEEN)
& np.isfinite(localMeanX)
& (localMeanX != hp.UNSEEN)
)
validY = (
np.isfinite(yArray)
& (yArray != hp.UNSEEN)
& np.isfinite(localMeanY)
& (localmeanY != hp.UNSEEN)
)
validBoth = validX & validY
xResidual = np.full(
npix,
hp.UNSEEN,
dtype=np.float64,
)
yResidual = np.full(
npix,
hp.UNSEEN,
dtype = np.float64,
)
localNeighborY = np.full(
npix,
hp.UNSEEN,
dtype=np.float64,
)
validIds = self.keepids[validBoth[self.keepIds]]
localIndices = np.where(
validBoth[self.keepIds]
)[0]
if len(validIds) == 0:
raise ValueError(
"No pixels remain after applying local validity requirements."
)
localSpares = self.w.sparse[
localIndices
][:, localIndices]
localW = libpysal.weights.W.from_sparse(
localSparse
)
localW.transform = "R"
yValid = yResidual[validIds]
neighborY = localW.sparse @ yValid
localNeighborY[validIds] = np.asarray(
neighborY
).ravel()
localMap = np.full(
npix,
hp.UNSEEN,
dtype=np.float64,
)
xHigh = xResidual > 0
xLow = xResidual < 0
yHigh = localneighborY > 0
yLow = local NeighborY < 0
hhMask = validBoth & xHigh & yHigh
lhMask = validBoth & xLow & yHigh
llMask = validBoth & xLow & yLow
hlMask = validBoth & xHigh & yLow
localMap[hhMask] = 1
localMap[lhMask] = 2
localMap[llMask] = 3
localMap[hlMask] = 4
return {
"localMap": localMap,
"localMeanX": localMeanX,
"localMeanY": localMeanY,
"xResidual": xResidual,
"yResidual": yResidual,
"localNeighborY": localNeighborY,
"localKeepIds": validIds,
"localW": localW,
"radiusDeg": radiusDeg,
"minValidFraction": minValidFraction,
}
def lisa_bv_local_report(
self,
xArray,
yArray,
label,
radiusDeg=2.0,
minValidFraction=0.5,
permutations=999
):
"""Calculate a locally normalized bivariate LISA_statistic.
Each map is locally mean-subtracted before calculating Moran's I.
The local mean is calculated within an angular radius around each pixel. Pixels without
sufficient data are excluded.
"""
xArray = np.asarray(xArray, dtype=np.float64)
yArray = np.asarray(yArray, dtype=np.float64)
npix = hp.nside2npix(self.nside)
if len(xArray) != npix:
raise ValueError(
f"xArray has {len(xArray)} pixels, "
f"but nside = {self.nside} requires {npix} pixels."
)
if len(yArray) != npix:
raise ValueError(
f"yArray has {len(yArray)} pixels, "
f"but nside={self.nside} requires {npix} pixels."
)
localMeanX = self.compute_local_means(
xArray,
radiusDeg=radiusDeg,
minValidFraction=minValidFraction,
)
localMeanY = self.compute_local_means(
yArray,
radiusDeg=radiusDeg,
minValidFraction=minValidFraction,
)
xResidual = np.full(
npix,
hp.UNSEEN,
dtype=np.float64,
)
yResidual = np.full(
npix,
hp.UNSEEN,
dtype=np.float64,
)
validX = (
np.isfinite(xArray)
& (xArray!=hp.UNSEEN)
& (np.isfinite(localMeanX))
& (localMeanX != hp.UNSEEN)
)
validY = (
np.isfinite(yArray)
& (yArray != hp.UNSEEN)
& (np.isfinite(localMeanY))
& (localMeanY != hp.UNSEEN)
)
validBoth = validX & validY
xResidual[validBoth]= (
xArray[validBoth] - localMeanX[validBoth]
)
yResidual[validBoth] = (
yArray[validBoth] - localMeanY[validBoth]
)
localKeepIds = self.keepIds[
validBoth[self.keepIds]
]
localIndices = np.where(
validBoth[self.keepIds]
)[0]
if len(localKeepIds) == 0:
raise ValueError(
"No pixels remain after applying the local validity requirements."
)
localSparse = self.w.sparse[
localIndices
][:,localIndices]
localW = libpysal.weights.W.from_sparse(
localSparse
)
localW.transform="R"
xValid = xResidual[localKeepIds]
yValid = yResidual[localKeepIds]
lisaBV = Moran_Local_BV(
xValid,
yValid,
localW,
permutations=permutations,
)
lisaBV.significanceMethod = "pysal"
lisaBV.localMeanX = localMeanX
lisaBV.localMeanY = localMeanY
lisaBV.xResidual = xResidual
lisaBV.yResidual = yResidual
lisaBV.localKeepIds = localKeepIds
lisaBV.localW = localW
lisaBV.radiusDeg = radiusDeg
lisaBV.minValidFraction = minValidFraction
sig = lisaBV.p_sim < self.alpha
counts = {
q: int(np.sum(sig & (lisaBV.q == q)))
for q in [1, 2, 3, 4]
}
print(
f"{label:20s} "
f"significant={np.sum(sig):6d}/{len(xValid)} "
f"HH={counts[1]:5d} "
f"LH={counts[2]:5d} "
f"LL={counts[3]:5d} "
f"HL={counts[4]:5d}"
)
return (
lisaBV,
localKeepIds,
localMeanX,
localMeanY,
xResidual,
yResidual,
)
def generate_gaussian_surrogates(
self,
dataMap,
nSurrogates=100,
lmax=None,
):
"""
Generate Guassian HEALPix surrogate maps taht preserce the angular power spectrum of the
input map
"""
if lmax is None:
lmax = 3 * self.nside - 1
workingMap = dataMap.copy()
validMask = (
np.isfinite(workingMap)
& (workingMap != hp.UNSEEN)
)
fillValue = np.nanmean(workingMap[validMask])
workingMap[~validMask] = fillValue
cl = hp.anafast(
workingMap,
lmax=lmax,
)
surrogates = []
for i in range(nSurrogates):
surrogate=hp.synfast(
cl,
nside=self.nside,
lmax=lmax,
#verbose=False,
)
surrogate[~validMask] = hp.UNSEEN
surrogates.append(surrogate)
return surrogates
def _run_surrogate_lisa(
self,
xArray,
surrogateMap,
quadrants=(1, 4),
):
"""
Run a single surrogate LISA realization and extract the resulting
significant regions.
Returns
-------
lisa : Moran_Local_BV
LISA object for the surrogate realization.
regions : list
Output from extract_lisa_bv_regions().
"""
xValid = xArray[self.keepIds].astype(np.float64)
surrogateValid = surrogateMap[self.keepIds].astype(np.float64)
lisa = self.lisa_bv_report(
xArray,
surrogateMap,
significanceMethod="pysal",
)
regions = self.extract_lisa_bv_regions(
lisa,
quadrants=quadrants,
alpha=self.alpha,
printRes=False,
)
return lisa, regions
def compare_regions_to_surrogates(
self,
xArray,
realRegions,
surrogateMaps,
quadrants=(1,4),
alpha=None,
printResults=True,
):
"""
Compare real LISA region statistics against surrogate realizations. Serves as a regional
significance test, as opposed to the pixel level significance that the built in pySal function does"""
xValid = xArray[self.keepIds].astype(np.float64)
largestRegions = []
numberRegions = []
meanRegionSizes = []
medianRegionSizes = []
totalPixels = []
largestFractions = []
for surrogate in surrogateMaps:
lisa, regions = self._run_surrogate_lisa(
xArray,
surrogate,
quadrants=quadrants,
)
if len(regions) == 0:
largestRegions.append(0)
numberRegions.append(0)
meanRegionSizes.append(0)
medianRegionSizes.append(0)
totalPixels.append(0)
largestFractions.append(0)
continue
regionSizes = np.array(
[region["size"] for region in regions]
)
largestRegions.append(regionSizes.max())
numberRegions.append(len(regionSizes))
meanRegionSizes.append(regionSizes.mean())
medianRegionSizes.append(np.median(regionSizes))
totalPixels.append(regionSizes.sum())
largestFractions.append(regionSizes.max() / regionSizes.sum())
if len(realRegions) == 0:
realSizes = np.array([0])
else:
realSizes = np.array(
[region["size"] for region in realRegions]
)
realStatistics = {
"largestRegion": realSizes.max(),
"numberRegions": len(realSizes) if len(realRegions) > 0 else 0,
"meanRegionSize": realSizes.mean(),
"medianRegionSize": np.median(realSizes),
"totalPixels": realSizes.sum(),
"largestFraction": realSizes.max() / realSizes.sum()
if realSizes.sum() > 0 else 0,
}
surrogateStatistics = {
"largestRegion": np.array(largestRegions),
"numberRegions": np.array(numberRegions),
"meanRegionSize": np.array(meanRegionSizes),
"medianRegionSize": np.array(medianRegionSizes),
"totalPixels": np.array(totalPixels),
"largestFraction": np.array(largestFractions),
}
results = {}
for statistic in surrogateStatistics:
distribution = surrogateStatistics[statistic]
exceedances = np.sum(
distribution >= realStatistics[statistic]
)
pValue = (
exceedances + 1
) / (
len(distribution) + 1
)
results[statistic] = {
"real": realStatistics[statistic],
"mean": distribution.mean(),
"std": distribution.std(),
"p": pValue,
"distribution": distribution,
}
return results
def compare_quadrant_abundance(
self,
xArray,
yArray,
surrogateMaps,
):
"""
Compare abundance of sig. HH, LH, LL and HL pixels against Gaussian surrogate realizations.
"""
xValid = xArray[self.keepIds].astype(np.float64)
yValid = yArray[self.keepIds].astype(np.float64)
realLisa = Moran_Local_BV(
xValid,
yValid,
self.w,
permutations=self.permutations,
)
realSig = realLisa.p_sim < self.alpha
realCounts = {}
for quadrant in [1, 2, 3, 4]:
realCounts[quadrant] = np.sum(
realSig &
(realLisa.q == quadrant)
)
surrogateCounts = {
1: [],
2: [],
3: [],
4: [],
}
for surrogate in surrogateMaps:
lisa, _ = self._run_surrogate_lisa(
xArray,
surrogate,
)
sig = lisa.p_sim < self.alpha
for quadrant in [1,2,3,4]:
surrogateCounts[quadrant].append(
np.sum(
sig &
(lisa.q == quadrant)
)
)
quadrantNames = {
1: "HH",
2: "LH",
3: "LL",
4: "HL",
}
results = {}
for quadrant in [1,2,3,4]:
distribution = np.asarray(
surrogateCounts[quadrant]
)
exceedances = np.sum(
distribution >= realCounts[quadrant]
)
pValue = (
exceedances + 1
) / (
len(distribution) + 1
)
results[quadrantNames[quadrant]] = {
"real": int(realCounts[quadrant]),
"mean": distribution.mean(),
"std": distribution.std(),
"p": pValue,
"distribution": distribution,
}
return results
def report_quadrant_abundance(
self,
results,
):
"""
Print the quadrant abundance comparison.
"""
print("=" * 60)
print("Quadrant Abundance Test")
print("=" * 60)
for quadrant in ["HH", "LH", "LL", "HL"]:
r = results[quadrant]
print(quadrant)
print(f" Real Count: {r['real']}")
print(f" Surrogate Mean: {r['mean']:.1f}")
print(f" Surrogate Std: {r['std']:.1f}")
print(f" Empirical p: {r['p']:.4f}")
if r["p"] < self.alpha:
print(" Result: Significant")
else:
print(" Result: Not significant")
print("-" * 60)
def compare_region_hierarchy(
self,
xArray,
realRegions,
surrogateMaps,
quadrants=(1,4),
nRanks=10,
):
"""
Compare to the hierarchy of the largest regions against Gaussian surrogate realizations
"""
xValid = xArray[self.keepIds].astype(np.float64)
realSizes = np.sort(
[r["size"] for r in realRegions]
)[::-1]
realHierarchy = np.zeros(nRanks)
realHierarchy[:min(len(realSizes), nRanks)] = \
realSizes[:nRanks]
surrogateHierarchy = np.zeros(
(
len(surrogateMaps),
nRanks,
)
)
for i, surrogate in enumerate(surrogateMaps):
_, regions = self._run_surrogate_lisa(
xArray,
surrogate,
quadrants=quadrants,
)
if len(regions)==0:
continue
sizes = np.sort(
[r["size"] for r in regions]
)[::-1]
surrogateHierarchy[
i,
:min(
len(sizes),
nRanks,
)
] = sizes[:nRanks]
results = []
for rank in range(nRanks):
distribution = surrogateHierarchy[:, rank]
exceedances = np.sum(
distribution >= realHierarchy[rank]
)
pValue = (
exceedances + 1
) / (
len(distribution) + 1
)
results.append(
{
"rank": rank + 1,
"real": realHierarchy[rank],
"mean": distribution.mean(),
"std": distribution.std(),
"p": pValue,
"distribution": distribution,
}
)
return results
def report_region_hierarchy(
self,
results,
):
"""
Print the ranked region comparison.
"""
print("=" * 72)
print("Region Hierarchy Test")
print("=" * 72)
print(
f"{'Rank':>6}"
f"{'Real':>12}"
f"{'Mean':>12}"
f"{'Std':>12}"
f"{'p':>12}"
)
print("-" * 72)
for r in results:
print(
f"{r['rank']:6d}"
f"{r['real']:12.0f}"
f"{r['mean']:12.1f}"
f"{r['std']:12.1f}"
f"{r['p']:12.4f}"
)
def report_surrogate_statistics(
self,
results,
):
"""
Print a summary of surrogate region statistics.
"""
print("=" * 60)
print("Region Significance Test")
print("=" * 60)
statisticNames = {
"largestRegion": "Largest Region",
"numberRegions": "Number of Regions",
"meanRegionSize": "Mean Region Size",
"medianRegionSize": "Median Region Size",
"totalPixels": "Total Significant Pixels",
"largestFraction": "Largest Region Fraction",
}
for key, title in statisticNames.items():
statistic = results[key]
print(title)
print(f" Real: {statistic['real']:.3f}")
print(f" Surrogate Mean: {statistic['mean']:.3f}")
print(f" Surrogate Std: {statistic['std']:.3f}")
print(f" Empirical P: {statistic['p']:.4f}")
if statistic["p"] < self.alpha:
print("Result significant")
else:
print("Result insignificant")
print("-"*60)
def calculate_gaussian_significance(
self,
xArray,
yArray,
nSurrogates=100,
):
"""
Compute LISA significance by comparing the observed local Moran's I values against Gaussian
surrogate realizations of the neighbour-producing field.
"""
xValid = xArray[self.keepIds].astype(np.float64)
yValid = yArray[self.keepIds].astype(np.float64)
realLisa = Moran_Local_BV(
xValid,
yValid,
self.w,
permutations=0,
)
surrogateMaps = self.generate_gaussian_surrogates(
yArray,
nSurrogates=nSurrogates,
)
surrogateI = np.zeros(
(
nSurrogates,
len(xValid),
)
)
for i, surrogate in enumerate(surrogateMaps):
surrogateValid = surrogate[self.keepIds]
surrogateLisa = Moran_Local_BV(
xValid,
surrogateValid,
self.w,
permutations=0,
)
surrogateI[i] = surrogateLisa.Is
pValues = np.mean(
np.abs(surrogateI)
>= np.abs(realLisa.Is),
axis=0,
)
realLisa.p_sim=pValues
realLisa.significanceMethod="gaussian"
return realLisa
def lisa_bv_report(
self,
xArray,
yArray,
label="",
permutations=None,
significanceMethod="pysal",
nSurrogates=100,
):
"""
Runs a bivariate Local Moran's I analysis and reports significant quadrant counts.
Returns
-------
esda.moran.Moran_Local_BV
LISA-BV result object.
"""
if self.w is None or self.keepIds is None:
self.build_healpix_weights()
if permutations is None:
permutations = self.permutations
xValid = xArray[self.keepIds].astype(np.float64)
yValid = yArray[self.keepIds].astype(np.float64)
if significanceMethod == "pysal":
xValid = xArray[self.keepIds].astype(np.float64)
yValid = yArray[self.keepIds].astype(np.float64)
lisaBV = Moran_Local_BV(
xValid,
yValid,
self.w,
permutations=self.permutations,
)
lisaBV.significanceMethod = "pysal"
elif significanceMethod == "gaussian":
lisaBV = self.calculate_gaussian_significance(
xArray,
yArray,
nSurrogates=nSurrogates,
)
else:
raise ValueError(
f"Unknown significance method: {significanceMethod}"
)
sig = lisaBV.p_sim < self.alpha
counts = {
q: int(np.sum(sig & (lisaBV.q == q)))
for q in [1,2,3,4]
}
print(
f"{label:20s} significant={np.sum(sig):6d}/{len(xValid)} "
f"HH={counts[1]:5d} LH={counts[2]:5d} "
f"LL={counts[3]:5d} HL={counts[4]:5d}"
)
return lisaBV
def extract_lisa_bv_regions(
self,
lisaBV,
alpha=None,
quadrants=(1,2,3,4),
printRes=False,
):
"""
Find contiguous spatial clusters o significant LISA-BV pixels, report each region's quadrant composition, and
calculate a spherical-mean center for each region in GAL coords (l,b):
QUADRANT CODES:
1 = HH
2 = LH
3 = LL
4 = HL
"""
if self.w is None or self.keepIds is None:
raise RuntimeError(
"Spatial weights have not been built."
"Call build_healpix_weights() first."
)
if alpha is None:
alpha = self.alpha
quadrantNames = {
1: "HH",
2: "LH",
3: "LL",
4: "HL",
}
sig = (
(lisaBV.p_sim < alpha)
& np.isin(lisaBV.q, quadrants)
)
sigLocalIdx = np.where(sig)[0]
sigQuadCodes = lisaBV.q[sigLocalIdx]
subAdj = self.w.sparse[sigLocalIdx][:, sigLocalIdx]
nComponents, labels = connected_components(
subAdj,
directed=False,
)
regionPixIds = self.keepIds[sigLocalIdx]
results=[]
for regionId in range(nComponents):
mask = labels == regionId
pixIds = regionPixIds[mask]
localIndices = sigLocalIdx[mask]
quadCodes = sigQuadCodes[mask]
vecs = np.array(
hp.pix2vec(
self.nside,
pixIds,
)
).T
meanVec = vecs.mean(axis=0)
meanVec /= np.linalg.norm(meanVec)
theta, phi = hp.vec2ang(meanVec)
centerLat = 90.0 - np.degrees(theta)
centerLon = np.degrees(phi)
if centerLon > 180:
centerLon -= 360
angularDistances = np.degrees(
np.arccos(
np.clip(
vecs @ meanVec,
-1,
1,
)
)
)
counts = {
quadrantNames[q]: int(np.sum(quadCodes == q))
for q in quadrants
}
presentTypes = [
t for t, c in counts.items()
if c > 0
]
if len(presentTypes) == 1:
regionType = presentTypes[0]
else:
dominant = max(counts, key=counts.get)
regionType = (
f"mixed (dominant {dominant}: {counts})"
)
vecs = np.array(
hp.pix2vec(self.nside, pixIds)
).T
meanVec = vecs.mean(axis=0)
meanVec /= np.linalg.norm(meanVec)
theta, phi = hp.vec2ang(meanVec)
centerLat = 90.0 - np.degrees(theta[0])
centerLon = np.degrees(phi[0])
if centerLon > 180:
centerLon -= 360
results.append(
{
"regionId": regionId,
"size": len(pixIds),
"type": regionType,
"composition": counts,
"centerGalLat": centerLat,
"centerGalLon": centerLon,
"pixIds": pixIds,
"centerVec": meanVec,
"angularDistances": angularDistances,
"pixIds": pixIds,
"localIndices": localIndices,
}
)
results.sort(
key=lambda region: -region["size"]
)
if printRes:
for region in results[:15]:
print(
f"region {region['regionId']:3d} "
f"size={region['size']:5d} "
f"type={region['type']:25s}"
f"\ncenter=("
f"l={region['centerGalLon']:7.2f}, "
f"b={region['centerGalLat']:7.2f})"
f"\n"
f"----------------------------------------------------------------"
)
return results
def swap_lisa_bv_quadrants(
self,
lisaMap,
):
"""
Return a copy of a LISA cluster map with the Low-High and
High-Low quadrants exchanged. This is useful when plotting
the reversed variable ordering in a bivariate Moran analysis.
"""
swappedMap = lisaMap.copy()
lhMask = lisaMap == 2
hlMask = lisaMap == 4
swappedMap[lhMask] = 4
swappedMap[hlMask] = 2
return swappedMap
def build_lisa_bv_map(self, lisaBV, alpha=None):
"""
Maps Moran_Local_BV results onto a full-sky HEALPix array.
"""
if alpha is None:
alpha = self.alpha
npix = hp.nside2npix(self.nside)
quadMap = np.full(npix, hp.UNSEEN)
sig = lisaBV.p_sim < alpha
code = np.where(sig, lisaBV.q, 0)
quadMap[self.keepIds] = code
return quadMap
def plot_lisa_bv_map(
self,
lisaMap,
title="LISA Bivariate Cluster Map",
reverse=False,
showLegend=True,
):
"""
Plot a LISA-BV cluster map using the built-in categorical colormap.
"""
if reverse:
lisaMap = self.swap_lisa_bv_quadrants(lisaMap)
hp.newvisufunc.projview(
lisaMap,
cmap=self.lisaCmap,
norm=self.lisaNorm,
min=0,
max=4,
coord=["G"],
graticule=True,
graticule_labels=True,
title=title,
cbar=False,
)
if showLegend:
labels = self.lisaLabels.copy()
if reverse:
labels[2], labels[4] = labels[4], labels[2]
colors = self.lisaCmap.colors
handles = [
plt.Rectangle((0, 0), 1, 1, color=color)
for color in colors
]
plt.legend(
handles,
labels,
loc="lower center",
bbox_to_anchor=(0.5, -0.30),
ncol=3,
frameon=True,
)
def plot_region(
self,
region
):
pixIds = region["pixIds"]
npix = hp.nside2npix(self.nside)
highlightMap = np.full(npix, hp.UNSEEN)
highlightMap[pixIds] = 1
highlightMap[highlightMap < -1e30] = 0
hp.newvisufunc.projview(
highlightMap,
coord=["G"],
graticule=True,
graticule_labels=True,
title=f"Region {region['regionId']} ({len(pixIds)} pixels)",
cbar=False,
cmap='plasma',
min = np.nanpercentile(highlightMap, 0),
max = np.nanpercentile(highlightMap, 100)
)
def describe_region(
self,
regionId,
regionPixIds,
regionLabels,
):
"""Return the HEALPix pixel IDs and galactic coordinates for a region,
along with a printed angular extent summary.
"""
mask = regionLabels == regionId
pixIds = regionPixIds[mask]
theta, phi = hp.pix2ang(
self.nside,
pixIds,
)
galLatRegion = 90.0 - np.degrees(theta)
galLonRegion = np.degrees(phi)
galLonRegion = np.where(
galLonRegion > 180,
galLonRegion - 360,
galLonRegion,
)
print(
f" galactic latitude: "
f"{galLatRegion.min():.1f} to "
f"{galLatRegion.max():.1f} deg "
f"(center ~{galLatRegion.mean():.1f})"
)
return (
pixIds,
galLatRegion,
galLonRegion,
)
print(f"Region {regionId}: {len(pixIds)} pixels")
def compute_region_morphology(
self,
region,
):
"""
Compute basic morphology statistics for a single region.
Paramaters
----------
region: dict
output from extract_lisa_bv__regions()
returns
-------
dict
Morphological properties of the region
"""
pixIds = region["pixIds"]
centerVec = region["centerVec"]
angularDistances = region["angularDistances"]
extentDeg = angularDistances.max()
meanRadius = angularDistances.mean()
radiusStd = angularDistances.std()
pixelArea = hp.nside2pixarea(
self.nside,
degrees=True,
)
regionArea = len(pixIds) * pixelArea
circleArea = np.pi * extentDeg ** 2
filllingFactor = regionArea / circleArea
localIndices = region["localIndices"]
subgraph = self.w.sparse[
localIndices
][:, localIndices]
graphDistances = shortest_path(
subgraph,
directed=False,
unweighted=True,
)
finiteDistances = graphDistances[
np.isfinite(graphDistances)
]
graphDiameter = finiteDistances.max()
meanDegree = np.asarray(
subgraph.sum(axis=1)
).mean()
return{
"extentDeg": extentDeg,
"meanRadius": meanRadius,
"radiusStd": radiusStd,
"fillingFactor": fillingFactor,
"graphDiameter": graphDiameter,
"meanDegree": meanDegree,
}
def compare_region_persistence(
self,
xArray,
yArray,
quadrants=(1,4),
alphas=None,
):
"""
Measure how the detected LISA regoins persist as the PySal significance threshold is varied.
params:
-------
xArray: ndarray
primary data array
yArray: ndarray
neighbour-producing data array
quadrants: tuple, optional
LISA quadrants to include
alphas: sequence, optional
Significance thresholds to test.
defaults to [0.20, 0.10, 0.05, 0.02, 0.01].
Returns
-------
dict
Region statistics as a function of alpha.
"""
if alphas is None:
alphas = [
0.20,
0.10,
0.05,
0.02,
0.01
]
xValid = xArray[self.keepIds].astype(np.float64)
yValid = yArray[self.keepIds].astype(np.float64)
lisa = Moran_Local_BV(
xValid,
yValid,
self.w,
permutations=self.permutations,
)
largestRegion = []
numberRegions = []
totalPixels = []
for alpha in alphas:
regions = self.extract_lisa_bv_regions(
lisa,
quadrants=quadrants,
alpha=alpha,
printRes=False,
)
numberRegions.append(
len(regions)
)
if len(regions) == 0:
largestRegions.append(0)
totalPixels.append(0)
continue
sizes = np.array(
[
region["size"]
for region in regions
]
)
largestRegion.append(
sizes.max()
)
totalPixels.append(
sizes.sum()
)
return{
"alpha": np.asarray(alphas),
"largestRegion":np.asarray(
largestRegion
),
"numberRegions": np.asarray(
numberRegions
),
"totalPixels": np.asarray(
totalPixels
),
}
def report_region_persistence(
self,
persistance,
):
"""
Print a summary of region persistence across significance thresholds.
"""
print("=" * 70)
print("Region persistence")
print("=" * 70)
print(
f"{'Alpha':>8}"
f"{'Largest':>12}"
f"{'Regions':>12}"
f"{'Pixels':>12}"
)
print("-" * 70)
for alpha, largest, nRegions, pixels in zip(
persistence["alpha"],
persistence["largestRegion"],
persistence["numberRegions"],
persistence["totalPixels"],
):
print(
f"{alpha:8.2f}"
f"{largest:12d}"
f"{nRegions:12d}"
f"{pixels:12d}"
)
def plot_region_cutout(
self,
dataMap,
region,
paddingFactor=2.0,
title=None,
unit="",
cmap=None,
):
"""
Plot a gnomonic projection centered on a region's spherical mean position.
"""
centerVec = hp.ang2vec(
region["centerGalLon"],
region["centerGalLat"],
lonlat=True,
)
pixVecs = np.array(
hp.pix2vec(
self.nside,
region["pixIds"],
)
).T
angDist = np.degrees(
np.arccos(
np.clip(
pixVecs @ centerVec,
-1,
1,
)
)
)
radiusDeg = angDist.max() * paddingFactor
resoArcmin = max(
radiusDeg * 60 / 100,
0.5,
)
xsize = int(
2 * radiusDeg * 60 / resoArcmin
)
plotTitle = (
title
or f"Region {region['regionId']} "
f"({region['type']}, {region['size']} px)"
)
hp.gnomview(
dataMap,
rot=(
region["centerGalLon"],
region["centerGalLat"],
0,
),
coord="G",
xsize=xsize,
reso=resoArcmin,
title=plotTitle,
unit=unit,
cmap=cmap,
)
hp.graticule()
def extract_region_cutout_data(
self,
dataDict,
region,
paddingFactor=2.0,
):
"""
Return pixels and data values contained in a region cutout.
"""
centerVec = hp.ang2vec(
region["centerGalLon"],
region["centerGalLat"],
lonlat=True,
)
pixVecs = np.array(
hp.pix2vec(
self.nside,
region["pixIds"],
)
).T
angDist = np.arccos(
np.clip(
pixVecs @ centerVec,
-1,
1,
)
)
radiusRad = angDist.max() * paddingFactor
cutoutPixIds = hp.query_disc(
self.nside,
centerVec,
radiusRad,
inclusive=True,
)
inRegionMask = np.isin(
cutoutPixIds,
region["pixIds"],
)
cutoutData = {
name: arr[cutoutPixIds]
for name, arr in dataDict.items()
}
return {
"pixIds": cutoutPixIds,
"inRegionMask": inRegionMask,
"data": cutoutData,
"radiusDeg": np.degrees(radiusRad),
}
def plot_footprint_vs_region(self, region, title=None):
"""
Overlay the analysis footprint against a region's cutout, to check
whether the region's boundary is tracing the mask edge rather than
genuine bivariate structure.
"""
footprintMap = self.keepMask.astype(float)
plotTitle = title or f"Footprint mask vs Region {region['regionId']} boundary"
self.plot_region_cutout(footprintMap, region, title=plotTitle)
def check_footprint_edge_overlap(self, region):
"""
Reports what fraction of a region's boundary pixels have at least
one masked-out neighbor — i.e., sit directly against the footprint
edge, versus being separated from it by valid, non-significant pixels.
"""
npix = hp.nside2npix(self.nside)
regionMask = np.zeros(npix, dtype=bool)
regionMask[region["pixIds"]] = True
neighbors = hp.get_all_neighbours(self.nside, region["pixIds"]).T # (nPix, 8)
boundaryPix = []
onFootprintEdge = []
for pixId, nbrs in zip(region["pixIds"], neighbors):
nbrs = nbrs[nbrs != -1]
outsideRegion = nbrs[~regionMask[nbrs]] # neighbors not part of this region
if len(outsideRegion) > 0:
boundaryPix.append(pixId)
if np.any(~self.keepMask[outsideRegion]): # any of those neighbors masked out?
onFootprintEdge.append(pixId)
boundaryPix = np.array(boundaryPix)
onFootprintEdge = np.array(onFootprintEdge)
fraction = len(onFootprintEdge) / len(boundaryPix) if len(boundaryPix) > 0 else 0.0
print(f"{len(onFootprintEdge)}/{len(boundaryPix)} boundary pixels have a masked-out neighbor "
f"({fraction:.1%})")
return {"boundaryPix": boundaryPix, "onFootprintEdge": onFootprintEdge, "fraction": fraction}
def _get_gal_lat(self):
"""
Lazily computes and caches galactic latitude per pixel for self.nside.
"""
if not hasattr(self, "_galLat") or self._galLat is None:
npix = hp.nside2npix(self.nside)
theta, _ = hp.pix2ang(self.nside, np.arange(npix))
self._galLat = 90.0 - np.degrees(theta)
return self._galLat
def detrend_by_latitude(self, dataArray, nBins=18):
"""
Removes a latitude trend from a full-sky array using median-binned
galactic latitude shells (matching the 10-deg shell stratification
used elsewhere in this project). Only pixels in self.keepMask
contribute to the trend or receive a residual.
Returns
-------
residual : np.ndarray
Full-sky array, hp.UNSEEN outside self.keepMask.
trend : np.ndarray
Full-sky array of the fitted per-shell median, hp.UNSEEN outside self.keepMask.
"""
if self.keepIds is None:
raise RuntimeError(
"Spatial weights have not been built."
"Call build_healpix_weights() first."
)
galLat = self._get_gal_lat()
npix = hp.nside2npix(self.nside)
workingMap = dataArray.copy().astype(np.float64)
validMask = np.zeros(npix, dtype=bool)
validMask[self.keepIds] = np.isfinite(workingMap[self.keepIds]) & (workingMap[self.keepIds] != hp.UNSEEN)
trend = np.full(npix, hp.UNSEEN)
residual = np.full(npix, hp.UNSEEN)
binEdges = np.linspace(-90, 90, nBins + 1)
for i in range(nBins):
inBin = np.where(
validMask
& (galLat >= binEdges[i])
& (galLat < binEdges[i + 1])
)[0]
if len(inBin) > 0:
binMedian = np.median(workingMap[inBin])
trend[inBin] = binMedian
residual[inBin] = workingMap[inBin] - binMedian
return residual, trend
def generate_latitude_aware_surrogates(self, dataMap, nSurrogates=100, lmax=None, nBins=18):
"""
Generates Gaussian surrogates that preserve the input map's
latitude trend exactly (real trend, not resimulated) while
randomizing the residual structure via power-spectrum-matched
synfast. Use in place of generate_gaussian_surrogates() when
testing for structure beyond a known latitude confound.
"""
if lmax is None:
lmax = 3 * self.nside - 1
residual, trend = self.detrend_by_latitude(dataMap, nBins=nBins)
npix = hp.nside2npix(self.nside)
validMask = residual != hp.UNSEEN
workingResidual = residual.copy()
fillValue = np.mean(workingResidual[validMask])
workingResidual[~validMask] = fillValue
cl = hp.anafast(workingResidual, lmax=lmax)
surrogates = []
for i in range(nSurrogates):
surrogateResidual = hp.synfast(cl, nside=self.nside, lmax=lmax)
surrogate = np.full(npix, hp.UNSEEN)
surrogate[validMask] = surrogateResidual[validMask] + trend[validMask]
surrogates.append(surrogate)
return surrogates