Merge pull request #90 from kateharborne/dev-Voronoi

dev-Voronoi
kateharborne · Oct 26, 2023 · 2df3fc5 · 2df3fc5
2 parents 76fd5c2 + e51e73e
commit 2df3fc5
Show file tree

Hide file tree

Showing 30 changed files with 1,497 additions and 292 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: SimSpin
 Type: Package
 Title: SimSpin - A package for the kinematic analysis of galaxy simulations
-Version: 2.6.0
+Version: 2.7.0
 Author: Katherine Harborne
 Co-author: Alice Serene
 Maintainer: <katherine.harborne@icrar.org>
@@ -32,7 +32,7 @@ Imports:
     stringr,
     stats,
     utils
-RoxygenNote: 7.2.2
+RoxygenNote: 7.2.3
 Suggests: 
     testthat,
     knitr,

diff --git a/NEWS.md b/NEWS.md
@@ -1,6 +1,6 @@
-# SimSpin v2.6.0 News
+# SimSpin v2.7.0 News
 
-### Last edit: 03/08/23
+### Last edit: 25/10/23
 
 Below is a table containing a summary of all changes made to SimSpin, since the date this file was created on 26/08/2021.
 
@@ -16,6 +16,7 @@ All changes are noted in the changelog table below.
 
 | Date     	| Summary of change                                                                                                                                                                                                                                                                                                                                                                                                        	| Version 	| Commit                                   	| Author            |
 |----------	|--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------	|---------	|------------------------------------------	| ----------------- |
+| 25/10/23 | Adding functionality to voronoi bin individual pixels until the contain some minimum number of particles. This is important for reducing the effects of shot noise on the measured kinematics, especially in dispersion supported systems. We add a new function `voronoi`, which performs the tessellation algorithm as used in the [Vorbin python package](https://pypi.org/project/vorbin/) (Cappellari & Copin, 2003). When requested, this will add a new image to the `raw_images` output list called `voronoi_bins`, detailing the bin ids assigned to each pixel in the field-of-view of a given observation. | 2.7.0 |  | Kate Harborne |
 | 03/08/23 | Updating the information stored in the SimSpin files following a suggestion from SimSpin v2.5.0 paper review. Do not need to store the full spectrum for each particle. Just the weights and id's of the spectra within the template set. Reduces the simspin file sizes by factor of 100. Old SimSpin files will still work with the code, but will now issue a warning when using a file older than v2.6.0. Tagging for release. Thank you anonymous reviewer!!! | 2.6.0 | ecb9d38031b40e63216d6200409b24ece8fadb69 | Kate Harborne |
 | 19/07/23 | Tagging a stable release of the code to align with the submission of the SimSpin v2.5.0 paper to PASA on the 28th June 2023. | 2.5.0 | cc3b2018527995dfaf09130b1afbe4341566e4f7 | Kate Harborne |
 | 15/06/23 | Updating the noise implementation for `build_datacube` using the median flux pixel to scale the level of noise such that the S/N does not go extreme towards the peak flux or minimum flux in the image. | 2.4.10 | 3bf98efea656776436fca3384cda1411f8f9e642 | Kate Harborne |

diff --git a/R/build_datacube.R b/R/build_datacube.R
@@ -48,6 +48,11 @@
 #' and velocity gridding on. Default is 1.
 #'@param mass_flag Boolean flag that, when set to TRUE, will compute properties
 #' using a mass weighting rather than a luminosity weighting. Default is FALSE.
+#'@param voronoi_bin Boolean flag that, when set to TRUE, will bin pixels into
+#' voronoi tessellated cells that contain a minimum number of particles per
+#' pixel, specified by \code{vorbin_limit}. Default is FALSE.
+#'@param vorbin_limit Integer float that describes the minimum number of
+#' particles per pixel within a given bin, only used if \code{voronoi_bin = T}.
 #'@return Returns a list containing four elements:
 #'\enumerate{
 #' \item \code{spectral_cube} or \code{velocity_cube} - a 3D array containing
@@ -80,7 +85,8 @@ build_datacube = function(simspin_file, telescope, observing_strategy,
                           telescope_name="SimSpin",
                           observer_name="Anonymous",
                           split_save=F,
-                          cores=1, mass_flag = F){
+                          cores=1, mass_flag = F,
+                          voronoi_bin=F, vorbin_limit=10){
 
   if (missing(method)){
     if ("method" %in% names(telescope)){
@@ -234,6 +240,19 @@ build_datacube = function(simspin_file, telescope, observing_strategy,
   # which particles sit in each spaxel?
   part_in_spaxel = galaxy_data[, list(val=list(ID), .N), by = "pixel_pos"]
 
+  if (voronoi_bin){ # Returning the binned pixels based on some voronoi limit
+    if (verbose){cat("Binning spaxels into voronoi bins... \n")}
+
+    part_in_spaxel = voronoi(part_in_spaxel=part_in_spaxel, obs=observation,
+                             particle_limit = as.numeric(vorbin_limit),
+                             roundness_limit = 0.3, uniform_limit = 0.8,
+                             verbose = verbose)
+
+    observation$particle_limit = as.numeric(vorbin_limit)
+
+  }
+
+
   # SPECTRAL mode method =======================================================
   if (observation$method == "spectral"){
 
@@ -289,7 +308,8 @@ build_datacube = function(simspin_file, telescope, observing_strategy,
       dispersion_image = array(data = output[[4]], dim = c(observation$sbin, observation$sbin)),
       age_image = array(data = output[[5]], dim = c(observation$sbin, observation$sbin)),
       metallicity_image = array(data = output[[6]], dim = c(observation$sbin, observation$sbin)),
-      particle_image = array(data = output[[7]], dim = c(observation$sbin, observation$sbin))
+      particle_image = array(data = output[[7]], dim = c(observation$sbin, observation$sbin)),
+      voronoi_bins = array(data = output[[8]], dim = c(observation$sbin, observation$sbin))
       )
 
     output = list("spectral_cube"    = cube,
@@ -358,7 +378,8 @@ build_datacube = function(simspin_file, telescope, observing_strategy,
       age_image = array(data = output[[6]], dim = c(observation$sbin, observation$sbin)),
       metallicity_image = array(data = output[[7]], dim = c(observation$sbin, observation$sbin)),
       mass_image = array(data = output[[8]], dim = c(observation$sbin, observation$sbin)),
-      particle_image = array(data = output[[9]], dim = c(observation$sbin, observation$sbin))
+      particle_image = array(data = output[[9]], dim = c(observation$sbin, observation$sbin)),
+      voronoi_bins = array(data = output[[10]], dim = c(observation$sbin, observation$sbin))
       )
     observed_images = list(
       flux_image = array(data = output[[3]], dim = c(observation$sbin, observation$sbin)),
@@ -465,7 +486,8 @@ build_datacube = function(simspin_file, telescope, observing_strategy,
       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)),
-      particle_image = array(data = output[[8]], dim = c(observation$sbin, observation$sbin))
+      particle_image = array(data = output[[8]], dim = c(observation$sbin, observation$sbin)),
+      voronoi_bins = array(data = output[[9]], dim = c(observation$sbin, observation$sbin))
       )
     observed_images = list(
       flux_image = array(data = output[[2]], dim = c(observation$sbin, observation$sbin)),
@@ -538,6 +560,10 @@ build_datacube = function(simspin_file, telescope, observing_strategy,
 
   }
 
+  if (!voronoi_bin){ # if vorbin has not been requested, don't supply it in the output
+    output$raw_images = output$raw_images[-which(names(output$raw_images) == "voronoi_bins")]
+  }
+
   # Trimming off extra zeros from images outside the aperture of the telescope
   aperture_region = matrix(data = observation$aperture_region, nrow = observation$sbin, ncol = observation$sbin)
 

diff --git a/R/observation.R b/R/observation.R
@@ -84,6 +84,7 @@ observation = function(telescope, observing_strategy, method){
                 lum_dist        = lum_dist,
                 method          = method,
                 origin          = paste0("SimSpin_v", packageVersion("SimSpin")),
+                particle_limit  = 0,
                 pointing_kpc    = xy_kpc(observing_strategy$pointing),
                 pointing_deg    = xy_deg(observing_strategy$pointing),
                 pixel_region    = aperture_region * pixel_index,