Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 3 additions & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
Package: magicaxis
Type: Package
Title: Pretty Scientific Plotting with Minor-Tick and Log Minor-Tick Support
Version: 2.6.2
Date: 2026-07-28
Version: 2.7.0
Date: 2026-07-29
Authors@R: person(given = "Aaron",
family = "Robotham",
role = c("aut", "cre"),
Expand All @@ -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
2 changes: 1 addition & 1 deletion NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -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")
72 changes: 72 additions & 0 deletions R/magMWdust.R
Original file line number Diff line number Diff line change
@@ -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
#
#=========================================

utils::globalVariables("SFD_dust")

# 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))
}
23 changes: 15 additions & 8 deletions R/magproj.R
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand All @@ -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){
Expand Down Expand Up @@ -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{
Expand Down
Binary file added data/SFD_dust.rda
Binary file not shown.
99 changes: 99 additions & 0 deletions man/magMWdust.Rd
Original file line number Diff line number Diff line change
@@ -0,0 +1,99 @@
\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

#the below will take about a minute to generate
\dontrun{
magproj(0,0,type='n')
magMWdust(type = "pl") #add MW dust map as polygons
}
}
\keyword{plot}
5 changes: 4 additions & 1 deletion man/magproj.Rd
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ magproj(long, lat, type = "b", plottext, longlim = c(-180, 180), latlim = c(-90,
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, ...)
labloc = c(90, -45), labeltype = "deg", crunch = FALSE, lab.cex = 1, ...)

magprojgrid(nlat = 6, nlong = 6, prettybase = 30, box = TRUE, ...)

Expand Down Expand Up @@ -93,6 +93,9 @@ Should the labels be drawn using degrees (deg) or colon delimited sexigesimal (s
}
\item{crunch}{
If set to FALSE the full output of \code{\link[celestial]{deg2hms}} and \code{\link[celestial]{deg2dms}} is printed. If set to TRUE a simplified output is used, where only the hours and degrees parts are extracted and appended with a 'h' and a degree symbol respectively.
}
\item{lab.cex}{
Character expansion factor for the coordinate labels. Default is 1. Passed as \option{cex} to the internal \code{\link{text}} calls that draw the longitude and latitude labels.
}
\item{\dots}{
For magproj, Extra options that are either passed to \code{\link{points}} (\option{type='p'}), \code{\link{lines}} (\option{type='l'}), \code{\link{polygon}} (\option{type='pl'}), \code{\link{text}} (\option{type='t'}), or \code{\link{polygon}} (\option{type='b'}). For magprojgrid dots are pased to \code{\link{lines}} for drawing the grid lines. For magprojlabels dots are passed to \code{\link{text}} for adding text labels.
Expand Down
Loading