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

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

Downloading https://github.com/thomaswmorris/maria-data/raw/master/world_heightmap.h5: 100%|████████████████| 7.34M/7.34M [00:00<00:00, 81.6MB/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()

Downloading https://github.com/thomaswmorris/maria-data/raw/master/maps/cluster1.fits: 100%|████████████████| 4.20M/4.20M [00:00<00:00, 13.7MB/s]

ProjectedMap:
  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
  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 rad, 0 rad, 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.

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=50)

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

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

print(sim)

Downloading https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/spectra/am/v2/chajnantor.h5: 100%|████████████████| 14.3M/14.3M [00:00<00:00, 31.5MB/s]
Downloading https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/weather/era5/chajnantor.h5: 100%|████████████████| 8.00M/8.00M [00:00<00:00, 22.6MB/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.22°)
├ '2d'
└ ProjectedMap:
    shape(stokes, nu, y, x): (1, 1, 1024, 1024)
    stokes: I
    nu: [150.] GHz
    t: naive
    z: naive
    quantity: rayleigh_jeans_temperature
    units: K_RJ
      min: -9.574e-04
      max: -1.246e-06
    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 rad, 0 rad, 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:

tods = sim.run()

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

2025-08-24 18:57:51.622 INFO: Simulating observation 1 of 2
Constructing atmosphere: 100%|████████████████| 10/10 [00:15<00:00,  1.59s/it]
Generating turbulence: 100%|████████████████| 10/10 [00:00<00:00, 16.92it/s]
Sampling turbulence: 100%|████████████████| 10/10 [00:10<00:00,  1.00s/it]
Computing atmospheric emission: 100%|████████████████| 2/2 [00:02<00:00,  1.23s/it, band=f150]
Sampling map: 100%|████████████████| 2/2 [00:19<00:00,  9.70s/it, band=f150, channel=(0 Hz, inf Hz), stokes=I]
Generating noise: 100%|████████████████| 2/2 [00:01<00:00,  1.16it/s, band=f150]
2025-08-24 18:59:14.527 INFO: Simulated observation 1 of 2 in 82.89 s
2025-08-24 18:59:14.528 INFO: Simulating observation 2 of 2
Constructing atmosphere: 100%|████████████████| 10/10 [00:16<00:00,  1.69s/it]
Generating turbulence: 100%|████████████████| 10/10 [00:00<00:00, 16.98it/s]
Sampling turbulence: 100%|████████████████| 10/10 [00:09<00:00,  1.02it/s]
Computing atmospheric emission: 100%|████████████████| 2/2 [00:02<00:00,  1.20s/it, band=f150]
Sampling map: 100%|████████████████| 2/2 [00:20<00:00, 10.31s/it, band=f150, channel=(0 Hz, inf Hz), stokes=I]
Generating noise: 100%|████████████████| 2/2 [00:01<00:00,  1.27it/s, band=f150]
2025-08-24 19:00:41.392 INFO: Simulated observation 2 of 2 in 86.85 s

[TOD(shape=(1798, 30000), fields=['atmosphere', 'map', 'noise'], units='K_RJ', start=2024-08-07 01:41:14.979 +00:00, duration=600.0s, sample_rate=50.0Hz, metadata={'atmosphere': True, 'sim_time': <Arrow [2025-08-24T18:58:45.649658+00:00]>, 'altitude': 5065.0, 'region': 'chajnantor', 'pwv': 0.5, 'base_temperature': 272.7}), TOD(shape=(1798, 30000), fields=['atmosphere', 'map', 'noise'], units='K_RJ', start=2024-08-07 01:51:52.479 +00:00, duration=600.0s, sample_rate=50.0Hz, metadata={'atmosphere': True, 'sim_time': <Arrow [2025-08-24T19:00:10.323825+00:00]>, 'altitude': 5065.0, 'region': 'chajnantor', 'pwv': 0.925, 'base_temperature': 272.683})]

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={
        "window": {"name": "tukey", "kwargs": {"alpha": 0.1}},
        "remove_spline": {"knot_spacing": 30, "remove_el_gradient": True},
        "remove_modes": {"modes_to_remove": [0]},
    },
    map_postprocessing={
        "gaussian_filter": {"sigma": 1},
    },
    units="mK_RJ",
)

mapper.add_tods(tods)

output_map = mapper.run()

Mapping band f090: 100%|██████████| 2/2 [00:17<00:00,  8.51s/it, stokes=I, tod=2/2]
2025-08-24 19:01:12.011 INFO: Ran mapper for band f090 in 17.02 s.
Mapping band f150: 100%|██████████| 2/2 [00:15<00:00,  7.72s/it, stokes=I, tod=2/2]
2025-08-24 19:01:27.460 INFO: Ran mapper for band f150 in 15.44 s.

We can see the recovered map with

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