add tests for Fejer1 and CC

simonsobs · Dec 16, 2023 · 89b9cda · 89b9cda
1 parent 6309e85
commit 89b9cda
Show file tree

Hide file tree

Showing 7 changed files with 48 additions and 9 deletions.
diff --git a/src/enmap.jl b/src/enmap.jl
@@ -217,7 +217,7 @@ function read_map(path::String; hdu::Int=1, sel=(), wcs::Union{WCSTransform,Noth
             @assert wcs0.ctype[1] == "RA---CAR" 
             @assert wcs0.ctype[2] == "DEC--CAR"
             if isfejer1(wcs0)
-                wcs = convert(Fejer1{Float64}, wcs0)
+                wcs = convert(CarFejer1{Float64}, wcs0)
             elseif isclenshawcurtis(wcs0)
                 wcs = convert(CarClenshawCurtis{Float64}, wcs0)
             else

diff --git a/src/projections/car_proj.jl b/src/projections/car_proj.jl
@@ -38,20 +38,29 @@ function Base.convert(::Type{CarFejer1{T}}, w0::WCSTransform) where T
 end
 
 function isfejer1(w0::WCSTransform)
-    cdelt, crpix, crval = getcdelt(w0), getcrpix(w0), getcrval(w0)
+    Δα, Δδ = getcdelt(w0)
+    skyloc = pix2sky(ones(Int, w0.naxis), w0, ones(w0.naxis))  # get first pixel
+    δ = rad2deg(skyloc[2])
     # reference DEC is an integer number of pixels from each pole minus half a pixel
-    half_from_pole = 90 - cdelt[2]/2
+    half_from_pole = 90 - Δδ/2
     npix1, npix2 = (half_from_pole - δ) / Δδ, (-half_from_pole - δ) / Δδ
-    return isapprox(δ + npix1 * Δδ, 90.0) && isapprox(δ + npix2 * Δδ, -90.0)
+    return isapprox(npix1, round(npix1)) && isapprox(npix2, round(npix2))
 end
 
-function isclenshawcurtis(w0::WCSTransform, tol=1e-6)
-    Δδ, j, δ = map(x->x[2], (getcdelt(w0), getcrpix(w0), getcrval(w0)))
+function isclenshawcurtis(w0::WCSTransform)
+    Δα, Δδ = getcdelt(w0)
+    skyloc = pix2sky(ones(Int, w0.naxis), w0, ones(w0.naxis))  # get first pixel
+    δ = rad2deg(skyloc[2])
     # reference DEC is an integer number of pixels from each pole
     npix1, npix2 = (90 - δ) / Δδ, (-90 - δ) / Δδ
-    return isapprox(δ + npix1 * Δδ, 90.0) && isapprox(δ + npix2 * Δδ, -90.0)
+    return isapprox(npix1, round(npix1)) && isapprox(npix2, round(npix2))
 end
 
+isclenshawcurtis(::CarClenshawCurtis) = true 
+isclenshawcurtis(x) = false 
+isfejer1(::CarFejer1) = true 
+isfejer1(x) = false 
+
 function Base.convert(::Type{WCSTransform}, w0::AbstractCARWCS)
     return WCSTransform(2;
         ctype = ["RA---CAR", "DEC--CAR"],

diff --git a/test/data/cc.fits b/test/data/cc.fits
diff --git a/test/data/fejer1.fits b/test/data/fejer1.fits
diff --git a/test/data_gen/gen_sht_test_data.py b/test/data_gen/gen_sht_test_data.py
@@ -8,9 +8,18 @@ def gen_simple(imap, powerlaw_exponent):
         for j in range(imap.shape[1]):
             imap[i,j] = (i * imap.shape[1] + j + 1)**powerlaw_exponent
 
-shape, wcs = enmap.fullsky_geometry(10.0*utils.degree)
+shape, wcs = enmap.fullsky_geometry(10.0*utils.degree, variant='CC')
 imap = enmap.zeros(shape, wcs=wcs)
 gen_simple(imap, 2)
+enmap.write_map("../data/cc.fits", imap)
+
+shape, wcs = enmap.fullsky_geometry(10.0*utils.degree, variant='fejer1')
+imap = enmap.zeros(shape, wcs=wcs)
+gen_simple(imap, 2)
+enmap.write_map("../data/fejer1.fits", imap)
+
+
+# %%
 
 alms = curvedsky.map2alm(imap, lmax=18, method="auto")
 np.savetxt("data/simple_analytic_sht.txt", np.column_stack([alms.real, alms.imag]), fmt='%.60g')

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -10,4 +10,4 @@ include("test_transforms.jl")
 include("test_distance_transform.jl")
 include("test_io.jl")
 include("test_plot.jl")
-include("test_utils.jl")
+include("test_utils.jl")
diff --git a/test/test_io.jl b/test/test_io.jl
@@ -14,3 +14,24 @@
     end
     # todo: add tests on IAU conversion w/ and w/o sel
 end
+
+##
+@testset "Enmap ClenshawCurtis vs Fejer1" begin
+
+    imap1 = read_map("data/cc.fits"; trim=false)
+    @test Pixell.isfejer1(imap1.wcs) == false
+    @test Pixell.isclenshawcurtis(imap1.wcs)
+
+    imap2= read_map("data/fejer1.fits"; trim=false)
+    @test Pixell.isfejer1(imap2.wcs)
+    @test Pixell.isclenshawcurtis(imap2.wcs) == false
+
+    imap1 = read_map("data/cc.fits"; trim=true)
+    @test Pixell.isfejer1(imap1.wcs) == false
+    @test Pixell.isclenshawcurtis(imap1.wcs)
+
+    imap2= read_map("data/fejer1.fits"; trim=true)
+    @test Pixell.isfejer1(imap2.wcs)
+    @test Pixell.isclenshawcurtis(imap2.wcs) == false
+
+end