diff --git a/maria/instrument/__init__.py b/maria/instrument/__init__.py index 83e8d79..f3324a7 100644 --- a/maria/instrument/__init__.py +++ b/maria/instrument/__init__.py @@ -84,7 +84,13 @@ def from_config(cls, config): return cls(bands=dets.bands, dets=dets, **config) def __post_init__(self): - self.field_of_view = np.round(self.dets.sky_x.ptp(), 3) + self.field_of_view = np.round(np.degrees(self.dets.sky_x.ptp()), 3) + if self.field_of_view < 0.5 / 60: + self.units = "arcsec" + elif self.field_of_view < 0.5: + self.units = "arcmin" + else: + self.units = "degrees" def __repr__(self): nodef_f_vals = ( @@ -168,7 +174,10 @@ def physical_fwhm(self, z): # """ # return construct_beam_filter(self.physical_fwhm(z), res, beam_profile=beam_profile, buffer=buffer) - def plot_dets(self, units="deg"): + def plot_dets(self, units=None): + if units is None: + units = self.units + fig, ax = plt.subplots(1, 1, figsize=(5, 5), dpi=160) legend_handles = [] diff --git a/maria/instrument/detectors.py b/maria/instrument/detectors.py index 11a1509..bc08910 100644 --- a/maria/instrument/detectors.py +++ b/maria/instrument/detectors.py @@ -173,8 +173,8 @@ def generate_dets( ) band_dets.loc[:, "band_name"] = band - band_dets.loc[:, "sky_x"] = array_offset[0] + detector_offsets[:, 0] - band_dets.loc[:, "sky_y"] = array_offset[1] + detector_offsets[:, 1] + band_dets.loc[:, "sky_x"] = np.radians(array_offset[0] + detector_offsets[:, 0]) + band_dets.loc[:, "sky_y"] = np.radians(array_offset[1] + detector_offsets[:, 1]) band_dets.loc[:, "baseline_x"] = baseline_offset[0] + baselines[:, 0] band_dets.loc[:, "baseline_y"] = baseline_offset[1] + baselines[:, 1] @@ -362,7 +362,7 @@ def n(self): return len(self.df) @property - def offset(self): + def offsets(self): return np.c_[self.sky_x, self.sky_y] @property diff --git a/maria/map/Planck_Parchment_RGB.txt b/maria/map/Planck_Parchment_RGB.txt new file mode 100644 index 0000000..53d036b --- /dev/null +++ b/maria/map/Planck_Parchment_RGB.txt @@ -0,0 +1,256 @@ + 0 0 255 + 0 2 255 + 0 5 255 + 0 8 255 + 0 10 255 + 0 13 255 + 0 16 255 + 0 18 255 + 0 21 255 + 0 24 255 + 0 26 255 + 0 29 255 + 0 32 255 + 0 34 255 + 0 37 255 + 0 40 255 + 0 42 255 + 0 45 255 + 0 48 255 + 0 50 255 + 0 53 255 + 0 56 255 + 0 58 255 + 0 61 255 + 0 64 255 + 0 66 255 + 0 69 255 + 0 72 255 + 0 74 255 + 0 77 255 + 0 80 255 + 0 82 255 + 0 85 255 + 0 88 255 + 0 90 255 + 0 93 255 + 0 96 255 + 0 98 255 + 0 101 255 + 0 104 255 + 0 106 255 + 0 109 255 + 0 112 255 + 0 114 255 + 0 117 255 + 0 119 255 + 0 122 255 + 0 124 255 + 0 127 255 + 0 129 255 + 0 132 255 + 0 134 255 + 0 137 255 + 0 139 255 + 0 142 255 + 0 144 255 + 0 147 255 + 0 150 255 + 0 152 255 + 0 155 255 + 0 157 255 + 0 160 255 + 0 162 255 + 0 165 255 + 0 167 255 + 0 170 255 + 0 172 255 + 0 175 255 + 0 177 255 + 0 180 255 + 0 182 255 + 0 185 255 + 0 188 255 + 0 190 255 + 0 193 255 + 0 195 255 + 0 198 255 + 0 200 255 + 0 203 255 + 0 205 255 + 0 208 255 + 0 210 255 + 0 213 255 + 0 215 255 + 0 218 255 + 0 221 255 + 6 221 254 + 12 221 253 + 18 222 252 + 24 222 251 + 30 222 250 + 36 223 249 + 42 223 248 + 48 224 247 + 54 224 246 + 60 224 245 + 66 225 245 + 72 225 244 + 78 225 243 + 85 226 242 + 91 226 241 + 97 227 240 + 103 227 239 + 109 227 238 + 115 228 237 + 121 228 236 + 127 229 236 + 133 229 235 + 139 229 234 + 145 230 233 + 151 230 232 + 157 230 231 + 163 231 230 + 170 231 229 + 176 232 228 + 182 232 227 + 188 232 226 + 194 233 226 + 200 233 225 + 206 233 224 + 212 234 223 + 218 234 222 + 224 235 221 + 230 235 220 + 236 235 219 + 242 236 218 + 248 236 217 + 255 237 217 + 255 235 211 + 255 234 206 + 255 233 201 + 255 231 196 + 255 230 191 + 255 229 186 + 255 227 181 + 255 226 176 + 255 225 171 + 255 223 166 + 255 222 161 + 255 221 156 + 255 219 151 + 255 218 146 + 255 217 141 + 255 215 136 + 255 214 131 + 255 213 126 + 255 211 121 + 255 210 116 + 255 209 111 + 255 207 105 + 255 206 100 + 255 205 95 + 255 203 90 + 255 202 85 + 255 201 80 + 255 199 75 + 255 198 70 + 255 197 65 + 255 195 60 + 255 194 55 + 255 193 50 + 255 191 45 + 255 190 40 + 255 189 35 + 255 187 30 + 255 186 25 + 255 185 20 + 255 183 15 + 255 182 10 + 255 181 5 + 255 180 0 + 255 177 0 + 255 175 0 + 255 172 0 + 255 170 0 + 255 167 0 + 255 165 0 + 255 162 0 + 255 160 0 + 255 157 0 + 255 155 0 + 255 152 0 + 255 150 0 + 255 147 0 + 255 145 0 + 255 142 0 + 255 140 0 + 255 137 0 + 255 135 0 + 255 132 0 + 255 130 0 + 255 127 0 + 255 125 0 + 255 122 0 + 255 120 0 + 255 117 0 + 255 115 0 + 255 112 0 + 255 110 0 + 255 107 0 + 255 105 0 + 255 102 0 + 255 100 0 + 255 97 0 + 255 95 0 + 255 92 0 + 255 90 0 + 255 87 0 + 255 85 0 + 255 82 0 + 255 80 0 + 255 77 0 + 255 75 0 + 251 73 0 + 247 71 0 + 244 69 0 + 240 68 0 + 236 66 0 + 233 64 0 + 229 62 0 + 226 61 0 + 222 59 0 + 218 57 0 + 215 55 0 + 211 54 0 + 208 52 0 + 204 50 0 + 200 48 0 + 197 47 0 + 193 45 0 + 190 43 0 + 186 41 0 + 182 40 0 + 179 38 0 + 175 36 0 + 172 34 0 + 168 33 0 + 164 31 0 + 161 29 0 + 157 27 0 + 154 26 0 + 150 24 0 + 146 22 0 + 143 20 0 + 139 19 0 + 136 17 0 + 132 15 0 + 128 13 0 + 125 12 0 + 121 10 0 + 118 8 0 + 114 6 0 + 110 5 0 + 107 3 0 + 103 1 0 + 100 0 0 diff --git a/maria/map/__init__.py b/maria/map/__init__.py index b01d8df..bc5aaca 100644 --- a/maria/map/__init__.py +++ b/maria/map/__init__.py @@ -1,18 +1,31 @@ +import os from dataclasses import dataclass, fields from operator import attrgetter from typing import List, Tuple import astropy as ap +import matplotlib as mpl import matplotlib.pyplot as plt import numpy as np import scipy as sp from astropy.io import fits from astropy.wcs import WCS +from matplotlib.colors import ListedColormap from tqdm import tqdm from .. import utils from ..instrument import beams +here, this_filename = os.path.split(__file__) + +# from https://gist.github.com/zonca/6515744 +cmb_cmap = ListedColormap( + np.loadtxt(f"{here}/Planck_Parchment_RGB.txt") / 255.0, name="cmb" +) +cmb_cmap.set_bad("white") + +mpl.colormaps.register(cmb_cmap) + @dataclass class Map: @@ -32,6 +45,7 @@ class Map: frame: str = "ra_dec" units: str = "K_RJ" data: np.array = None # 3D instrument + weight: np.array = None def __post_init__(self): assert len(self.freqs) == self.n_freqs @@ -41,6 +55,9 @@ def __post_init__(self): self.width = self.res * self.n_x self.height = self.res * self.n_y + if self.weight is None: + self.weight = np.ones(self.data.shape) + def __repr__(self): nodef_f_vals = ( (f.name, attrgetter(f.name)(self)) @@ -88,7 +105,9 @@ def y_side(self): y = self.res * np.arange(self.n_y) return y - y.mean() - def plot(self, cmap="plasma", units="degrees", **kwargs): + def plot( + self, cmap="cmb", rel_vmin=0.001, rel_vmax=0.999, units="degrees", **kwargs + ): for i_freq, freq in enumerate(self.freqs): header = fits.header.Header() @@ -135,7 +154,13 @@ def plot(self, cmap="plasma", units="degrees", **kwargs): if units == "arcsec": map_extent = 3600 * np.degrees(map_extent_radians) - vmin, vmax = np.nanpercentile(self.data, q=[1, 99]) + d = self.data.ravel() + w = self.weight.ravel() + + sort = np.argsort(d) + d, w = d[sort], w[sort] + + vmin, vmax = np.interp([rel_vmin, rel_vmax], np.cumsum(w) / np.sum(w), d) map_im = ax.imshow( self.data.T, diff --git a/maria/map/mappers.py b/maria/map/mappers.py index 494187c..8a07604 100644 --- a/maria/map/mappers.py +++ b/maria/map/mappers.py @@ -158,7 +158,8 @@ def run(self): self.raw_map_sums = {band: np.zeros((self.n_x, self.n_y)) for band in self.band} self.raw_map_cnts = {band: np.zeros((self.n_x, self.n_y)) for band in self.band} - self.map_data = np.zeros((len(self.band), self.n_x, self.n_y)) + self.map_data = np.zeros((len(self.band), self.n_y, self.n_x)) + self.map_weight = np.zeros((len(self.band), self.n_y, self.n_x)) for iband, band in enumerate(np.unique(self.band)): self.band_data[band] = {} @@ -258,8 +259,11 @@ def run(self): mask, band_map_denom, 1 ) + self.map_weight[iband] = band_map_denom + self.map = Map( data=self.map_data, + weight=self.map_weight, freqs=np.array([self.band_data[band]["band_center"] for band in self.band]), width=np.degrees(self.width) if self.degrees else self.width, height=np.degrees(self.height) if self.degrees else self.height, diff --git a/maria/sim/base.py b/maria/sim/base.py index 30e6579..a560e48 100644 --- a/maria/sim/base.py +++ b/maria/sim/base.py @@ -1,5 +1,6 @@ import os +import dask.array as da import numpy as np from .. import utils @@ -127,14 +128,14 @@ def __init__( def _run(self): raise NotImplementedError() - def run(self): + def run(self, dtype=np.float32): self.data = {} # Simulate all the junk self._run() tod = TOD( - data=self.tod_data, + data=da.from_array(self.tod_data.astype(dtype)), dets=self.instrument.dets.df, coords=self.det_coords, )