Skip to content

Commit

Permalink
Merge pull request #58 from kateharborne/dev-obs_summary
Browse files Browse the repository at this point in the history 
dev-obs_summary
  • Loading branch information
@kateharborne
kateharborne committed Jul 28, 2022
2 parents be0c281 + a33d9f0 commit c7114e4
Show file tree
Hide file tree
Showing 11 changed files with 643 additions and 81 deletions.
2 changes: 1 addition & 1 deletion DESCRIPTION
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
Package: SimSpin
Type: Package
Title: SimSpin - A package for the kinematic analysis of galaxy simulations
Version: 2.3.10
Version: 2.3.11
Author: Katherine Harborne
Maintainer: <katherine.harborne@icrar.org>
Description: The purpose of this package is to provide a set of tools to "observe" simulations of galaxies as if you were an observer using an integral field unit (IFU). The galaxy model can be observed at a chosen inclination and an IFU observation can be produced in a selection of different spatial shapes (square, circular, and hexagonal) to account for the variety of CCD arrangements used by current telescopes. This data cube can be used to calculate the value an observer would measure for the spin parameter compared to the value that is measured directly from the simulation.
Expand Down
10 changes: 7 additions & 3 deletions NEWS.md
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
# SimSpin v2.3.10 News
# SimSpin v2.3.11 News

### Author: Kate Harborne

### Last edit: 27/07/22
### Last edit: 28/07/22

Below is a table containing a summary of all changes made to SimSpin, since the date this file was created on 26/08/2021.

Expand All @@ -18,7 +18,11 @@ All changes are noted in the changelog table below.

| Date | Summary of change | Version | Commit |
|----------|----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|---------|------------------------------------------|
| 27/07/22 | *Minor bugfix.* Update to `make_simspin_file` to increase the level of compression for output files. This will have no change on functions down the line. Also, fixed the units for `plot_h4` which was writing "km/s" when there are no units for h4. | 2.3.10 | |
| 28/07/22 | *Minor bugfix*. | 2.3.11 | |
| | (1) Update to `write_simspin_FITS`. The `observation` summary output by `build_datacube` was previously never saved to the FITS file output. However, this makes it difficult to reproduce if the only saved output is a FITS file. With this update, the list `observation` is now written to a table in EXTNUM = 3 for all FITS files produced. | | |
| | (2) Update to `method = "gas"`. Particle images will now be output when method is `gas`, along with both the raw and observed mass, velocity and dispersion maps. This mode was not updated with the `spectral` and `velocity` modes in v2.3.8 upgrade. This has now been remedied. | | |
| | (3) Addition of plotting functions for the gas mode specific images: `plot_SFR` and `plot_OH`. | | |
| 27/07/22 | *Minor bugfix.* Update to `make_simspin_file` to increase the level of compression for output files. This will have no change on functions down the line. Also, fixed the units for `plot_h4` which was writing "km/s" when there are no units for h4. | 2.3.10 | be0c281cb69f7352d3a88f59f7777bbffdc149d6 |
| 22/07/22 | **New feature!** Update to add plotting functions. A suite of functions with standard colour maps for each image type: `plot_flux`, `plot_mass`, `plot_velocity`, `plot_dispersion`, `plot_h3`, `plot_h4`, `plot_age`, `plot_metallicity` and `plot_particles`. Axis names and colour bars can be modified by the user. All new functions are tested and documented. | 2.3.9 | 2aceb04d6e40196168de13c4bbd450280de46c4d |
| 21/07/22 | *Minor bugfix*. With update to v2.3.5, flux conservation was implemented in mode="velocity". During this upgrade, filters were properly implemented meaning that re-run observations will now produce flux maps = 0 where previously they had flux. This occurs when the filter wavelength range does not overlap with the observed telescope wavelength range. | 2.3.8 | e043f8302ea8f5d0ac189a44edf33b6ceb866f50 |
| | (1) Errors now issued by `telescope` if you request a filter that will not allow any flux through the telescope wavelength range. | | |
Expand Down
16 changes: 15 additions & 1 deletion R/build_datacube.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -382,7 +382,8 @@ build_datacube = function(simspin_file, telescope, observing_strategy,
dispersion_image = array(data = output[[4]], dim = c(observation$sbin, observation$sbin)),
SFR_image = array(data = output[[5]], dim = c(observation$sbin, observation$sbin)),
metallicity_image = array(data = output[[6]], dim = c(observation$sbin, observation$sbin)),
OH_image = array(data = output[[7]], dim = c(observation$sbin, observation$sbin))
OH_image = array(data = output[[7]], dim = c(observation$sbin, observation$sbin)),
particle_image = array(data = output[[8]], dim = c(observation$sbin, observation$sbin))
)

output = list("velocity_cube" = cube,
Expand All @@ -400,12 +401,25 @@ build_datacube = function(simspin_file, telescope, observing_strategy,
output$observed_images$mass_image = array(0.0, dim = dims[c(1,2)])
output$observed_images$velocity_image = array(0.0, dim = dims[c(1,2)])
output$observed_images$dispersion_image = array(0.0, dim = dims[c(1,2)])
output$observed_images$h3_image = array(0.0, dim = dims[c(1,2)])
output$observed_images$h4_image = array(0.0, dim = dims[c(1,2)])

for (c in 1:dims[1]){
for (d in 1:dims[2]){
output$observed_images$mass_image[c,d] = sum(output$velocity_cube[c,d,])
output$observed_images$velocity_image[c,d] = .meanwt(observation$vbin_seq, output$velocity_cube[c,d,])
output$observed_images$dispersion_image[c,d] = sqrt(.varwt(observation$vbin_seq, output$velocity_cube[c,d,], output$observed_images$velocity_image[c,d]))

h3h4 = tryCatch({stats::optim(par = c(0,0),
fn = .losvd_fit,
x = observation$vbin_seq,
losvd = (output$velocity_cube[c,d,]/(max(output$velocity_cube[c,d,], na.rm=T))),
vel = output$observed_images$velocity_image[c,d],
sig = output$observed_images$dispersion_image[c,d],
method="BFGS", control=list(reltol=1e-9))$par},
error = function(e){c(0,0)})
output$observed_images$h3_image[c,d] = h3h4[1]
output$observed_images$h4_image[c,d] = h3h4[2]
}
}
}
Expand Down
138 changes: 138 additions & 0 deletions R/plot_images.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -619,6 +619,144 @@ plot_particles <- function(particle_image, fig = c(0,1,0,1), new=F,

}

#' Plotting pretty gas SFR images
#'
#' A function to produce a plot of the gas star formation rate image produced by
#' \code{build_datacube()} with associated colour bar and labels.
#'
#' @param SFR_image Numeric array containing the particle per pixel image
#' @param fig Numeric array of length 4 describing the boundary of the image
#' @param new Boolean. Should the image be added to the existing plot? Default
#' is FALSE.
#' @param units String describing the units of the values contained
#' in the image.
#' @param main Image title, default "".
#' @param radii list - containing a, b, and ang (if wishing to plot half-mass
#' radii ellipse).
#' @param na.color String. Colour given to NA values in the image.
#' @param zlim Numeric array of length 2. Describing the numeric range of
#' colours in the image. Default is NA, in which the range will be described by
#' the minimum and maximum values in the image.
#' @param legend Boolean to determine if the colour bar axis should be printed
#' at the bottom of the image. Default is T.
#' @param titleshift Numeric. Describes the distance between the colour bar and
#' the units. Default is -4.
#' @param labN Numeric. Describes the minimum number of numeric labels added to
#' the colour bar. Default is 5.
#' @param ... Further variables passed to magimage. See
#' \code{\link[magicaxis]{magimage}} for further details.
#' @return Returns an image to the plotting window of the input
#' \code{build_datacube} image.
#' @examples
#' ss_pd_eagle = system.file("extdata", "SimSpin_example_EAGLE.hdf5",
#' package = "SimSpin")
#' ss_eagle = make_simspin_file(ss_pd_eagle, write_to_file = FALSE)
#' cube = build_datacube(simspin_file = ss_eagle,
#' telescope = telescope(type="SAMI"),
#' observing_strategy = observing_strategy(),
#' method = "gas")
#' plot_SFR(cube$raw_images$SFR_image)

plot_SFR <- function(SFR_image, fig = c(0,1,0,1), new=F,
units = expression("sSFR, M"["sol"]*"/Gyr"), main="", radii = NA,
na.color = "white", zlim = NA, legend=T,
titleshift = -4, labN=5, ...){

part_map = SFR_image
im_dim = dim(part_map)%/%2
part_val = c(floor(min(part_map, na.rm=T)/2), max(part_map, na.rm=T))

if (all(part_val == 0)){
stop("Image contains only '0'. No image can be produced in this case. \n
Please check your build_datacube function and try again.")
}

part_map_cols = cmocean::cmocean("deep", version = "2.0", start = .1, end=1)(50)

par(pty="s", fig=fig, xpd=FALSE, ps=12, cex=1, new=new); options(scipen = 1)
.image_nan(z = part_map, zlim = if(is.na(zlim[1])){part_val}else{zlim}, col = part_map_cols, na.color = na.color, xaxt="n",
yaxt="n", ann=FALSE, magmap=FALSE, family="mono", font=1, main=main, ...)
if (!is.na(radii[1])){
plotrix::draw.ellipse(im_dim[1], im_dim[2], radii$a, radii$b, radii$ang, border = "red", density = NULL)
}
if (legend){
.magcolbar(position = "bottom", range = if(is.na(zlim[1])){part_val}else{zlim}, scale = c(1, 1/20),
col = part_map_cols, orient = "h", inset = -1/20, labN=labN, title = units,
titleshift = titleshift)
}

}

#' Plotting pretty log10(O/H) images
#'
#' A function to produce a plot of the gas log10(O/H) image produced by
#' \code{build_datacube()} with associated colour bar and labels.
#'
#' @param OH_image Numeric array containing the metallicity image.
#' @param fig Numeric array of length 4 describing the boundary of the image
#' @param new Boolean. Should the image be added to the existing plot? Default
#' is FALSE.
#' @param units String describing the units of the values contained
#' in the image.
#' @param main Image title, default "".
#' @param radii list - containing a, b, and ang (if wishing to plot half-mass
#' radii ellipse).
#' @param na.color String. Colour given to NA values in the image.
#' @param zlim Numeric array of length 2. Describing the numeric range of
#' colours in the image. Default is NA, in which the range will be described by
#' the minimum and maximum values in the image.
#' @param legend Boolean to determine if the colour bar axis should be printed
#' at the bottom of the image. Default is T.
#' @param titleshift Numeric. Describes the distance between the colour bar and
#' the units. Default is -4.
#' @param labN Numeric. Describes the minimum number of numeric labels added to
#' the colour bar. Default is 5.
#' @param ... Further variables passed to magimage. See
#' \code{\link[magicaxis]{magimage}} for further details.
#' @return Returns an image to the plotting window of the input
#' \code{build_datacube} image.
#' @examples
#' ss_pd_eagle = system.file("extdata", "SimSpin_example_EAGLE.hdf5",
#' package = "SimSpin")
#' ss_eagle = make_simspin_file(ss_pd_eagle, write_to_file = FALSE)
#' cube = build_datacube(simspin_file = ss_eagle,
#' telescope = telescope(type="SAMI"),
#' observing_strategy = observing_strategy(),
#' method = "gas")
#' plot_OH(cube$raw_images$OH_image)

plot_OH <- function(OH_image, fig = c(0,1,0,1), new=F,
units = expression("log10(O/H) + 12"), main="",
na.color = "white", zlim = NA, legend=T, radii = NA,
titleshift = -4, labN=5, ...){

met_map = OH_image
im_dim = dim(met_map)%/%2
met_val = c(floor(min(met_map, na.rm=T)/2), max(met_map, na.rm=T))

if (all(met_val == 0)){
stop("Image contains only '0'. No image can be produced in this case. \n
Please check your build_datacube function and try again.")
}

met_map_cols = cmocean::cmocean("dense", version = "2.0", start = .1, end=1)(50)

par(pty="s", fig=fig, xpd=FALSE, ps=12, cex=1, new=new); options(scipen = 1)
.image_nan(z = met_map, zlim = if(is.na(zlim[1])){met_val}else{zlim},
col = met_map_cols, na.color = na.color, xaxt="n",
yaxt="n", ann=FALSE, magmap=FALSE, family="mono", font=1, main=main, ...)
if (!is.na(radii[1])){
plotrix::draw.ellipse(im_dim[1], im_dim[2], radii$a, radii$b, radii$ang, border = "red", density = NULL)
}
if (legend){
.magcolbar(position = "bottom", range = if(is.na(zlim[1])){met_val}else{zlim}, scale = c(1, 1/20),
col = met_map_cols, orient = "h", inset = -1/20, labN=labN, title = units,
titleshift = titleshift)
}

}


.image_nan <- function(z, zlim, col, na.color='gray', ...){
# Function for plotting an image with NAs as a set colour
# z :==: array corresponding to the image to be plotted
Expand Down
15 changes: 11 additions & 4 deletions R/utilities.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1442,10 +1442,13 @@ globalVariables(c(".N", ":=", "Age", "Carbon", "CellSize", "Density", "Hydrogen"
SFR_map = array(data = 0.0, dim = observation$sbin^2)
Z_map = array(data = 0.0, dim = observation$sbin^2)
OH_map = array(data = 0.0, dim = observation$sbin^2)
part_map = array(data=0, dim = observation$sbin^2)

for (i in 1:(dim(part_in_spaxel)[1])){

num_part = part_in_spaxel$N[i] # number of particles in spaxel
part_map[part_in_spaxel$pixel_pos[i]] = num_part

galaxy_sample = galaxy_data[ID %in% part_in_spaxel$val[[i]]]

# adding the "gaussians" of each particle to the velocity bins
Expand All @@ -1461,7 +1464,7 @@ globalVariables(c(".N", ":=", "Age", "Carbon", "CellSize", "Density", "Hydrogen"

}

return(list(vel_spec, mass_map, vel_los, dis_los, SFR_map, Z_map, OH_map))
return(list(vel_spec, mass_map, vel_los, dis_los, SFR_map, Z_map, OH_map, part_map))
}

.gas_velocity_spaxels_mc = function(part_in_spaxel, observation, galaxy_data, simspin_data, verbose, cores){
Expand All @@ -1473,14 +1476,17 @@ globalVariables(c(".N", ":=", "Age", "Carbon", "CellSize", "Density", "Hydrogen"
SFR_map = array(data = 0.0, dim = observation$sbin^2)
Z_map = array(data = 0.0, dim = observation$sbin^2)
OH_map = array(data = 0.0, dim = observation$sbin^2)
part_map = array(data=0, dim = observation$sbin^2)

doParallel::registerDoParallel(cores)

i = integer()
output = foreach(i = 1:(dim(part_in_spaxel)[1]), .combine='.comb', .multicombine=TRUE,
.init=list(list(), list(), list(), list(), list(), list(), list())) %dopar% {
.init=list(list(), list(), list(), list(), list(), list(), list(), list())) %dopar% {

num_part = part_in_spaxel$N[i]
part_map = num_part

galaxy_sample = galaxy_data[ID %in% part_in_spaxel$val[[i]]]

# adding the "gaussians" of each particle to the velocity bins
Expand All @@ -1492,7 +1498,7 @@ globalVariables(c(".N", ":=", "Age", "Carbon", "CellSize", "Density", "Hydrogen"
Z_map = log10(mean(galaxy_sample$Metallicity)/0.0127)
OH_map = log10(mean(galaxy_sample$Oxygen/galaxy_sample$Hydrogen))+12

result = list(vel_spec, mass_map, vel_los, dis_los, SFR_map, Z_map, OH_map)
result = list(vel_spec, mass_map, vel_los, dis_los, SFR_map, Z_map, OH_map, part_map)
return(result)
closeAllConnections()
}
Expand All @@ -1504,8 +1510,9 @@ globalVariables(c(".N", ":=", "Age", "Carbon", "CellSize", "Density", "Hydrogen"
SFR_map[part_in_spaxel$pixel_pos] = matrix(unlist(output[[5]]))
Z_map[part_in_spaxel$pixel_pos] = matrix(unlist(output[[6]]))
OH_map[part_in_spaxel$pixel_pos] = matrix(unlist(output[[7]]))
part_map[part_in_spaxel$pixel_pos] = matrix(unlist(output[[8]]))

return(list(vel_spec, mass_map, vel_los, dis_los, SFR_map, Z_map, OH_map))
return(list(vel_spec, mass_map, vel_los, dis_los, SFR_map, Z_map, OH_map, part_map))

}

Expand Down

0 comments on commit c7114e4

Please sign in to comment.