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)
[2]:
array = {"field_of_view": 0.5,
         "shape": "circle",
         "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  2116  29.85’      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”
../_images/tutorials_custom-map-simulations_3_1.png

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

[3]:
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 *= 1e1

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, 184MB/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: -3.845e-05
    max: -5.005e-08
  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, psi): (0°, 0°, 0°)
  memory: 16.78 MB
../_images/tutorials_custom-map-simulations_5_2.png
[4]:
site = maria.get_site("llano_de_chajnantor", altitude=5065)

print(site)
site.plot()
Site:
  region: chajnantor
  location: 23°01’45.84”S 67°45’17.28”W
  altitude: 5.065 km
  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, 196MB/s]
../_images/tutorials_custom-map-simulations_6_2.png
[5]:
plan = maria.Plan(
    start_time="2024-08-06T03:00:00",
    scan_pattern="daisy",
    scan_options={"radius": 0.5, "speed": 0.1},  # in degrees
    duration=900,  # in seconds
    sample_rate=50,  # in Hz
    scan_center=(291.156, -31.23),
    frame="ra_dec")

print(plan)
plan.plot()
Plan:
  start_time: 2024-08-06 03:00:00.000 +00:00
  duration: 900 s
  sample_rate: 50 Hz
  center:
    ra: 19ʰ24ᵐ37.44ˢ
    dec: -31°13’48.00”
  scan_pattern: daisy
  scan_radius: 0.9989°
  scan_kwargs: {'radius': 0.5, 'speed': 0.1}
../_images/tutorials_custom-map-simulations_7_1.png
[6]:
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, 199MB/s]
Downloading https://github.com/thomaswmorris/maria-data/raw/master/atmosphere/weather/era5/chajnantor.h5: 100%|██████████| 8.00M/8.00M [00:00<00:00, 157MB/s]
Constructing atmosphere: 100%|██████████| 10/10 [00:57<00:00,  5.79s/it]
Simulation
├ Instrument(1 array)
│ ├ arrays:
│ │             n     FOV baseline        bands polarized
│ │  array1  2116  29.85’      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: 23°01’45.84”S 67°45’17.28”W
│   altitude: 5.065 km
│   seasonal: True
│   diurnal: True
├ Plan:
│   start_time: 2024-08-06 03:00:00.000 +00:00
│   duration: 900 s
│   sample_rate: 50 Hz
│   center:
│     ra: 19ʰ24ᵐ37.44ˢ
│     dec: -31°13’48.00”
│   scan_pattern: daisy
│   scan_radius: 0.9989°
│   scan_kwargs: {'radius': 0.5, 'speed': 0.1}
├ Atmosphere(10 processes with 10 layers):
│ ├ spectrum:
│ │   region: chajnantor
│ └ weather:
│     region: chajnantor
│     altitude: 5.065 km
│     time: Aug 5 23:07:29 -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: -1.915e-04
      max: -2.492e-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, psi): (0°, 0°, 0°)
    memory: 16.78 MB
[7]:
tod = sim.run()

print(tod)
tod.plot()
Generating turbulence: 100%|██████████| 10/10 [00:00<00:00, 12.49it/s]
Sampling turbulence: 100%|██████████| 10/10 [00:21<00:00,  2.13s/it]
Computing atmospheric emission: 100%|██████████| 2/2 [00:15<00:00,  7.73s/it, band=f150]
Sampling map: 100%|██████████| 2/2 [00:28<00:00, 14.45s/it, channel=[ 0. inf] Hz]
Generating noise: 100%|██████████| 2/2 [00:03<00:00,  1.53s/it, band=f150]
TOD(shape=(2116, 45000), fields=['atmosphere', 'map', 'noise'], units='pW', start=2024-08-06 03:14:59.979 +00:00, duration=900.0s, sample_rate=50.0Hz, metadata={'atmosphere': True, 'sim_time': <Arrow [2025-05-30T18:50:54.006404+00:00]>, 'altitude': 5065.0, 'region': 'chajnantor', 'pwv': 0.5, 'base_temperature': 272.523})
../_images/tutorials_custom-map-simulations_9_2.png
[8]:
from maria.mappers import BinMapper

mapper = BinMapper(
    center=(291.156, -31.23),
    frame="ra_dec",
    width=1.,
    height=1.,
    resolution=1. / 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.41it/s, band=f090, stokes=I]
2025-05-30 18:51:20.143 INFO: Ran mapper for band f090 in 20.49 s.
Mapping band f150: 100%|██████████| 1/1 [00:00<00:00,  2.41it/s, band=f150, stokes=I]
2025-05-30 18:51:39.763 INFO: Ran mapper for band f150 in 19.62 s.

We can see the recovered map with

[9]:
output_map.plot(nu_index=[0, 1])
../_images/tutorials_custom-map-simulations_12_0.png