Degrading existing fracdet maps
Closed this issue · 5 comments
Say I have a nside=4096 fracdet map and I want to degrade it to an nside=32 fracdet map. I don’t want to use any of the reduction specified by reduce_array (i.e., the mean, median, etc.), rather, I want to sum all the valid pixels (i.e., pixels not marked as hp.UNSEEN) and divide by the total number of sparse pixels. I think this is an easy operation to do (I just need to add sum to reduce_array, np.nansum exists, so this is trivial), but I thought there might already be another way to do this.
I would have expected one of the following to work.
import numpy as np
import pylab as plt
import healpy as hp
import healsparse as hsp
nside = 64
nside2 = 16
npix = hp.nside2npix(nside)
pix = np.arange(npix)
fracdet = pix/float(npix)
fracdet[pix >= npix/4] = 0.0
fracdet[pix >= npix/2] = hp.UNSEEN
# Original fracdet map
hp.mollview(fracdet,nest=True,title='fracdet (nside=%s)'%nside)
# Downsampled map
hsp_map = hsp.HealSparseMap(nside_coverage=nside2, healpix_map=fracdet)
hp.mollview(hsp_map.fracdet_map(nside2).generate_healpix_map(nside2),
title='fracdet (nside=%s; zero & UNSEEN)'%nside2,nest=True)
# Change the fracdet
fracdet[fracdet == hp.UNSEEN] = 0.0
hsp_map = hsp.HealSparseMap(nside_coverage=nside2, healpix_map=fracdet)
hp.mollview(hsp_map.fracdet_map(nside2).generate_healpix_map(nside2),
title='fracdet (nside=%s; zero only)'%nside2,nest=True)
# Zero sentinel
hsp_map = hsp.HealSparseMap(nside_coverage=nside2, healpix_map=fracdet,sentinel=0.0)
hp.mollview(hsp_map.fracdet_map(nside2).generate_healpix_map(nside2),
title='fracdet (nside=%s; zero only; zero sentinel)'%nside2,nest=True)
This does seem work
# Zero sentinel; degrade
hsp_map = hsp.HealSparseMap(nside_coverage=nside2, healpix_map=fracdet,sentinel=0.0)
hp.mollview(hsp_map.degrade(nside2).generate_healpix_map(nside2),
title='fracdet (nside=%s; reduce)'%nside2,nest=True)
I think I can explain what's going on. And I'm not sure how to document this clearly, which is another issue.
The main issue is that a fracdet map should have 0 as filler and not UNSEEN. Because the fractional coverage of an uncovered pixel is actually 0 and not undefined. The next confusing thing is that in order to change the resolution of fracdet map you want to degrade it (as you have shown) and not take the fracdet of a fractdet map. But the fracdet of a fracdet map is actually 100% as you've shown, because every pixel is defined in a fracdet map, it just happens to be 0!
One of your maps, the "zero sentinel" map, actually doesn't make sense because hsp_map = hsp.HealSparseMap(nside_coverage=nside2, healpix_map=fracdet,sentinel=0.0) will not replace the sentinel in the input map with 0, it just says that uncovered pixels should be 0 -- and since you input a map with UNSEEN pixels then those all count as valid values (that do not equal to 0). The sentinel here does not say "change the values", it says "this is the sentinel value in the map that I'm inputting".
So the last one that works it may be partially a coincidence. If you replace the UNSEEN with 0, and degrade, you will get the values that you want.
Does that help?
Thinks Eli, I agree that this is primarily a documentation problem. Maybe have fracdet_map warn if the input map contains values other than 1 or sentinel? Maybe more explicit documentation or an example?
One clarifying point, prior to the creation of the third map I replaced all UNSEEN values with zeros, i.e.,
fracdet[fracdet == hp.UNSEEN] = 0.0
So the fracdet map for the 3rd, 4th, and 5th tests no longer contains UNSEEN values. When I set zero as the sentinel, the behavior makes sense, as fracdet_map is counting every non-zero pixel as covered (i.e., fracdet value of 1). The degraded map functions in the 5th test functions as expected.