Add dust map extension for SkyCoords

Project.toml

@@ -12,18 +12,22 @@ UnitfulAstro = "6112ee07-acf9-5e0f-b108-d242c714bf9f"
 
 [weakdeps]
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+SkyCoords = "fc659fc5-75a3-5475-a2ea-3da92c065361"
 
 [extensions]
 MakieExt = "Makie"
+SkyCoordsExt = "SkyCoords"
 
 [compat]
 BSplineKit = "0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19"
 DataDeps = "0.7"
 FITSIO = "0.13.0, 0.14, 0.15, 0.16.1, 0.17"
 Makie = "0.21.18, 0.22, 0.23, 0.24"
+SkyCoords = "1"
 Unitful = "0.17.0, 1"
 UnitfulAstro = "0.3.0, 0.4, 1"
 julia = "1.6"
 
 [extras]
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+SkyCoords = "fc659fc5-75a3-5475-a2ea-3da92c065361"

docs/Project.toml

@@ -4,13 +4,15 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DustExtinction = "fb44c06c-c62f-5397-83f5-69249e0a3c8e"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
+SkyCoords = "fc659fc5-75a3-5475-a2ea-3da92c065361"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulAstro = "6112ee07-acf9-5e0f-b108-d242c714bf9f"
 
 [compat]
 Documenter = "1"
 LaTeXStrings = "1"
 Measurements = "2"
+SkyCoords = "1"
 Unitful = "1"
 UnitfulAstro = "1"
 

docs/make.jl

@@ -1,15 +1,17 @@
 using Documenter
 using DustExtinction
 using CairoMakie
+using SkyCoords
 
 CairoMakie.activate!(type="png", px_per_unit=3)
 
 DocMeta.setdocmeta!(DustExtinction, :DocTestSetup, :(using DustExtinction); recursive = true)
 MakieExt = Base.get_extension(DustExtinction, :MakieExt)
+SkyCoordsExt = Base.get_extension(DustExtinction, :SkyCoordsExt)
 include("pages.jl")
 
 makedocs(;
-    modules = [DustExtinction, MakieExt],
+    modules = [DustExtinction, MakieExt, SkyCoordsExt],
     sitename = "DustExtinction.jl",
     format = Documenter.HTML(;
         prettyurls = get(ENV, "CI", nothing) == "true",

docs/src/dust_maps.md

@@ -43,11 +43,13 @@ dplot() # hide
 
 ## Advanced Usage
 
-Our dust maps also have native support for `Unitful.jl` and `Measurements.jl`.
+Our dust maps also have native support for `Unitful.jl`, `Measurements.jl`, and `SkyCoords.jl`.
 
 ```jldoctest
 julia> using Measurements, Unitful
 
+julia> using SkyCoords: GalCoords
+
 julia> using Unitful: °
 
 julia> l = 45°; b = 0°;
@@ -62,6 +64,8 @@ julia> l = l ± 0.1°; b = b ± 0.3°;
 julia> dustmap(l, b)
 6.4 ± 5.7 mag
 
+julia> dustmap(GalCoords(l, b))
+6.4 ± 5.7
 ```
 
 ## API/Reference

ext/SkyCoordsExt.jl

@@ -0,0 +1,11 @@
+module SkyCoordsExt
+    import DustExtinction: SFD98Map
+    using SkyCoords: AbstractSkyCoords, GalCoords, lon, lat
+
+    function (dustmap::SFD98Map)(s::AbstractSkyCoords)
+        s2 = convert(GalCoords, s)
+        l, b = lon(s2), lat(s2)
+        return dustmap(l, b)
+    end
+
+end # module

src/dust_maps.jl

@@ -68,22 +68,39 @@ end
 """
     (dustmap::SFD98Map)(l::Real, b::Real)
     (dustmap::SFD98Map)(l::Quantity, b::Quantity)
+    (dustmap::SFD98Map)(s::SkyCoords.AbstractSkyCoords)
 
 Get E(B-V) value from a `SFD98Map` instance at galactic coordinates (`l`, `b`),
 given in radians. Uses bilinear interpolation between pixel values. If `l` and
 `b` are `Unitful.Quantity` they will be converted to radians and the output
-will be given as `UnitfulAstro.mag`.
+will be given as `UnitfulAstro.mag`. If a `SkyCoords.AbstractSkyCoords` is passed, 
+it will be converted to galactic coordinates (requires Julia >= v1.9).
 
 # Example
 
-```jldoctest
+```jldoctest sfd98map
 julia> using DustExtinction
 
 julia> m = SFD98Map();
 
 julia> m(1, 2)
 0.013439524544325624
+```
+
+And now we can use a SkyCoords type as input:
+
+```jldoctest sfd98map
+julia> using SkyCoords
+
+julia> s = GalCoords(1, 2);
+
+julia> m(s) == m(1, 2)
+true
+```
+
+Use broadcasting to get E(B-V) values for multiple coordinates:
 
+```jldoctest sfd98map
 julia> l = 0:0.5:2; b = 0:0.5:2;
 
 julia> m.(l, b)
@@ -95,6 +112,19 @@ julia> m.(l, b)
   0.01862100327420125
 ```
 
+And now we can add angle units:
+
+```jldoctest sfd98map
+julia> using Unitful
+
+julia> m.(l * u"rad", b * u"rad")
+5-element Vector{Gain{Unitful.LogInfo{:Magnitude, 10, -2.5}, :?, Float64}}:
+ 99.69757461547852 mag
+  0.10180447359074371 mag
+  0.019595484241066132 mag
+  0.010238757633890877 mag
+  0.01862100327420125 mag
+```
 """
 function (dustmap::SFD98Map)(l::Real, b::Real)
     if b >= 0

test/Project.toml

@@ -3,11 +3,13 @@ DataDeps = "124859b0-ceae-595e-8997-d05f6a7a8dfe"
 Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SkyCoords = "fc659fc5-75a3-5475-a2ea-3da92c065361"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulAstro = "6112ee07-acf9-5e0f-b108-d242c714bf9f"
 
 [compat]
 DataDeps = "0.7"
+SkyCoords = "1"
 Unitful = "0.17.0, 1"
 UnitfulAstro = "0.3.0, 0.4, 1"

test/dust_maps.jl

@@ -45,6 +45,16 @@ using DataDeps
         @test Measurements.values.(output) ≈ refebv rtol = 0.02
     end
 
+    if VERSION ≥ v"1.9"
+        @testset "SkyCoords" begin
+            for i in eachindex(refebv)
+                @test dustmap(GalCoords(ref_l[i], ref_b[i])) ≈ refebv[i] rtol = 0.02
+                @test dustmap(convert(ICRSCoords, GalCoords(ref_l[i], ref_b[i]))) ≈ refebv[i] rtol = 0.02
+            end
+            @test dustmap.(convert.(Ref(ICRSCoords), GalCoords.(ref_l, ref_b))) ≈ refebv rtol = 0.02
+        end
+    end
+
     @testset "Unitful" begin
         ref_l_u = ref_l * u"rad"
         ref_b_u = ref_b * u"rad"

test/runtests.jl

@@ -1,6 +1,6 @@
 using DustExtinction
 using DustExtinction: bounds, checkbounds
-using Test, Measurements, Unitful, UnitfulAstro, Random
+using Test, Measurements, SkyCoords, Unitful, UnitfulAstro, Random
 
 @testset "DustExtinction" begin
     Random.seed!(9994445781)