Custom map simulations¶
In this tutorial we will build a simulation from scratch.
We start by defining a Band that will determine our array’s sensitivity to different spectra. We then generate an array by specifying a field of view, which will be populated by evenly-spaced beams of the given band.
[1]:
import maria
from maria.instrument import Band
f090 = Band(
center=90e9, # in Hz
width=20e9, # in Hz
NET_RJ=40e-6, # in K sqrt(s)
knee=1e0, # in Hz
gain_error=5e-2)
f150 = Band(
center=150e9,
width=30e9,
NET_RJ=60e-6,
knee=1e0,
gain_error=5e-2)
We next define an array config, which specifies how detectors will be distributed on the focal plane. In this case, we supply our two bands as the bands argument, which will generate an array of multichroic detectors (for monochroic detectors, we would supply e.g. "bands": [f090]). The resolution of each detector is determined by frequency of the band and the primary_size parameter. The number of detectors is determined by filling up the specified field_of_view with detector
beams, with a relative spacing determined by the beam_spacing parameter.
[2]:
array = {"field_of_view": 0.5,
"beam_spacing": 1.5,
"primary_size": 25,
"bands": [f090, f150]}
instrument = maria.get_instrument(array=array)
print(instrument)
instrument.plot()
Instrument(1 array)
├ arrays:
│ n FOV baseline bands polarized
│ array1 1798 0.4957° 0 m [f090,f150] False
│
└ bands:
name center width η NEP NET_RJ NET_CMB FWHM
0 f090 90 GHz 20 GHz 0.5 5.445 aW√s 40 uK_RJ√s 49.13 uK_CMB√s 0.5832’
1 f150 150 GHz 30 GHz 0.5 12.25 aW√s 60 uK_RJ√s 104 uK_CMB√s 0.3499’
The Site defines the observing location, as well as the weather conditions. maria knows about a bunch of astronomical observing sites (to see them, run print(maria.site_data)); we can instantiate them using the get_site function. We can modify the site by passing kwargs to the get_site function, like changing its altitude (which will affect the vertical profile of different atmospheric parameters).
[3]:
site = maria.get_site("llano_de_chajnantor", altitude=5065)
print(site)
site.plot()
Site:
region: chajnantor
timezone: America/Santiago
location:
longitude: 67°45’17.28” W
latitude: 23°01’45.84” S
altitude: 5065 m
seasonal: True
diurnal: True
2026-01-21 18:10:15.143 INFO: Fetching https://github.com/thomaswmorris/maria-data/raw/master/world_heightmap.h5
Downloading: 100%|████████████████| 7.34M/7.34M [00:00<00:00, 103MB/s]
Here the fetch function downloads a map to the path map_filename, but map_filename can be any .h5 or .fits file of an image that corresponds to the maria map convention (see Maps). Additional kwargs added to the maria.map.load overwrite the metadata of the loaded map.
[4]:
from maria.io import fetch
map_filename = fetch("maps/cluster1.fits")
input_map = maria.map.load(
filename=map_filename,
nu=150e9,
center=(291.156, -31.23))
input_map.data *= 5e1
print(input_map)
input_map.to("K_RJ").plot()
2026-01-21 18:10:17.132 INFO: Fetching https://github.com/thomaswmorris/maria-data/raw/master/maps/cluster1.fits
Downloading: 100%|████████████████| 4.20M/4.20M [00:00<00:00, 26.6MB/s]
ProjectionMap:
shape(nu, y, x): (1, 1024, 1024)
stokes: naive
nu: [150.] GHz
t: naive
z: naive
quantity: spectral_flux_density_per_pixel
units: Jy/pixel
min: -1.923e-04
max: -2.502e-07
rms: 2.900e-05
center:
ra: 19ʰ24ᵐ37.44ˢ
dec: -31°13’48.00”
size(y, x): (1°, 1°)
resolution(y, x): (3.516”, 3.516”)
beam(maj, min, rot): [[0. 0. 0.]] rad
memory: 16.78 MB
We plan an observation using the Planner, which ensures that a given target as seen by a given site will be high enough above the horizon.
[5]:
from maria import Planner
planner = Planner(start_time="2024-08-06T03:00:00",
target=input_map,
site=site,
constraints={"el": (70, 90)})
plans = planner.generate_plans(total_duration=1200, # in seconds
max_chunk_duration=600, # in seconds
scan_pattern="daisy",
scan_options={"radius": input_map.width.deg / 3},
sample_rate=10)
print(plans)
plans[0].plot()
PlanList(2 plans, 1200 s):
start_time duration target(ra,dec) center(az,el)
chunk
0 2024-08-07 01:31:15.000 +00:00 600 s (291.2°, -31.23°) (119.7°, 71.13°)
1 2024-08-07 01:41:52.500 +00:00 600 s (291.2°, -31.23°) (122.6°, 73.23°)
[6]:
sim = maria.Simulation(
instrument,
plans=plans,
site=site,
atmosphere="2d",
atmosphere_kwargs={"weather": {"pwv": 0.5}},
map=input_map)
print(sim)
2026-01-21 18:10:36.862 INFO: Fetching https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/spectra/am/v3/chajnantor.h5
Downloading: 100%|████████████████| 32.9M/32.9M [00:00<00:00, 91.9MB/s]
2026-01-21 18:10:39.205 INFO: Fetching https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/weather/era5/chajnantor.h5
Downloading: 100%|████████████████| 8.00M/8.00M [00:00<00:00, 33.5MB/s]
Simulation
├ Instrument(1 array)
│ ├ arrays:
│ │ n FOV baseline bands polarized
│ │ array1 1798 0.4957° 0 m [f090,f150] False
│ │
│ └ bands:
│ name center width η NEP NET_RJ NET_CMB FWHM
│ 0 f090 90 GHz 20 GHz 0.5 5.445 aW√s 40 uK_RJ√s 49.13 uK_CMB√s 0.5832’
│ 1 f150 150 GHz 30 GHz 0.5 12.25 aW√s 60 uK_RJ√s 104 uK_CMB√s 0.3499’
├ Site:
│ region: chajnantor
│ timezone: America/Santiago
│ location:
│ longitude: 67°45’17.28” W
│ latitude: 23°01’45.84” S
│ altitude: 5065 m
│ seasonal: True
│ diurnal: True
├ PlanList(2 plans, 1200 s):
│ start_time duration target(ra,dec) center(az,el)
│ chunk
│ 0 2024-08-07 01:31:15.000 +00:00 600 s (291.2°, -31.23°) (119.7°, 71.13°)
│ 1 2024-08-07 01:41:52.500 +00:00 600 s (291.2°, -31.23°) (122.6°, 73.23°)
├ '2d'
└ ProjectionMap:
shape(stokes, nu, t, y, x): (1, 1, 1, 1024, 1024)
stokes: I
nu: [150.] GHz
t: [1.76901901e+09]
z: naive
quantity: spectral_flux_density_per_pixel
units: Jy/pixel
min: -1.923e-04
max: -2.502e-07
rms: 2.900e-05
center:
ra: 19ʰ24ᵐ37.44ˢ
dec: -31°13’48.00”
size(y, x): (1°, 1°)
resolution(y, x): (3.516”, 3.516”)
beam(maj, min, rot): [[[[0. 0. 0.]]]] rad
memory: 16.78 MB
We run the simulation, which spits out a TOD object (which stands for time-ordered data). We can then plot the TOD:
[7]:
tods = sim.run()
print(tods)
tods[0].plot()
2026-01-21 18:10:41.169 INFO: Simulating observation 1 of 2
Constructing atmosphere: 100%|████████████████| 8/8 [00:04<00:00, 1.63it/s]
Generating turbulence: 100%|████████████████| 8/8 [00:01<00:00, 7.66it/s]
Sampling turbulence: 100%|████████████████| 8/8 [00:09<00:00, 1.23s/it]
Computing atmospheric emission: 100%|████████████████| 2/2 [00:02<00:00, 1.28s/it, band=f150]
Sampling map: 100%|████████████████| 2/2 [00:06<00:00, 3.32s/it, band=f150, channel=(0 Hz, inf Hz)]
Generating noise: 100%|████████████████| 2/2 [00:00<00:00, 3.44it/s, band=f150]
2026-01-21 18:11:14.397 INFO: Simulated observation 1 of 2 in 33.21 s
2026-01-21 18:11:14.398 INFO: Simulating observation 2 of 2
Constructing atmosphere: 100%|████████████████| 8/8 [00:04<00:00, 1.66it/s]
Generating turbulence: 100%|████████████████| 8/8 [00:01<00:00, 7.64it/s]
Sampling turbulence: 100%|████████████████| 8/8 [00:09<00:00, 1.15s/it]
Computing atmospheric emission: 100%|████████████████| 2/2 [00:02<00:00, 1.16s/it, band=f150]
Sampling map: 100%|████████████████| 2/2 [00:05<00:00, 2.61s/it, band=f150, channel=(0 Hz, inf Hz)]
Generating noise: 100%|████████████████| 2/2 [00:00<00:00, 6.61it/s, band=f150]
2026-01-21 18:11:45.152 INFO: Simulated observation 2 of 2 in 30.74 s
[TOD(shape=(1798, 6000), fields=['atmosphere', 'map', 'noise'], units='K_RJ', start=2024-08-07 01:41:14.899 +00:00, duration=599.9s, sample_rate=10.0Hz, metadata={'atmosphere': True, 'sim_time': <Arrow [2026-01-21T18:11:09.104107+00:00]>, 'altitude': 5065.0, 'region': 'chajnantor', 'pwv': 0.5, 'base_temperature': 272.7}), TOD(shape=(1798, 6000), fields=['atmosphere', 'map', 'noise'], units='K_RJ', start=2024-08-07 01:51:52.399 +00:00, duration=599.9s, sample_rate=10.0Hz, metadata={'atmosphere': True, 'sim_time': <Arrow [2026-01-21T18:11:39.660099+00:00]>, 'altitude': 5065.0, 'region': 'chajnantor', 'pwv': 0.925, 'base_temperature': 272.683})]
[8]:
maria.undebug()
[9]:
input_map.center
[9]:
(291.2°, -31.23°)
We can then map the TOD using the built-in mapper.
[10]:
from maria.mappers import BinMapper
mapper = BinMapper(
# center=input_map.center,
# frame="ra/dec",
# width=input_map.width,
# height=input_map.height,
# resolution=input_map.width / 256,
tod_preprocessing={
"remove_spline": {"knot_spacing": 60, "remove_el_gradient": True},
"remove_modes": {"modes_to_remove": 1},
},
map_postprocessing={
"gaussian_filter": {"sigma": 1},
},
units="mK_RJ",
tods=tods,
)
mapper.add_tods(tods)
output_map = mapper.run()
2026-01-21 18:11:53.177 INFO: Inferring center {'ra': '19ʰ24ᵐ38.02ˢ', 'dec': '-31°14’8.10”'} for mapper.
2026-01-21 18:11:53.196 INFO: Inferring mapper width 1.134° for mapper from observation patch.
2026-01-21 18:11:53.197 INFO: Inferring mapper height 1.134° to match supplied width.
2026-01-21 18:11:53.214 INFO: Inferring mapper resolution 0.5314’ for mapper from observation patch.
2026-01-21 18:11:53.216 INFO: Inferring mapper stokes parameters 'I' for mapper.
Preprocessing TODs: 100%|████████████████| 2/2 [00:02<00:00, 1.19s/it]
Preprocessing TODs: 100%|████████████████| 2/2 [00:02<00:00, 1.12s/it]
Mapping: 100%|██████████| 4/4 [00:05<00:00, 1.39s/it, tod=1/4]
We can see the recovered map with
[11]:
output_map.plot(nu_index=[0, 1])
