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.

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.

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  29.74’      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  34.99”
   1  f150  150 GHz  30 GHz  0.5  12.25 aW√s  60 uK_RJ√s    104 uK_CMB√s     21”

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).

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: 5.065 km
  seasonal: True
  diurnal: True

2026-04-26 19:48:30.512 INFO: Fetching https://github.com/thomaswmorris/maria-data/raw/master/world_heightmap.h5
Downloading: 100%|██████████| 7.34M/7.34M [00:00<00:00, 118MB/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.

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-04-26 19:48:32.145 INFO: Fetching https://github.com/thomaswmorris/maria-data/raw/master/maps/cluster1.fits
Downloading: 100%|██████████| 4.20M/4.20M [00:00<00:00, 69.0MB/s]

ProjectionMap:
  data(1, 1024, 1024):
    min: -1.923e-04
    max: -2.502e-07
    units: Jy/pixel
    quantity: spectral_flux_density_per_pixel
  nu(1):
    values: [150.] GHz
  y(1024):
    height: 1°
    res: 3.516”
  x(1024):
    width: 1°
    res: 3.516”
  frame: ra/dec
  center:
    ra: 19ʰ24ᵐ37.44ˢ
    dec: -31°13’48.00”
  beam(maj, min, psi): (0 rad, 0 rad, 0 rad)
  memory: 8.389 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.

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.7°, 73.23°)

sim = maria.Simulation(
    instrument,
    plans=plans,
    site=site,
    atmosphere="2d",
    atmosphere_kwargs={"weather": {"pwv": 0.5}},
    map=input_map)

print(sim)

2026-04-26 19:48:49.019 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, 149MB/s]
2026-04-26 19:48:50.234 INFO: Fetching https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/weather/era5/chajnantor.h5
Downloading: 100%|██████████| 4.05M/4.05M [00:00<00:00, 96.6MB/s]

Simulation
├ Instrument(1 array)
│ ├ arrays:
│ │             n     FOV baseline        bands polarized
│ │  array1  1798  29.74’      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  34.99”
│    1  f150  150 GHz  30 GHz  0.5  12.25 aW√s  60 uK_RJ√s    104 uK_CMB√s     21”
├ Site:
│   region: chajnantor
│   timezone: America/Santiago
│   location:
│     longitude: 67°45’17.28” W
│     latitude:  23°01’45.84” S
│     altitude: 5.065 km
│   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.7°, 73.23°)
├ '2d'
└ ProjectionMap:
    data(1, 1, 1, 1024, 1024):
      min: -1.923e-04
      max: -2.502e-07
      units: Jy/pixel
      quantity: spectral_flux_density_per_pixel
    stokes(1):
      components: I
    nu(1):
      values: [150.] GHz
    t(1):
      values: [1.77723291e+09] s
    y(1024):
      height: 1°
      res: 3.516”
    x(1024):
      width: 1°
      res: 3.516”
    frame: ra/dec
    center:
      ra: 19ʰ24ᵐ37.44ˢ
      dec: -31°13’48.00”
    beam(maj, min, psi): (0 rad, 0 rad, 0 rad)
    memory: 8.389 MB

We run the simulation, which spits out a TOD object (which stands for time-ordered data). We can then plot the TOD:

tods = sim.run()

print(tods)
tods[0].plot()

2026-04-26 19:48:52.420 INFO: Simulating observation 1 of 2
Constructing atmosphere: 100%|██████████| 8/8 [00:05<00:00,  1.46it/s]
Generating turbulence: 100%|██████████| 8/8 [00:01<00:00,  7.37it/s]
Sampling turbulence: 100%|██████████| 8/8 [00:10<00:00,  1.33s/it]
Computing atmospheric emission: 100%|██████████| 2/2 [00:02<00:00,  1.32s/it, band=f150]
Sampling map: 100%|██████████| 2/2 [00:08<00:00,  4.25s/it, band=f150, channel=(0 Hz, inf Hz)]
Generating noise: 100%|██████████| 2/2 [00:00<00:00,  3.11it/s, band=f150]
2026-04-26 19:49:29.828 INFO: Simulated observation 1 of 2 in 37.39 s
2026-04-26 19:49:29.828 INFO: Simulating observation 2 of 2
Constructing atmosphere: 100%|██████████| 8/8 [00:04<00:00,  1.61it/s]
Generating turbulence: 100%|██████████| 8/8 [00:01<00:00,  7.27it/s]
Sampling turbulence: 100%|██████████| 8/8 [00:09<00:00,  1.23s/it]
Computing atmospheric emission: 100%|██████████| 2/2 [00:02<00:00,  1.38s/it, band=f150]
Sampling map: 100%|██████████| 2/2 [00:06<00:00,  3.22s/it, band=f150, channel=(0 Hz, inf Hz)]
Generating noise: 100%|██████████| 2/2 [00:00<00:00,  6.66it/s, band=f150]
2026-04-26 19:50:03.541 INFO: Simulated observation 2 of 2 in 33.69 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-04-26T19:49:24.052765+00:00]>, 'altitude': 5065.0, 'region': 'chajnantor', 'pwv': 0.5, 'base_temperature': 272.616}), 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-04-26T19:49:57.681509+00:00]>, 'altitude': 5065.0, 'region': 'chajnantor', 'pwv': 0.911, 'base_temperature': 272.594})]

maria.undebug()

input_map.center

(291.2°, -31.23°)

We can then map the TOD using the built-in mapper.

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-04-26 19:50:10.510 INFO: Inferring center {'ra': '19ʰ24ᵐ37.37ˢ', 'dec': '-31°14’20.14”'} for mapper.
2026-04-26 19:50:10.529 INFO: Inferring mapper width 1.127° for mapper from observation patch.
2026-04-26 19:50:10.530 INFO: Inferring mapper height 1.127° to match supplied width.
2026-04-26 19:50:10.559 INFO: Inferring mapper resolution 31.69” for mapper from observation patch.
2026-04-26 19:50:10.564 INFO: Inferring mapper stokes parameters 'I' for mapper.
Preprocessing TODs: 100%|██████████| 2/2 [00:02<00:00,  1.20s/it]
Preprocessing TODs: 100%|██████████| 2/2 [00:02<00:00,  1.11s/it]
Mapping: 100%|██████████| 4/4 [00:09<00:00,  2.47s/it, tod=1/4]

We can see the recovered map with

output_map.plot(nu_index=[0, 1])