Skip to content

Commit

Permalink
Merge pull request #25 from Joshiwavm/finish_tutorials
Browse files Browse the repository at this point in the history 
Notebook works
  • Loading branch information
@thomaswmorris
thomaswmorris authored Aug 25, 2023
2 parents 5222af7 + 164a6e4 commit 832f044
Show file tree
Hide file tree
Showing 6 changed files with 298 additions and 192 deletions.
296 changes: 208 additions & 88 deletions docs/source/tutorials/mock-observations.ipynb
View file

Large diffs are not rendered by default.

58 changes: 31 additions & 27 deletions maria/mappers.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,9 +37,11 @@ class BaseMapper:
def __init__(self, **kwargs):

self.tods = []
self.map_res = kwargs.get("resolution", np.radians(1/60))
self.map_res = kwargs.get("map_res", np.radians(1/60))
self.map_width = kwargs.get("map_width", np.radians(5))
self.map_height = kwargs.get("map_height", np.radians(5))

self.header = fits.header.Header()
# self.header = fits.header.Header()

@property
def maps(self):
Expand Down Expand Up @@ -94,50 +96,52 @@ def get_map_center_lonlat(self):
mean_unit_colat, mean_unit_lon = np.r_[hp.vec2ang(mean_unit_vec)]

return mean_unit_lon, np.pi/2 - mean_unit_colat



def save_maps(self, filepath):

self.header = self.tods[0].header
self.header['comment'] = 'Made Synthetic observations via maria code'
self.header['comment'] = 'Changed resolution and size of the output map'
self.header['comment'] = 'Overwrote resolution and size of the output map'
self.header['CDELT1'] = np.rad2deg(self.map_res)
self.header['CDELT2'] = np.rad2deg(self.map_res)
self.header['CRPIX1'] = self.maps[list(self.maps.keys())[0]].shape[0]
self.header['CRPIX2'] = self.maps[list(self.maps.keys())[0]].shape[1]

self.header['comment'] = 'Changed pointing location of the output map'
self.header['CRPIX1'] = self.maps[list(self.maps.keys())[0]].shape[0]/2
self.header['CRPIX2'] = self.maps[list(self.maps.keys())[0]].shape[1]/2

center_lon, center_lat = self.get_map_center_lonlat
self.header['CRVAL1'] = np.rad2deg(self.tods[0].cntr[0])
self.header['CRVAL2'] = np.rad2deg(self.tods[0].cntr[1])

self.header['CRVAL1'] = center_lon # ra
self.header['CRVAL2'] = center_lat # dec
self.header['CTYPE1'] = 'RA---SIN'
self.header['CUNIT1'] = 'deg '
self.header['CTYPE2'] = 'DEC--SIN'
self.header['CUNIT2'] = 'deg '
self.header['comment'] = 'Overwrote pointing location of the output map'

self.header['comment'] = 'Changed spectral position of the output map'
self.header['comment'] = 'Overwrote spectral position of the output map'
self.header['CTYPE3'] = 'FREQ '
self.header['CUNIT3'] = 'Hz '
self.header['CRPIX3'] = 1.000000000000E+00

self.header['BTYPE'] = 'Intensity'
# if self.map_data["inbright"] == 'Jy/pixel':
# self.header['BUNIT'] = 'Jy/beam '
# else:
self.header['BUNIT'] = 'Kelvin RJ'
if self.tods[0].unit == 'Jy/pixel':
self.header['BUNIT'] = 'Jy/pixel '
else:
self.header['BUNIT'] = 'Kelvin RJ'

for i, key in enumerate(self.maps.keys()):

# what is this?
# what is this? --> Frequency information in the header
self.header['CRVAL3'] = self.nom_freqs[key]
# self.header['CDELT3'] = self.array.dets[i][1]

save_map = self.maps[list(self.maps.keys())[i]]

# if self.map_data["inbright"] == 'Jy/pixel':
# save_map *= utils.KbrightToJyPix(self.header['CRVAL3'],
# self.header['CDELT1'],
# self.header['CDELT2']
# )
fits.writeto( filename = filepath,
if self.tods[0].unit == 'Jy/pixel':
save_map *= utils.KbrightToJyPix(self.header['CRVAL3'],
np.deg2rad(self.header['CDELT1']),
np.deg2rad(self.header['CDELT2'])
)

fits.writeto( filename = filepath.split('.fits')[0] + '_'+list(self.maps.keys())[i]+'.fits',
data = save_map,
header = self.header,
overwrite = True
Expand All @@ -148,12 +152,12 @@ def save_maps(self, filepath):

class BinMapper(BaseMapper):

def __init__(self, map_width=np.radians(5), map_height=np.radians(5), map_res=np.radians(10/60), **kwargs):
def __init__(self, **kwargs):
super().__init__(**kwargs)

self._nmtr = kwargs.get("n_modes_to_remove", 0)
self.x_bins = np.arange(-0.5*map_width, 0.5*map_width, map_res)
self.y_bins = np.arange(-0.5*map_height, 0.5*map_height, map_res)
self.x_bins = np.arange(-0.5*self.map_width, 0.5*self.map_width, self.map_res)
self.y_bins = np.arange(-0.5*self.map_height, 0.5*self.map_height, self.map_res)
self.n_x, self.n_y = len(self.x_bins) - 1, len(self.y_bins) - 1

def run(self):
Expand Down
11 changes: 10 additions & 1 deletion maria/sim.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ def __init__(self,
noise_model=None,
**kwargs):

#_validate_kwargs(kwargs)
_validate_kwargs(kwargs)

if isinstance(array, str):
array = get_array(array, **kwargs)
Expand Down Expand Up @@ -88,6 +88,15 @@ def run(self):
tod.el = self.pointing.el
tod.ra = self.pointing.ra
tod.dec = self.pointing.dec
tod.cntr = self.pointing.scan_center

if self.map_sim is not None:
tod.unit = self.map_sim.input_map.units
tod.header = self.map_sim.input_map.header
else:
tod.unit = 'K'
tod.header = fits.header.Header()


# number of bands are lost here
tod.data = np.zeros((self.array.n_dets, self.pointing.n_time))
Expand Down
67 changes: 17 additions & 50 deletions maria/sky.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -134,69 +134,36 @@ def __init__(self,
self.input_map_file = map_file
hudl = ap.io.fits.open(map_file)

map_res_radians = np.radians(kwargs.get("map_res", 1/1000))
map_center_radians = np.radians(kwargs.get("map_center", (10.5, 4)))

self.input_map = Map(data=hudl[0].data[None],
freqs=np.atleast_1d(kwargs.get("map_freqs", 90e9)),
res=map_res_radians,
center=map_center_radians,
frame=kwargs.get("map_frame", "ra_dec"),
inbright=kwargs.get("inbright", None),
header=hudl[0].header,
units="K")

freqs = []
for ke in self.array.detectors.keys():
freqs.append(self.array.detectors[ke][0])

self.input_map = Map(data = hudl[0].data[None],
header = hudl[0].header,
freqs = np.atleast_1d(kwargs.get("map_freqs", freqs)),
res = np.radians(kwargs.get("map_res", 1/1000)),
center = np.radians(kwargs.get("map_center", (10.5, 4))),
frame = kwargs.get("map_frame", "ra_dec"),
inbright = kwargs.get("map_inbright", None),
units = kwargs.get("map_units", "K")
)

self.map_X, self.map_Y = utils.lonlat_to_xy(self.RA, self.DEC, *self.input_map.center)

# # Prep mapmaker
# # -------------
# self.input_map_file = map_file
# hudl = ap.io.fits.open(self.input_map_file)
# map_image = hudl[0].data

# if map_image.ndim == 2:
# map_image = map_image[None]

# *_, map_nx, map_ny = map_image.shape

# #map_image = map_image.reshape(-1, map_nx, map_ny)

# if len(self.array.dets) != map_image.shape[0] | map_image.shape[0] != 1:
# raise InvalidNBandsError(len(self.array.dets))

# map is a degree wide by default
# DEFAULT_MAP_CENTER = (self.RA.mean(), self.DEC.mean())
# DEFAULT_MAP_RES = 1 / np.maximum(map_nx, map_ny)

# map_center = kwargs.get("map_center", DEFAULT_MAP_CENTER)
# map_res = kwargs.get("map_res", DEFAULT_MAP_RES)
# map_units = kwargs.get("map_units", 'KRJ')
# map_inbright = kwargs.get("map_inbright", None)

self.input_map.header['HISTORY'] = 'History_input_adjustments'
self.input_map.header['comment'] = 'Changed input CDELT1 and CDELT2'
self.input_map.header['CDELT1'] = self.input_map.res
self.input_map.header['CDELT2'] = self.input_map.res
self.input_map.header['comment'] = 'Changed surface brightness units to ' + self.input_map.units
self.input_map.header['comment'] = 'Repositioned the map on the sky'

# self.input_map_data = {"images": map_image,
# "shape": map_image.shape,
# "header": self.input_map.header,
# "center": map_center,
# "res": map_res,
# "inbright": map_inbright,
# "units": map_units,
# }

if self.input_map.inbright is not None:
self.input_map.data *= self.input_map.inbright / np.nanmax(self.input_map.data)
self.input_map.header[""] = "Amplitude is rescaled."
self.input_map.header['comment'] = "Amplitude is rescaled."

if self.input_map.units == 'Jy/pixel':
for i, nu in enumerate(self.input_map.freqs):
self.input_map.data[i] = self.input_map.data[i] / utils.KbrightToJyPix(nu, self.input_map.res, self.input_map.res)
self.input_map.data[i] = self.input_map.data[i] / utils.KbrightToJyPix(nu,
self.input_map.res,
self.input_map.res)

def run(self, **kwargs):

Expand Down
57 changes: 31 additions & 26 deletions maria/tests/test_mock_obs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,36 +8,41 @@ def test_we_observe():

sim = Simulation(

# Mandatory minimal weither settings
# ---------------------
array = 'MUSTANG-2',
pointing = 'DAISY-2deg',
site = 'GBT',

# True sky input
# ---------------------
map_file = "./maps/protocluster.fits", #input files must be a fits file
# map_file = "/Users/jvanmarr/Documents/GitHub/maria/maps/ACT0329_CL0035.097.z_036.00GHz.fits", #input files must be a fits file
map_center = (4, 10.5), # RA & Dec in degree

# Defeault Observational setup
# ----------------------------
integration_time = 600, # seconds
scan_center = (4, 10.5), # degrees
pointing_frame = "ra_dec", # frame
scan_radius = 1, # How large the scanning pattern is in degree

# Additional inputs:
# ----------------------
quantiles = {'column_water_vapor' : 0.5}, # Weather conditions specific for that site
map_units = 'Jy/pixel', # Kelvin Rayleigh Jeans (KRJ, defeault) or Jy/pixel
# map_inbright = -5.37 * 1e3 * 0.000113, # In units of the map_units key
map_res = 0.5 / 1000, # degree, overwrites header information
# Mandatory minimal weither settings
# ---------------------
array = 'MUSTANG-2', # Array type
pointing = 'DAISY-2deg', # Scanning strategy
site = 'GBT', # Site
atm_model = 'linear_angular', # The atmospheric model, set to None if you want a noiseless observation.
# atm_model = None, # The atmospheric model, set to None if you want a noiseless observation.

# True sky input
# ---------------------
map_file = "./maps/protocluster.fits", # Input files must be a fits file.
# map_file can also be set to None if are only interested in the noise
map_center = (4.5, 10), # RA & Dec in degree

# Defeault Observational setup
# ----------------------------
integration_time = 600, # seconds
scan_center = (4.5, 10), # degrees
pointing_frame = "ra_dec", # frame
scan_radius = 0.05, # How large the scanning pattern is in degree

# Additional inputs:
# ----------------------
quantiles = {'column_water_vapor' : 0.5}, # Weather conditions specific for that site
map_units = 'Jy/pixel', # Kelvin Rayleigh Jeans (K, defeault) or Jy/pixel
# map_inbright = -6e-6, # Linearly scale the map to have this peak value.
map_res = 0.1 / 1000, # degree, overwrites header information
)


tod = sim.run()

mapper = mappers.BinMapper(map_res=np.radians(0.1/60))
mapper = mappers.BinMapper(map_height=np.radians(10/60),
map_width=np.radians(10/60),
map_res=np.radians(0.4/1000))
mapper.add_tods(tod)
mapper.run()

Expand Down
1 change: 1 addition & 0 deletions requirements.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,3 +10,4 @@ scipy
tables
tqdm
weathergen >= 2.2.2
cmocean

0 comments on commit 832f044

Please sign in to comment.