add magMWdust function

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: magicaxis
 Type: Package
 Title: Pretty Scientific Plotting with Minor-Tick and Log Minor-Tick Support
-Version: 2.6.0
-Date: 2026-05-14
+Version: 2.7.0
+Date: 2026-06-04
 Authors@R: person(given = "Aaron",
                       family = "Robotham",
                       role = c("aut", "cre"),
@@ -14,3 +14,4 @@ Imports: grDevices, graphics, stats, celestial (>= 1.4.1), MASS, plotrix, sm, ma
 Remotes: asgr/ParmOff
 VignetteBuilder: knitr
 Config/testthat/edition: 3
+LazyData: true

NAMESPACE

@@ -14,5 +14,5 @@ importFrom("RANN", "nn2")
 importFrom("plotrix", "color.legend", "draw.ellipse")
 importFrom("sm", "sm.density")
 importFrom("MASS", "kde2d")
-importFrom("utils", "str")
+importFrom("utils", "str", "setTxtProgressBar", "txtProgressBar")
 importFrom("ParmOff", "ParmOff")

R/magMWdust.R

@@ -0,0 +1,72 @@
+#=========================================
+#
+# File Name : magMWdust.R
+# Created By : awright
+# Creation Date : 03-06-2026
+# Last Modified : Wed Jun  3 13:40:13 2026
+#
+#=========================================
+
+
+# Define a plotting helper that overlays Milky Way dust on a projected sky plot.
+magMWdust <- function(dust.data = NULL, dlon = NULL, dlat = NULL, type = "p", pch = 16, pt.cex = 0.5, opacity.range = c(0, 0.5), whiteblack.percentile = c(0.5, 0.95), stretch = "lin", min.opacity.plot = 0.01, show.status = TRUE, ...) {
+
+  # Restrict the drawing style to points or polygons.
+  if (!type %in% c("pl", "p")) stop("magMWdust function expects type of 'p' (for points) or 'pl' (for polygons) only")
+  # If none provided, read the dust map data, which is a dlon=dlat=1 sampling  
+  if (is.null(dust.data)) { 
+    # Define the dlon and dlat values
+    dlon <- dlat <- 1
+    # Lazy load the SFD_dust data
+    dust_all<-SFD_dust
+  } else { 
+    if (!is.data.frame(dust.data)) {
+      stop("dust.data is not a data frame; load an example with data(SFD_dust)") 
+    }
+    if (!all(c("ebv","ra","dec")%in%colnames(dust.data))) {
+      stop("dust.data is missing required components; load an example with data(SFD_dust)") 
+    }
+    dust_all<-dust.data 
+    if (type=='pl') { 
+      if (is.null(dlon)) stop("dlon must be provided when providing input dust.data and using type == 'pl'") 
+      if (is.null(dlat)) stop("dlat must be provided when providing input dust.data and using type == 'pl'") 
+      if (!is.numeric(dlon)) stop("dlon must be numeric") 
+      if (!is.numeric(dlat)) stop("dlat must be numeric") 
+    }
+  }
+  # Map dust values onto an opacity scale for plotting.
+  dust_all$map <- magicaxis::magmap(dust_all$ebv, range = opacity.range, hicut = whiteblack.percentile[2], locut = whiteblack.percentile[1], stretch = stretch)$map
+  # Drop grid cells that would be too faint to plot usefully.
+  dust <- dust_all[which(dust_all$map > min.opacity.plot), ]
+  # Stop on empty result 
+  if (nrow(dust)==0) stop("threshold for min.opacity.plot causes no data to be plotted. Reduce value to produce a result")
+
+  # Draw filled polygons when polygon mode has been requested.
+  if (type == "pl") {
+    # Open a progress bar for the per-cell polygon loop.
+    if (interactive() & isTRUE(show.status)) { 
+      pb <- txtProgressBar(style = 3, min = 1, max = nrow(dust))
+    }
+    # Iterate over each retained sky cell.
+    for (i in 1:nrow(dust)) {
+      # Project and draw the four corners of the current sky cell.
+      magicaxis::magproj(c(dust$ra[i]  - dlon/2, dust$ra[i]  - dlon/2, dust$ra[i]  + dlon/2, dust$ra[i]  + dlon/2), 
+                         c(dust$dec[i] - dlat/2, dust$dec[i] + dlat/2, dust$dec[i] + dlat/2, dust$dec[i] - dlat/2), 
+                         type = type, add = TRUE, col = hsv(v = 0, alpha = dust$map[i]), border=NA, ...)
+      # Advance the progress bar after drawing the current polygon.
+      if (interactive() & isTRUE(show.status)) {
+        setTxtProgressBar(pb, i)
+      } 
+    }
+    # Close the progress bar when the polygon layer is complete.
+    if (interactive() & isTRUE(show.status)) {
+      close(pb)
+    }
+  } else {
+    # Draw the retained dust grid cells as projected points.
+    magicaxis::magproj(dust$ra, dust$dec, type = "p", add = TRUE, pch = pch, cex = pt.cex, col = hsv(v = 0, alpha = dust$map), ...)
+  }
+  # Return invisibly because this function is used for its plotting side effects.
+  return(invisible(NULL))
+}
+

R/magproj.R

@@ -15,7 +15,7 @@ magproj=function(long, lat, type='b', plottext, longlim=c(-180,180), latlim=c(-9
                  projection="aitoff", parameters=NULL, centre=c(0,0), add=FALSE,
                  fliplong=FALSE, nlat=6, nlong=6, prettybase=30, labels=TRUE, grid=TRUE,
                  grid.col='grey', grid.lty=2, auto=FALSE, upres=100, box=TRUE, labloc=c(90,-45),
-                 labeltype='deg', crunch=FALSE, ...){
+                 labeltype='deg', crunch=FALSE, lab.cex=1, ...){
 
   if(is.matrix(long) | is.data.frame(long)){
     lat = long[, 2]
@@ -42,7 +42,14 @@ magproj=function(long, lat, type='b', plottext, longlim=c(-180,180), latlim=c(-9
   }
 
   long= (long+(180-orientation[2])) %% 360 - (180-orientation[2])
-  lat= (lat+90) %% 180 - 90
+  .lat_wrap <- function(x) { 
+    x[abs(x)==180]<-0
+    x=(x+180)%%360-180
+    x[which(x>= +90)]<- +90-x[which(x>= +90)]%%90
+    x[which(x<= -90)]<- -90+abs(x[which(x<= -90)])%%90
+    return=x
+  }
+  lat= .lat_wrap(lat)
 
   if(add==FALSE){
     if(auto==TRUE){
@@ -97,16 +104,16 @@ magproj=function(long, lat, type='b', plottext, longlim=c(-180,180), latlim=c(-9
       latpretty=latgrid$tickat
       latpretty=latpretty[latpretty>latlim[1] & latpretty<latlim[2]]
       temp=mapproject(longpretty, rep(labloc[2],length(longpretty)))
-      if(labeltype=='deg'){text(temp,labels = longpretty %% 360)}
+      if(labeltype=='deg'){text(temp,labels = longpretty %% 360,cex = lab.cex)}
       if(labeltype=='sex'){
-        if(crunch==FALSE){text(temp,labels = deg2hms(longpretty %% 360,type='cat'))}
-        if(crunch==TRUE){text(temp,labels = paste(deg2hms(longpretty %% 360,type='mat')[,1],'h',sep=''))}
+        if(crunch==FALSE){text(temp,labels = deg2hms(longpretty %% 360,type='cat'),cex = lab.cex)}
+        if(crunch==TRUE){text(temp,labels = paste(deg2hms(longpretty %% 360,type='mat')[,1],'h',sep=''),cex = lab.cex)}
       }
       temp=mapproject(rep(labloc[1],length(latpretty)), latpretty)
-      if(labeltype=='deg'){text(temp,labels = latpretty)}
+      if(labeltype=='deg'){text(temp,labels = latpretty,cex = lab.cex)}
       if(labeltype=='sex'){
-        if(crunch==FALSE){text(temp,labels = deg2dms(latpretty,type='cat'))}
-        if(crunch==TRUE){text(temp,labels = paste(deg2dms(latpretty,type='mat')[,1],'\u00B0',sep=''))}
+        if(crunch==FALSE){text(temp,labels = deg2dms(latpretty,type='cat'),cex = lab.cex)}
+        if(crunch==TRUE){text(temp,labels = paste(deg2dms(latpretty,type='mat')[,1],'\u00B0',sep=''),cex = lab.cex)}
       }
     }
   }else{

man/magMWdust.Rd

@@ -0,0 +1,100 @@
+\name{magMWdust}
+\alias{magMWdust}
+\alias{SFD_dust}
+\title{
+Plot a Milky Way dust overlay
+}
+\description{
+Adds a projected Milky Way dust layer to an existing sky plot using a lookup-table 
+constructed from the SFD dust map on an equatorial grid, and drawing the result as 
+points or filled polygons.
+}
+\usage{
+magMWdust(dust.data = NULL, 
+  dlon = NULL, dlat = NULL, 
+  type = "p", pch = 16, pt.cex = 0.5,
+  opacity.range = c(0, 0.5),
+  whiteblack.percentile = c(0.5, 0.95),
+  stretch = "lin",
+  min.opacity.plot = 0.01, show.status = TRUE, ...)
+}
+\arguments{
+  \item{dust.data}{
+The data set to use for drawing the dust map. If \code{NULL} (the default), the function 
+lazy-loads the SFD dust map data provided with the package (\code{data("SFD_dust")}). The 
+code expects this data to be a data frame containing columns "ra", "dec", and "ebv". 
+}
+  \item{dlon}{
+Step size between element of longitude. Required when specifying dust.data directly, 
+and uses \code{dlon=1} when lazy-loading the internal dust map. 
+}
+  \item{dlat}{
+Step size between element of latitude. Required when specifying dust.data directly, 
+and uses \code{dlat=1} when lazy-loading the internal dust map. 
+}
+  \item{type}{
+Character scalar giving the drawing mode: \code{"p"} for points or \code{"pl"}
+for polygons.
+}
+  \item{pch}{
+Plotting symbol used when \code{type = "p"}.
+}
+  \item{pt.cex}{
+Point expansion factor used when \code{type = "p"}.
+}
+  \item{opacity.range}{
+Two-element numeric vector giving the output opacity range passed to
+\code{\link{magmap}}.
+}
+  \item{whiteblack.percentile}{
+Two-element numeric vector giving the lower and upper percentiles used when
+mapping dust values to opacity.
+}
+  \item{stretch}{
+Stretch mode passed through to \code{\link{magmap}}.
+}
+  \item{min.opacity.plot}{
+Minimum mapped opacity required for a grid cell to be drawn.
+}
+  \item{show.status}{
+Option to show a progress bar when drawing polygons. 
+}
+  \item{\dots}{
+Additional arguments passed to \code{\link{magproj}}.
+}
+}
+\details{
+The function plots \code{E(B-V)} values from a provided dust map data.frame, or if none is provided uses a lazy-loaded
+SFD dust map which is sampled on a 1x1 deg grid. This lookup table is then mapped to points (or polygons, see below) 
+with opacity determined by \code{\link{magmap}}, using the requested white/black points defined using percentiles 
+(\code{whiteblack.percentile}) of the map, and throwing away cells whose opacity is below the requested 
+minimum (\code{min.opacity.plot}).
+
+When \code{type = "pl"} each retained grid cell is drawn as a projected
+four-corner polygon. Otherwise the retained cells are drawn as projected
+points. Polygon drawing is very slow when drawing to X11 (many minutes) except 
+with very low resolutions. Polygon drawing may nonetheless be preferable when 
+writing to file, where the operation is much faster and reduces aliasing effects 
+caused by over- (or under-)lapping points. 
+}
+\value{
+Returns \code{NULL} invisibly. The function is called for its plotting side
+effects.
+}
+\references{
+  Schlegel, Finkbeiner, & Davis, 1998, ApJ, 500, 2, 525
+}
+\author{
+  Angus H Wright
+}
+\seealso{
+\code{\link{magmap}}, \code{\link{magproj}}
+}
+\examples{
+magproj(0,0,type='n') #Construct a base magproj plot 
+magMWdust() #add MW dust map as points 
+\dontrun{
+magMWdust(type = "pl") #add MW dust map as polygons 
+}
+}
+\keyword{plot}