fullsky_geometry, made CarProjection abstract, broadcasting, tests

src/Pixell.jl

@@ -6,7 +6,9 @@ using FFTW
 using Printf
 
 include("enmap.jl")
+include("enmap_ops.jl")
 
-export Enmap
+export Enmap, CarClenshawCurtis
+export fullsky_geometry
 
 end

src/enmap.jl

@@ -1,10 +1,11 @@
 using Base: ViewIndex, @propagate_inbounds, AbstractCartesianIndex
 
-abstract type MapProjection end
-abstract type EquiCylProjection <: MapProjection end  # equidistant cylindrical projection
+abstract type AbstractMapProjection end
+abstract type EquiCylProjection <: AbstractMapProjection end  # equidistant cylindrical projection
+abstract type CarProjection <: EquiCylProjection end          # plate carrée
+
+struct CarClenshawCurtis <: CarProjection end      # plate carrée with pixels on poles and equator
 
-struct CarProjection <: EquiCylProjection end  # plate carrée
-struct CeaProjection <: EquiCylProjection end  # cylindrical equal area
 
 """
 Map type, contains an AbstractArray and a WCS object, but behaves like the
@@ -13,7 +14,7 @@ It only implements the subset of Base.Array operations which are common on maps.
 You should work with the data directly using `enmap_instance.data` if you need
 additional Array functions.
 """
-struct Enmap{T,N,AA<:AbstractArray,P<:MapProjection} <: AbstractArray{T,N}
+struct Enmap{T,N,AA<:AbstractArray,P<:AbstractMapProjection} <: AbstractArray{T,N}
     data::AA  # some kind of abstract array
     wcs::WCSTransform  # WCS object from WCS.jl
 end
@@ -22,9 +23,9 @@ end
 function Enmap(data::A, wcs, ::Type{PT}) where {A<:AbstractArray,PT}
     Enmap{eltype(A),ndims(A),A,PT}(data, wcs)
 end
-# create CAR maps by default
+# create CAR (Clenshaw-Curtis variant) maps by default
 function Enmap(data::A, wcs) where {A<:AbstractArray}
-    Enmap{eltype(A),ndims(A),A,CarProjection}(data, wcs)
+    Enmap{eltype(A),ndims(A),A,CarClenshawCurtis}(data, wcs)
 end
 
 Base.parent(x::Enmap) = x.data
@@ -87,7 +88,7 @@ struct NoWCS end
 combine(x::NoWCS, y) = y
 combine(x, y::NoWCS) = x
 combine(x::NoWCS, ::NoWCS) = x
-combine(x::Enmap, y::Enmap) = x  # TODO: check compatibility
+combine(x::WCSTransform, y::WCSTransform) = x  # TODO: check compatibility
 
 
 ####

src/enmap_ops.jl

@@ -0,0 +1,49 @@
+
+"""
+    fullsky_geometry([P=CarClenshawCurtis], res; shape = nothing, dims = ())
+
+# Arguments:
+- `proj=CarClenshawCurtis`: [optional] projection
+- `res`: resolution in radians. Passing a Number produces a square pixel.
+    Passing a tuple with (ΔRA, ΔDEC) produces a rectangular pixel.
+
+# Keywords
+- `shape::NTuple=nothing`: shape of the map. If not specified, will be computed.
+- `dims::NTuple=()`: additional dimensions to append to the shape, such as (3,) for IQU
+    to generate a map with `(nx, ny, 3)`.
+
+# Returns: 
+- `shape::Tuple, wcs::WCSTransform`: a tuple containing the shape of the map and the WCS
+
+# Examples
+```julia-repl
+julia> shape, wcs = fullsky_geometry(deg2rad(30/60))  # 30 arcmin pixel
+((720, 361), WCSTransform(naxis=2,cdelt=[-0.5, 0.5],crval=[0.25, 0.0],crpix=[360.5, 181.0]))
+```
+"""
+function fullsky_geometry(P::Type{<:CarClenshawCurtis}, res; shape = nothing, dims = ())
+    if isnothing(shape)
+        shape = (round.(Int, (2π, π) ./ res .+ (0, 1)))  # CAR has pixels on poles
+    end
+    nx, ny = shape
+    resx, resy = res
+
+    @assert abs(resx * nx - 2π) < 1e-8 "Horizontal resolution does not evenly divide the sky; this is required for SHTs."
+    @assert abs(resy * (ny - 1) - π) < 1e-8 "Vertical resolution does not evenly divide the sky; this is required for SHTs."
+
+    # Note the reference point is shifted by half a pixel to keep
+    # the grid in bounds, from ra=180+cdelt/2 to ra=-180+cdelt/2.
+    wcs = WCSTransform(2;
+        cdelt = [-360.0 / nx, 180.0 / (ny - 1)],
+        ctype = ["RA---CAR", "DEC--CAR"],
+        crpix = [floor(nx / 2) + 0.5, (ny + 1) / 2],
+        crval = [resy * 90 / π, 0])
+
+    return (nx, ny, dims...), wcs
+end
+
+fullsky_geometry(proj::Type{<:CarClenshawCurtis}, res::Number; shape = nothing, dims = ()) =
+    fullsky_geometry(proj, (res, res); shape = shape, dims = dims)
+
+fullsky_geometry(res; shape = nothing, dims = ()) =
+    fullsky_geometry(CarClenshawCurtis, res; shape = shape, dims = dims)

test/runtests.jl

@@ -1,6 +1,30 @@
 using Pixell
 using Test
 
-@testset "Pixell.jl" begin
-    # Write your tests here.
+
+@testset "Enmap geometry" begin
+    shape, wcs = fullsky_geometry(deg2rad(1 / 60))
+    @test wcs.cdelt ≈ [-0.016666666666666666, 0.016666666666666666]
+    @test wcs.crpix ≈ [10800.5, 5401.0]
+    @test wcs.crval ≈ [0.008333333333333333, 0.0]
+
+    shape, wcs = fullsky_geometry(deg2rad(1 / 61))
+    @test wcs.cdelt ≈ [-0.01639344262295082, 0.01639344262295082]
+    @test wcs.crpix ≈ [10980.5, 5491.0]
+    @test wcs.crval ≈ [0.00819672131147541, 0.0]
+
+    shape, wcs = fullsky_geometry(deg2rad(5); dims=(3,))
+    @test shape == (72, 37, 3)
+end
+
+
+@testset "Enmap broadcasting" begin
+    shape, wcs = fullsky_geometry(deg2rad(1); dims=(3,))
+    A, B = rand(shape...), rand(shape...)
+    ma = Enmap(A, wcs)
+    mb = Enmap(B, wcs)
+    @test A .+ B == ma .+ mb
+    @test A .+ B == ma .+ B
+    @test A .+ B == A .+ mb
+    @test A .+ B .* sin.(A.^2) == (ma .+ mb .* sin.(ma.^2))
 end