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√s  49.13 uK√s  0.5832’
   1  f150  150 GHz  30 GHz  0.5  12.25 aW√s  60 uK√s    104 uK√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
  location:
    lat: 23°01’45.84” S
    lon: 67°45’17.28” W
    alt: 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, 207MB/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, 177MB/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°, 0°, 0°)
  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,
                  el_bounds=(40, 90))

plan = planner.generate_plan(total_duration=1200,  # in seconds
                             scan_pattern="daisy",
                             scan_options={"radius": input_map.width.deg / 3},
                             sample_rate=50)

print(plan)
plan.plot()

Plan:
  duration: 1200 s
    start: 2024-08-06 03:00:00.000 +00:00
    end:   2024-08-06 03:20:00.000 +00:00
  location:
    lat: 23°01’45.84” S
    lon: 67°45’17.28” W
    alt: 5065 m
  sample_rate: 50 Hz
  center:
    ra:  19ʰ24ᵐ37.44ˢ
    dec: -31°13’48.00”
    az(mean): 198.9°
    el(mean): 81.24°
  scan_pattern: daisy
  scan_radius: 0.666°
  scan_kwargs: {'radius': np.float64(0.3333333333333333)}

import matplotlib.pyplot as plt

plt.scatter(0,0)
plt.scatter(2,2)

plt.annotate(xy=(0, 0), xycoords='data', xytext=(0.1, 0.1), textcoords="axes fraction", text="abc", bbox={"boxstyle": "round,pad=0.3", "fc": "w", "ec": "k", "alpha": 0.5, "lw": 1}, )

Text(0.1, 0.1, 'abc')

sim = maria.Simulation(
    instrument,
    plan=plan,
    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, 186MB/s]
Downloading https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/weather/era5/chajnantor.h5: 100%|██████████| 8.00M/8.00M [00:00<00:00, 158MB/s]
Constructing atmosphere: 100%|██████████| 10/10 [00:46<00:00,  4.64s/it]

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√s  49.13 uK√s  0.5832’
│    1  f150  150 GHz  30 GHz  0.5  12.25 aW√s  60 uK√s    104 uK√s      21”
├ Site:
│   region: chajnantor
│   location:
│     lat: 23°01’45.84” S
│     lon: 67°45’17.28” W
│     alt: 5065 m
│   seasonal: True
│   diurnal: True
├ Plan:
│   duration: 1200 s
│     start: 2024-08-06 03:00:00.000 +00:00
│     end:   2024-08-06 03:20:00.000 +00:00
│   location:
│     lat: 23°01’45.84” S
│     lon: 67°45’17.28” W
│     alt: 5065 m
│   sample_rate: 50 Hz
│   center:
│     ra:  19ʰ24ᵐ37.44ˢ
│     dec: -31°13’48.00”
│     az(mean): 198.9°
│     el(mean): 81.24°
│   scan_pattern: daisy
│   scan_radius: 0.666°
│   scan_kwargs: {'radius': np.float64(0.3333333333333333)}
├ Atmosphere(10 processes with 10 layers):
│ ├ spectrum:
│ │   region: chajnantor
│ └ weather:
│     region: chajnantor
│     altitude: 5065 m
│     time: Aug 5 23:09:59 -04:00
│     pwv[mean, rms]: (0.5 mm, 15 um)
└ 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°, 0°, 0°)
    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:

tod = sim.run()

print(tod)
tod.plot()

Generating turbulence: 100%|██████████| 10/10 [00:00<00:00, 10.25it/s]
Sampling turbulence: 100%|██████████| 10/10 [00:23<00:00,  2.34s/it]
Computing atmospheric emission: 100%|██████████| 2/2 [00:17<00:00,  8.68s/it, band=f150]
Sampling map: 100%|██████████| 2/2 [00:33<00:00, 16.81s/it, channel=(0 Hz, inf Hz)]
Generating noise: 100%|██████████| 2/2 [00:03<00:00,  1.78s/it, band=f150]

TOD(shape=(1798, 60000), fields=['atmosphere', 'map', 'noise'], units='K_RJ', start=2024-08-06 03:19:59.978 +00:00, duration=1200.0s, sample_rate=50.0Hz, metadata={'atmosphere': True, 'sim_time': <Arrow [2025-07-12T17:54:43.094920+00:00]>, 'altitude': 5065.0, 'region': 'chajnantor', 'pwv': 0.5, 'base_temperature': 272.519})

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(tod)

output_map = mapper.run()

Mapping band f090: 100%|██████████| 1/1 [00:00<00:00,  2.01it/s, band=f090, stokes=I]
2025-07-12 17:56:05.475 INFO: Ran mapper for band f090 in 13.8 s.
Mapping band f150: 100%|██████████| 1/1 [00:00<00:00,  2.01it/s, band=f150, stokes=I]
2025-07-12 17:56:18.748 INFO: Ran mapper for band f150 in 13.27 s.

We can see the recovered map with

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