Simulations natively in CAR

[zonca]

Should we reproject or create simulations in CAR natively?

[zonca]

One option would be: improved based on feedback below

Run PySM to alm and reproject

• Run PySM in HEALPix (optionally use higher resolution than normally required to help with reprojection error, 8192 for LAT, 1024 for SAT), galactic components are only available until NSIDE 4096, but would ud_grade to 8192 reduce reprojection error? probably not if we reproject in harmonics space
• Do bandpass integration in map domain
• Do map2alm, apply smoothing and rotation to equatorial
• From full sky spherical alm -> CAR
• Noise maps instead are created directly in CAR

[aiolasimone]

From ACT experience, that plan sounds good.

You can definitely go from HEALPix to CAR via alm, see this pixell function.

If your PySM sims are in galactic coordinates, you do need to rotate the alms from galactic to equatorial. The function that I have linked above supports gal->equ rotation.

On the long run having PySM run on equatorial (and CAR) directly may save a lot of time (though I am not super familiar with PySM).

Noise definitely in equatorial from the beginning, and sure I would project to CAR right away.

[ajvanengelen]

The lensing convergence map in the sims was converted to CAR in this way - the original map was in HEALPIX, then I projected to CAR using harmonic transforms. This allowed me to use the lensing routine in pixell.

Not sure how some of our small-scale science (e.g. measuring tSZ cluster profiles) is affected by this.

[jcolinhill]

I think Andrea's suggestion of use higher resolution than normally required would ideally be enough to avoid serious problems with small-scale science, but we might need higher resolution than WebSky currently provides (can WebSky be run at arbitrary Nside?)

[zonca]

we need a volunteer that explores different options, quantifies the errors and decides the best strategy.

[ajvanengelen]

For the time being, what we have in hand are the websky maps at nside = 4096, plus the Galactic maps also in Healpix,, so the baseline plan should probably be to reproject those. Note that some care will have to be applied with the two different pixel window functions.

[zonca]

yes, does someone volunteers to do the reprojection of the 3 releases?:

• https://github.com/simonsobs/map_based_simulations/tree/master/201904_highres_foregrounds_equatorial
• https://github.com/simonsobs/map_based_simulations/tree/master/201904_highres_foregrounds_variable_spectral_index
• https://github.com/simonsobs/map_based_simulations/tree/master/201905_extragalactic

please contribute the script to a "scripts" folder into mapsims, if you check the release folder, it has subfolders 4096 and 512, you can add CAR_? or something simiilar at that level. then modify the release README to describe the CAR maps.

For noise, can someone volunteer to make a pull request to https://github.com/simonsobs/mapsims/blob/master/mapsims/noise.py to generate noise in CAR? I'll then run the sims. You need to install pysm 3 from github.com/healpy/pysm (clone master and pip install . it)

[msyriac]

This plan (alm reproject foregrounds, directly project lensed sims, directly generate noise sims) generally sounds good to me. A few things:

1. for ACT work, I use this script for reprojection of Planck : https://github.com/ACTCollaboration/tile-c/blob/master/bin/planck/planck2car.py making sure to use the "map" option for diffuse maps and the "ivar" option for hit counts. We should get this adapted and cleaned up into mapsims if possible.
2. The hit count maps for noise will have to be reprojected with interpolation (not in alm space) using the "ivar" option above. We will have to think about what resolution to store these at, and what kind of flexibility we want to allow in possible CAR resolutions for the sims. (equivalent to choosing what nside we produce maps at)
3. an aspiration of mine is that the Websky sims will some day be natively generated in CAR pixels at 0.5 arcmin resolution. Without this, we will have ringing in the clusters and sources, which could affect work done by the SZ group. (Pinging @marcelo-alvarez and @mattyowl )

I cannot volunteer myself for any of this, this week. But might be able to next week. If someone else wants to pick it up, it would be appreciated.

[zonca]

@zonca to split this in 2

[marcelo-alvarez]

@msyriac Would it be useful to start with websky-tSZ and CIB in CAR (i.e. leave out kappa and kSZ) for now?

[msyriac]

We could just do something very simple where we loop through each sim in the final output directory and alm reproject and save it with _car.fits. This is fine since the signal is provided separately from the noise anyway. Is there any reason to prefer doing this at the "template" level? I haven't checked, but I would guess that the number of reprojections needed is not much larger.

[zonca]

as it is very simple, let's do this first, doing it at template level requires PySM 3 to be able to handle CAR pixelization, which is doable but it would require additional development time.

[zonca]

@aiolasimone, as discussed yesterday, it would be great if you can create a Jupyter Notebook that gets a full sky map in healpix, for example (at NERSC):

• /project/projectdirs/sobs/v4_sims/mbs/201906_highres_foregrounds_extragalactic_tophat/512/combined/0000/simonsobs_combined_uKCMB_laMFF2_nside512_0000.fits
• /project/projectdirs/sobs/v4_sims/mbs/201906_highres_foregrounds_extragalactic_tophat/4096/combined/0000/simonsobs_combined_uKCMB_laMFF2_nside4096_0000.fits

reprojects to CAR and create plots to compare the 2 (that is good for my reference so I can check results).

Doing this you should define what are good resolution parameters for CAR both at 4096 and 512.

[aiolasimone]

This should do it. We may want to have some of the AWGs to verify that the re-projection looks fine. https://gist.github.com/aiolasimone/5f0016dbc0fa160477eea7a45b6eec5b

[msyriac]

This looks good to me. With @mhasself and @amaurea , we may want to standardize exactly what CAR shape,wcs to use, based on what the map-maker will produce. This is closely related to the discussion in #18 whose results we should enshrine in the analysis interface document.

[zonca]

@mhasself @amaurea any feedback about CAR parameters?

[mhasself]

I can help with constructing valid WCS but some map user should confirm that this is the sort of thing you want, have used in the past, know to work, etc. The following is implemented in pixell.

For RA, we fix the system by insisting to have a column of pixels centered at RA=0. Then one gets to choose where to place the branch cut; call this RA_MAX (corresponding to the centers of the pixels of the left edge of the map). For resolution RES this leads to the form:

  naxis1 = round(360/RES)
  crpix1 = RA_MAX/RES + 0.5
  crval1 = RES/2
  cdelt1 = -RES

The half-pixel shift in crpix/crval is needed so that the column of pixels centered at RA_MAX are entirely inside the map.

For declination, centering a row of pixels at dec=0 leads naturally to also having whole rows of pixels centered on the poles. That's odd, but seems to be allowed, and is good for spherical harmonic transforms (?). That leads to:

  naxis2 = round(180/RES) + 1
  crpix2 = 90/RES + 1
  crval2 = 0.
  cdelt2 = RES

Often one will be dealing with smaller patches, which we will probably want to be pixel-compatible with the full sky version. I believe pixel-compatibility is guaranteed provided that, for a given resolution, crpix differs by an integer and crval differs by an integer multiple of RES.

[msyriac]

Ok, this doesn't quite match what we have in an AdvACT maps at 0.5 arcmin resolutions:

NAXIS1  =                43200                                                  
NAXIS2  =                10320                                                  
NAXIS3  =                    3                                                  
WCSAXES =                    2 / Number of coordinate axes                      
CRPIX1  =              21601.0 / Pixel coordinate of reference point            
CRPIX2  =               7561.0 / Pixel coordinate of reference point            
CDELT1  =  -0.0083333333333333 / [deg] Coordinate increment at reference point  
CDELT2  =   0.0083333333333333 / [deg] Coordinate increment at reference point  
CUNIT1  = 'deg'                / Units of coordinate increment and value        
CUNIT2  = 'deg'                / Units of coordinate increment and value        
CTYPE1  = 'RA---CAR'           / Right ascension, plate caree projection        
CTYPE2  = 'DEC--CAR'           / Declination, plate caree projection            
CRVAL1  =                  0.0 / [deg] Coordinate value at reference point      
CRVAL2  =                  0.0 

There is value in having the highest resolution maps be perfectly WCS compatible with the AdvACT maps, and as Matthew suggested elsewhere and in the discussion in https://github.com/simonsobs/map_based_simulations/issues/18 have the lower resolution maps just scale the number of pixels by integer values.

[marcelo-alvarez]

I think Andrea's suggestion of use higher resolution than normally required would ideally be enough to avoid serious problems with small-scale science, but we might need higher resolution than WebSky currently provides (can WebSky be run at arbitrary Nside?)

We have plans to re-run everything in websky at 8192 on a short time scale, hopefully that will be sufficient for the time being.

[mhasself]

@msyriac The AdvACT WCS is compatible with the kind I proposed... The generalization of my earlier statement is "pixel-compatibility is guaranteed for two CAR maps if the crpix differs by (n+f) and crval differs by (m+f)*RES, where n and m are integers."

However, WCS validity on the full RA range, using only 360/RES pixels, requires putting the reference point on a pixel edge, which means crpix1 should be something-point-five.

[zonca]

split implementation of CAR for noise to #28 for foregrounds I think I have all I need to create a first version of the rotation script, I'll ask for more feedback then.

[zonca]

PySM can now produce CAR maps