Skip to content

Commit

Permalink
Merge pull request #100 from kateharborne/dev-obs_mass
Browse files Browse the repository at this point in the history 
dev_obs-mass
  • Loading branch information
@kateharborne
kateharborne authored Mar 14, 2024
2 parents 0469ff4 + 00f0ea7 commit 2e878f7
Show file tree
Hide file tree
Showing 7 changed files with 70 additions and 46 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.8.2
Version: 2.8.3
Author: Katherine Harborne
Co-author: Alice Serene
Maintainer: <katherine.harborne@icrar.org>
Expand Down
7 changes: 4 additions & 3 deletions NEWS.md
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
# SimSpin v2.8.2 News
# SimSpin v2.8.3 News

### Last edit: 29/02/2024
### Last edit: 13/03/2024

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 @@ -16,7 +16,8 @@ All changes are noted in the changelog table below.

| Date | Summary of change | Version | Commit | Author |
|---------- |-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- |--------- |------------------------------------------ | ----------------- |
| 29/02/24 | *Bug fix.* Updating `make_simspin_file` to account for the `Time` at which the input simulation was extracted (for proper assignment of stellar ages relative to the snapshot, not z=0). Fixing issue #98.| 2.8.2 | | Kate Harborne |
| 13/03/24 | Adding observation-ally mapped mass and SFR images to the output `$observed_images` - the raw images are not blurred with the PSF, but for comparison with other observed properties, this PSF convolution is necessary. Now mass maps are output on all `method="velocity"` cubes as default and SFR maps on all `method="gas"` or `"sf gas"`. | 2.8.3 | | Kate Harborne |
| 29/02/24 | *Bug fix.* Updating `make_simspin_file` to account for the `Time` at which the input simulation was extracted (for proper assignment of stellar ages relative to the snapshot, not z=0). Fixing issue #98.| 2.8.2 | 0469ff454756e7b6195bf927b7c867fd02f20743 | Kate Harborne |
| 10/01/24 | Also updating treatment for `mass_flag` = TRUE. No longer will the `observed_images` outputs change for the flag, the flux map will always be included. However, the kinematic maps observed will be weighted differently (by mass instead of the flux), which is now recorded in the `observation` log and output to FITS with a new header unit denoting mass rather than flux weighting. | 2.8.1 | 352c5601b04767e10d2598358b164f1b203a8d54 | Kate Harborne |
| 10/01/24 | Standardising which cubes/images output are voronoi binned. Properties, such as the mass, flux and number of particles, are produced simply as the sum along that line-of-sight and as such should be presented at the pixel resolution of the instrument. As default, summed images are now given in this way, while any spectral/velocity cubes are binned along with any output average maps (kinematics such as velocity, dispersion, h3 and h4, metallicity and age). The same is done for gas maps, with all but the mass, SFR and particle number being given in binned form when the voronoi function is turned on. Also, added `raw_images$mass_image` to the output of spectral mode as this is now "free"" to produce. | 2.8.0 | 5af1ad3a90c7f84a55e6f9cb50492308fd422d7a | Kate Harborne |
| 17/11/23 | Bug fix. The mass/flux in voronoi binned pixels were displayed as the sum of all the particles that fell within that bin. Now updated so that this value is averaged across the area of the bin (else outer larger bins appear brighter than inner ones...). Fixing issue #91. | 2.7.1 | e93c763c40ec670da58b37a353586887d0b34143 | Kate Harborne |
Expand Down
19 changes: 16 additions & 3 deletions R/blur_datacube.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,9 +79,22 @@ blur_datacube = function(datacube_output){
}

# 2. Blurring the observed flux map
blur_image = array(data = 0.0, dim = cube_dims[c(1,2)])
blur_image = ProFit::profitBruteConv(datacube_output$observed_images$flux_image, observation$psf_kernel) * aperture_region
datacube_output$observed_images$flux_image = blur_image
if (observation$method == "velocity"){ # if your observation is of stars, blur the flux map
blur_image = array(data = 0.0, dim = cube_dims[c(1,2)])
blur_image = ProFit::profitBruteConv(datacube_output$observed_images$flux_image, observation$psf_kernel) * aperture_region
datacube_output$observed_images$flux_image = blur_image
}

blur_mass_image = array(data = 0.0, dim = cube_dims[c(1,2)])
blur_mass_image = ProFit::profitBruteConv(datacube_output$observed_images$mass_image, observation$psf_kernel) * aperture_region
datacube_output$observed_images$mass_image = blur_mass_image

if (observation$method == "gas" | observation$method == "sf gas"){ # if your observation is of gas, blur the sfr map
blur_sfr_image = array(data = 0.0, dim = cube_dims[c(1,2)])
blur_sfr_image = ProFit::profitBruteConv(datacube_output$observed_images$SFR_image, observation$psf_kernel) * aperture_region
datacube_output$observed_images$SFR_image = blur_sfr_image
}


# Returning output in same format as input
blur_output = list("velocity_cube" = blur_cube,
Expand Down
10 changes: 5 additions & 5 deletions R/build_datacube.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -432,7 +432,8 @@ build_datacube = function(simspin_file, telescope, observing_strategy,
dispersion_image = array(0.0, dim = c(observation$sbin, observation$sbin)),
h3_image = array(0.0, dim = c(observation$sbin, observation$sbin)),
h4_image = array(0.0, dim = c(observation$sbin, observation$sbin)),
residuals = array(0.0, dim = c(observation$sbin, observation$sbin))
residuals = array(0.0, dim = c(observation$sbin, observation$sbin)),
mass_image = array(data = summed_images$mass, dim = c(observation$sbin, observation$sbin))
)

output = list("velocity_cube" = cube,
Expand Down Expand Up @@ -532,12 +533,13 @@ build_datacube = function(simspin_file, telescope, observing_strategy,
voronoi_bins = array(data = output[[6]], dim = c(observation$sbin, observation$sbin))
)
observed_images = list(
flux_image = array(data = summed_images$mass, dim = c(observation$sbin, observation$sbin)),
mass_image = array(data = summed_images$mass, dim = c(observation$sbin, observation$sbin)),
velocity_image = array(0.0, dim = c(observation$sbin, observation$sbin)),
dispersion_image = array(0.0, dim = c(observation$sbin, observation$sbin)),
h3_image = array(0.0, dim = c(observation$sbin, observation$sbin)),
h4_image = array(0.0, dim = c(observation$sbin, observation$sbin)),
residuals = array(0.0, dim = c(observation$sbin, observation$sbin))
residuals = array(0.0, dim = c(observation$sbin, observation$sbin)),
SFR_image = array(data = summed_images$sfr, dim = c(observation$sbin, observation$sbin))
)

output = list("velocity_cube" = cube,
Expand Down Expand Up @@ -596,8 +598,6 @@ build_datacube = function(simspin_file, telescope, observing_strategy,
}
}

names(output$observed_images)[which(names(output$observed_images) == "flux_image")] = "mass_image"

if (verbose){cat("Done! \n")}

}
Expand Down
52 changes: 29 additions & 23 deletions R/write_simspin_FITS.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -532,39 +532,43 @@ write_simspin_FITS = function(output_file, simspin_datacube, object_name,
"EXTNAME"="Image extension name")

extnames = if (voronoi){
c("OBS_FLUX", "OBS_VEL", "OBS_DISP", "OBS_H3", "OBS_H4", "RESIDUAL", "RAW_FLUX", "RAW_MASS", "RAW_VEL", "RAW_DISP", "RAW_AGE", "RAW_Z", "NPART", "VORONOI")
c("OBS_FLUX", "OBS_VEL", "OBS_DISP", "OBS_H3", "OBS_H4", "RESIDUAL", "OBS_MASS", "RAW_FLUX", "RAW_MASS", "RAW_VEL", "RAW_DISP", "RAW_AGE", "RAW_Z", "NPART", "VORONOI")
} else {
c("OBS_FLUX", "OBS_VEL", "OBS_DISP", "OBS_H3", "OBS_H4", "RESIDUAL", "RAW_FLUX", "RAW_MASS", "RAW_VEL", "RAW_DISP", "RAW_AGE", "RAW_Z", "NPART")
c("OBS_FLUX", "OBS_VEL", "OBS_DISP", "OBS_H3", "OBS_H4", "RESIDUAL", "OBS_MASS", "RAW_FLUX", "RAW_MASS", "RAW_VEL", "RAW_DISP", "RAW_AGE", "RAW_Z", "NPART")
}

bunits = if (voronoi){
c("erg/s/cm**2", "km/s", "km/s", "unitless", "unitless", "percentage", "erg/s/cm**2", "Msol", "km/s", "km/s", "Gyr", "Z_solar", "Particle number", "Bin ID")
c("erg/s/cm**2", "km/s", "km/s", "unitless", "unitless", "percentage", "Msol", "erg/s/cm**2", "Msol", "km/s", "km/s", "Gyr", "Z_solar", "Particle number", "Bin ID")
} else {
c("erg/s/cm**2", "km/s", "km/s", "unitless", "unitless", "percentage", "erg/s/cm**2", "Msol", "km/s", "km/s", "Gyr", "Z_solar", "Particle number")
c("erg/s/cm**2", "km/s", "km/s", "unitless", "unitless", "percentage", "Msol", "erg/s/cm**2", "Msol", "km/s", "km/s", "Gyr", "Z_solar", "Particle number")
}

image_names = if (voronoi){
c("flux_image", "velocity_image", "dispersion_image", "h3_image", "h4_image", "residuals",
c("flux_image", "velocity_image", "dispersion_image", "h3_image", "h4_image",
"residuals", "mass_image",
"flux_image", "mass_image", "velocity_image", "dispersion_image", "age_image",
"metallicity_image", "particle_image", "voronoi_bins")
} else {
c("flux_image", "velocity_image", "dispersion_image", "h3_image", "h4_image", "residuals",
c("flux_image", "velocity_image", "dispersion_image", "h3_image", "h4_image",
"residuals", "mass_image",
"flux_image", "mass_image", "velocity_image", "dispersion_image", "age_image",
"metallicity_image", "particle_image")
}

rawobs = if (voronoi){
c("obs", "obs", "obs", "obs", "obs", "obs", "raw", "raw", "raw", "raw", "raw", "raw", "raw", "raw")
c("obs", "obs", "obs", "obs", "obs", "obs", "obs",
"raw", "raw", "raw", "raw", "raw", "raw", "raw", "raw")
} else {
c("obs", "obs", "obs", "obs", "obs", "obs", "raw", "raw", "raw", "raw", "raw", "raw", "raw")
c("obs", "obs", "obs", "obs", "obs", "obs", "obs",
"raw", "raw", "raw", "raw", "raw", "raw", "raw")
}

output_image_file_names = paste0(output_dir, "/", output_file_root, "_", rawobs, "_", image_names, ".FITS")

extnum = if (voronoi){
c(4,5,6,7,8,9,10,11,12,13,14,15,16,17)
c(4,5,6,7,8,9,10,11,12,13,14,15,16,17,18)
} else {
c(4,5,6,7,8,9,10,11,12,13,14,15,16)
c(4,5,6,7,8,9,10,11,12,13,14,15,16,17)
}

if (split_save){ # if writing each image to a seperate file
Expand All @@ -578,7 +582,7 @@ write_simspin_FITS = function(output_file, simspin_datacube, object_name,
image_keyvalues$BUNIT = bunits[i]
image_keyvalues$EXTNAME = extnames[i]

if (i < 7){ # 2. Write the image to this new file HDU 2
if (i < 8){ # 2. Write the image to this new file HDU 2
Rfits::Rfits_write_image(data = simspin_datacube$observed_images[[which(names(simspin_datacube$observed_images) == image_names[i])]],
filename = output_image_file_names[i], ext=2,
keyvalues = image_keyvalues, keycomments = image_keycomments,
Expand All @@ -598,7 +602,7 @@ write_simspin_FITS = function(output_file, simspin_datacube, object_name,
image_keyvalues$BUNIT = bunits[i]
image_keyvalues$EXTNAME = extnames[i]

if (i < 7){ # Write each subsequent image to the next HDU
if (i < 8){ # Write each subsequent image to the next HDU
Rfits::Rfits_write_image(data = simspin_datacube$observed_images[[which(names(simspin_datacube$observed_images) == image_names[i])]],
filename = cube_file_name, ext=extnum[i],
keyvalues = image_keyvalues, keycomments = image_keycomments,
Expand Down Expand Up @@ -718,51 +722,53 @@ write_simspin_FITS = function(output_file, simspin_datacube, object_name,

extnames =
if (voronoi){
c("OBS_MASS", "OBS_VEL", "OBS_DISP", "OBS_H3", "OBS_H4", "RESIDUAL",
c("OBS_MASS", "OBS_VEL", "OBS_DISP", "OBS_H3", "OBS_H4", "RESIDUAL", "OBS_SFR",
"RAW_MASS", "RAW_VEL", "RAW_DISP", "RAW_Z", "RAW_OH", "RAW_SFR", "NPART", "VORONOI")
} else {
c("OBS_MASS", "OBS_VEL", "OBS_DISP", "OBS_H3", "OBS_H4", "RESIDUAL",
c("OBS_MASS", "OBS_VEL", "OBS_DISP", "OBS_H3", "OBS_H4", "RESIDUAL", "OBS_SFR",
"RAW_MASS", "RAW_VEL", "RAW_DISP", "RAW_Z", "RAW_OH", "RAW_SFR", "NPART")
}

bunits =
if (voronoi){
c("Msol", "km/s", "km/s", "unitless", "unitless", "percentage",
c("Msol", "km/s", "km/s", "unitless", "unitless", "percentage", "Msol/year",
"Msol", "km/s", "km/s", "log10(Z/Z_solar)", "log10(O/H)+12", "Msol/year",
"Particle number", "Bin ID")
} else {
c("Msol", "km/s", "km/s", "unitless", "unitless", "percentage",
c("Msol", "km/s", "km/s", "unitless", "unitless", "percentage", "Msol/year",
"Msol", "km/s", "km/s", "log10(Z/Z_solar)", "log10(O/H)+12", "Msol/year",
"Particle number")
}

image_names =
if (voronoi){
c("mass_image", "velocity_image", "dispersion_image", "h3_image", "h4_image",
"residuals", "mass_image", "velocity_image", "dispersion_image", "metallicity_image",
"residuals", "SFR_image",
"mass_image", "velocity_image", "dispersion_image", "metallicity_image",
"OH_image", "SFR_image", "particle_image", "voronoi_bins")
} else {
c("mass_image", "velocity_image", "dispersion_image", "h3_image", "h4_image",
"residuals", "mass_image", "velocity_image", "dispersion_image", "metallicity_image",
"residuals", "SFR_image",
"mass_image", "velocity_image", "dispersion_image", "metallicity_image",
"OH_image", "SFR_image", "particle_image")
}

rawobs =
if (voronoi){
c("obs", "obs", "obs", "obs", "obs", "obs",
c("obs", "obs", "obs", "obs", "obs", "obs", "obs",
"raw", "raw", "raw", "raw", "raw", "raw", "raw", "raw")
} else {
c("obs", "obs", "obs", "obs", "obs", "obs",
c("obs", "obs", "obs", "obs", "obs", "obs", "obs",
"raw", "raw", "raw", "raw", "raw", "raw", "raw")
}

output_image_file_names = paste0(output_dir, "/", output_file_root, "_", rawobs, "_", image_names, ".FITS")

extnum =
if (voronoi){
c(4,5,6,7,8,9,10,11,12,13,14,15,16,17)
c(4,5,6,7,8,9,10,11,12,13,14,15,16,17,18)
} else {
c(4,5,6,7,8,9,10,11,12,13,14,15,16)
c(4,5,6,7,8,9,10,11,12,13,14,15,16,17)
}

if (split_save){ # if writing each image to a seperate file
Expand All @@ -775,7 +781,7 @@ write_simspin_FITS = function(output_file, simspin_datacube, object_name,
image_keyvalues$BUNIT = bunits[i]
image_keyvalues$EXTNAME = extnames[i]

if (i < 7){ # write observed mass, velocity and dispersion images to the file
if (i < 8){ # write observed mass, velocity and dispersion images to the file
Rfits::Rfits_write_image(data = simspin_datacube$observed_images[[which(names(simspin_datacube$observed_images) == image_names[i])]],
filename = output_image_file_names[i], ext=2,
keyvalues = image_keyvalues, keycomments = image_keycomments,
Expand Down
2 changes: 1 addition & 1 deletion README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@

<p>&nbsp;</p>

v2.8.2 - A package for producing mock observations:
v2.8.3 - A package for producing mock observations:

SimSpin allows you to take a simulation of a galaxy and produce a data cube in the style of an Integral Field Spectroscopy (IFS) instrument. You can find the live documentation for this code at the following [website](https://kateharborne.github.io/SimSpin/).

Expand Down
24 changes: 14 additions & 10 deletions tests/testthat/test_build_datacube.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -34,22 +34,24 @@ spectra_number_of_hdu_snfalse = 10
spectral_raw_vel_loc = 5

velocity_raw_images_size = 7
velocity_observed_images_size = 6
velocity_observed_images_size = 7

velocity_number_of_hdu_sntrue = 17
velocity_number_of_hdu_snfalse = 16
velocity_number_of_hdu_sntrue = 18
velocity_number_of_hdu_snfalse = 17
velocity_obs_vel_loc = 5
velocity_obs_h3_loc = 7
velocity_obs_h4_loc = 8
velocity_obs_res_loc = 9
velocity_raw_mass_loc = 11
velocity_variance_loc = 17
velocity_obs_mass_loc = 10
velocity_raw_mass_loc = 12
velocity_variance_loc = 18

gas_raw_images_size = 7
gas_observed_images_size = 6
gas_number_of_hdu_snfalse = 16
gas_number_of_hdu_sntrue = 17
gas_observed_images_size = 7
gas_number_of_hdu_snfalse = 17
gas_number_of_hdu_sntrue = 18
gas_obs_res_loc = 9
gas_obs_sfr_loc = 10

# Testing that build_datacube works in spectral mode ----

Expand Down Expand Up @@ -664,6 +666,7 @@ test_that("Data cubes can be written to a single files", {
expect_true(file.exists(paste0(temp_loc, "/ss_gadget_mft.FITS")))
expect_true(length(Rfits::Rfits_read(paste0(temp_loc, "/ss_gadget_mft.FITS"))) == velocity_number_of_hdu_sntrue)
expect_true(names(Rfits::Rfits_read(paste0(temp_loc, "/ss_gadget_mft.FITS")))[velocity_raw_mass_loc] == "RAW_MASS")
expect_true(names(Rfits::Rfits_read(paste0(temp_loc, "/ss_gadget_mft.FITS")))[velocity_obs_mass_loc] == "OBS_MASS")
expect_true(names(Rfits::Rfits_read(paste0(temp_loc, "/ss_gadget_mft.FITS")))[velocity_variance_loc] == "STAT")
expect_true(stringr::str_detect(Rfits::Rfits_read_header_raw(paste0(temp_loc, "/ss_gadget_mft.FITS")), "Msol"))

Expand Down Expand Up @@ -745,11 +748,12 @@ test_that("Data cubes can be written to a single files", {
observing_strategy = observing_strategy(dist_z = 0.03, inc_deg = 45, blur = T),
method="sf gas",
write_fits = T, output_location = paste0(temp_loc, "/ss_magneticum.FITS"),
split_save=F), built_cube_size)
split_save=F, voronoi_bin = T, vorbin_limit = 10), built_cube_size)

expect_true(file.exists(paste0(temp_loc, "/ss_magneticum.FITS")))
expect_true(length(Rfits::Rfits_read(paste0(temp_loc, "/ss_magneticum.FITS"))) == gas_number_of_hdu_snfalse)
expect_true(length(Rfits::Rfits_read(paste0(temp_loc, "/ss_magneticum.FITS"))) == (gas_number_of_hdu_snfalse+1))
expect_true(Rfits::Rfits_read(paste0(temp_loc, "/ss_magneticum.FITS"))[[gas_obs_res_loc]]$keyvalues$EXTNAME == "RESIDUAL")
expect_true(Rfits::Rfits_read(paste0(temp_loc, "/ss_magneticum.FITS"))[[gas_obs_sfr_loc]]$keyvalues$EXTNAME == "OBS_SFR")

expect_length(build_datacube(simspin_file = ss_hdf5,
telescope = telescope(type="IFU", lsf_fwhm = 3.6, signal_to_noise = 10, spatial_res = 0.25, fov=50),
Expand Down

0 comments on commit 2e878f7

Please sign in to comment.