fix units bug

thomaswmorris · Mar 21, 2024 · 2a78f31 · 2a78f31
1 parent b7cbc91
commit 2a78f31
Show file tree

Hide file tree

Showing 6 changed files with 305 additions and 10 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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,