|
| 1 | +# Requirements |
| 2 | + |
| 3 | +import warnings |
| 4 | +import numpy as np |
| 5 | +import matplotlib.pyplot as plt |
| 6 | +import matplotlib.colors as mcolors |
| 7 | +from scipy.stats import gaussian_kde |
| 8 | + |
| 9 | +__all__ = ['calculate_density', 'truncate_colormap', 'plot_density_scatter'] |
| 10 | + |
| 11 | +# Density Calculation and Sorting Function |
| 12 | +def calculate_density(x, y, bandwidth=None): |
| 13 | + """ |
| 14 | + Calculate point density for scatter plot data using Gaussian KDE. |
| 15 | + |
| 16 | + Points are sorted by density so that when plotted, denser regions |
| 17 | + appear on top (to emulate TOPCAT plotting behaviour). |
| 18 | +
|
| 19 | + Parameters: |
| 20 | + x : array-like |
| 21 | + X coordinates of the points. |
| 22 | + y : array-like |
| 23 | + Y coordinates of the points. |
| 24 | + bandwidth : float, optional |
| 25 | + Bandwidth for the KDE. If None, it will be estimated automatically. |
| 26 | +
|
| 27 | + Returns: |
| 28 | + x_sorted : array-like |
| 29 | + X coordinates sorted by density. |
| 30 | + y_sorted : array-like |
| 31 | + Y coordinates sorted by density. |
| 32 | + density: array-like |
| 33 | + Density values sorted from least to most dense. |
| 34 | +
|
| 35 | + Raises: |
| 36 | + ValueError |
| 37 | + If x and y have different lengths, if there are fewer than 2 points, or if bandwidth is non-positive. |
| 38 | +
|
| 39 | + Example Usage: |
| 40 | + x = np.random.randn(100) |
| 41 | + y = np.random.randn(100) |
| 42 | + x_sorted, y_sorted, density = calculate_density(x, y) |
| 43 | + """ |
| 44 | + |
| 45 | + x = np.asarray(x) |
| 46 | + y = np.asarray(y) |
| 47 | + |
| 48 | + # Check same length |
| 49 | + if len(x) != len(y): |
| 50 | + raise ValueError(f"x and y must have same length. Got x: {len(x)}, y: {len(y)}") |
| 51 | + |
| 52 | + # Check minimum points for KDE |
| 53 | + if len(x) < 2: |
| 54 | + raise ValueError(f"Need at least 2 points for density calculation. Got {len(x)}") |
| 55 | + |
| 56 | + if bandwidth is not None and bandwidth <= 0: |
| 57 | + raise ValueError(f"bandwidth must be positive. Got {bandwidth}") |
| 58 | + |
| 59 | + # Handle NaN/inf values |
| 60 | + mask = np.isfinite(x) & np.isfinite(y) |
| 61 | + if not mask.all(): |
| 62 | + n_bad = (~mask).sum() |
| 63 | + warnings.warn(f"Removing {n_bad} points with NaN or infinite values") |
| 64 | + x = x[mask] |
| 65 | + y = y[mask] |
| 66 | + |
| 67 | + xy = np.vstack([x, y]) |
| 68 | + z = gaussian_kde(xy, bw_method=bandwidth)(xy) # Evaluate density |
| 69 | + |
| 70 | + idx = z.argsort() |
| 71 | + x_sorted, y_sorted, density = x[idx], y[idx], z[idx] # Sort points by density |
| 72 | + |
| 73 | + return x_sorted, y_sorted, density |
| 74 | + |
| 75 | +def truncate_colormap(cmap, minval=0.4, maxval=0.9, n=256): |
| 76 | + """ |
| 77 | + Truncate a matplotlib colormap to a specified range. |
| 78 | +
|
| 79 | + Used in this script to avoid the light and dark extremes of colormaps |
| 80 | + to better emulate the TOPCAT density style. |
| 81 | + |
| 82 | + Parameters: |
| 83 | + cmap : matplotlib.colors.Colormap OR str |
| 84 | + The original colormap to be truncated. Can be a colormap object or a string name of a colormap. |
| 85 | + minval : float, default=0.4 |
| 86 | + The minimum value of the colormap to include (between 0 and 1). |
| 87 | + maxval : float, default=0.9 |
| 88 | + The maximum value of the colormap to include (between 0 and 1). |
| 89 | + n : int, default=256 |
| 90 | + The number of discrete colors to generate in the truncated colormap. |
| 91 | + |
| 92 | + Returns: |
| 93 | + LinearSegmentedColormap |
| 94 | + The truncated colormap. |
| 95 | +
|
| 96 | + Validates: |
| 97 | + - minval must be less than maxval. |
| 98 | + - minval and maxval must be between 0 and 1. |
| 99 | +
|
| 100 | + Example Usage: |
| 101 | + truncated_cmap = truncate_colormap('viridis', 0.2, 0.8) |
| 102 | + # or |
| 103 | + original_cmap = plt.get_cmap('viridis') |
| 104 | + truncated_cmap = truncate_colormap(original_cmap, 0.2, 0.8) |
| 105 | + """ |
| 106 | + |
| 107 | + # Allow string input (e.g., 'Reds' instead of plt.get_cmap('Reds')) |
| 108 | + if isinstance(cmap, str): |
| 109 | + cmap = plt.get_cmap(cmap) |
| 110 | + |
| 111 | + # Validate range |
| 112 | + if minval >= maxval: |
| 113 | + raise ValueError(f"minval must be less than maxval. Got minval={minval}, maxval={maxval}") |
| 114 | + |
| 115 | + if not (0 <= minval <= 1 and 0 <= maxval <= 1): |
| 116 | + raise ValueError(f"minval and maxval must be between 0 and 1. Got minval={minval}, maxval={maxval}") |
| 117 | + |
| 118 | + new_cmap = mcolors.LinearSegmentedColormap.from_list( |
| 119 | + 'truncated({n},{a:.2f},{b:.2f})'.format(n=cmap.name, a=minval, b=maxval), |
| 120 | + cmap(np.linspace(minval, maxval, n)) |
| 121 | + ) |
| 122 | + return new_cmap |
| 123 | + |
| 124 | +def plot_density_scatter(x, y, bandwidth=None, cmap='Reds', |
| 125 | + minval=0.4, maxval=0.9, n=256, ax=None, **scatter_kwargs): |
| 126 | + """ |
| 127 | + Create a density scatter plot similar to TOPCAT's density plot style. |
| 128 | +
|
| 129 | + This is a convenience function combining calculate_density() and |
| 130 | + truncate_colormap(). For more control, use those functions separately. |
| 131 | + By default, edgecolor is set to 'none' to match TOPCAT style. |
| 132 | +
|
| 133 | + Parameters: |
| 134 | + x : array-like |
| 135 | + X coordinates of the points. |
| 136 | + y : array-like |
| 137 | + Y coordinates of the points. |
| 138 | + bandwidth : float, optional |
| 139 | + Bandwidth for the KDE. If None, it will be estimated automatically. |
| 140 | + cmap : str or matplotlib.colors.Colormap, default='Reds' |
| 141 | + Colormap to use for density coloring. |
| 142 | + minval : float, default=0.4 |
| 143 | + Minimum value for truncating the colormap. |
| 144 | + maxval : float, default=0.9 |
| 145 | + Maximum value for truncating the colormap. |
| 146 | + n : int, default=256 |
| 147 | + The number of discrete colors to generate in the truncated colormap. |
| 148 | + ax : matplotlib.axes.Axes, optional |
| 149 | + Axes object to plot on. If None, uses current axes. |
| 150 | + **scatter_kwargs : keyword arguments |
| 151 | + Additional keyword arguments passed to plt.scatter(). |
| 152 | + |
| 153 | + Returns: |
| 154 | + scatter : matplotlib.collections.PathCollection |
| 155 | + The scatter plot object. |
| 156 | +
|
| 157 | + Example Usage: |
| 158 | + x = np.random.randn(1000) |
| 159 | + y = np.random.randn(1000) |
| 160 | + cmap='viridis' |
| 161 | + plot_density_scatter(x, y, cmap=cmap, minval=0.2, maxval=0.8) |
| 162 | + """ |
| 163 | + |
| 164 | + x_sorted, y_sorted, density = calculate_density(x, y, bandwidth) |
| 165 | + truncated_cmap = truncate_colormap(cmap, minval, maxval, n) |
| 166 | + |
| 167 | + if ax is None: |
| 168 | + ax = plt.gca() |
| 169 | + |
| 170 | + # Set default edgecolor if not specified by user |
| 171 | + scatter_kwargs.setdefault('edgecolor', 'None') |
| 172 | + |
| 173 | + scatter = ax.scatter(x_sorted, y_sorted, c=density, |
| 174 | + cmap=truncated_cmap, **scatter_kwargs) |
| 175 | + return scatter |
0 commit comments