add normalize dmass

Author: ibarraespinosa

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: rtorf
 Type: Package
 Title: 'rtorf' Tools for Observations, Receptors and Footprints
-Version: 2.5.0
+Version: 2.5.1
 Authors@R: c( 
   person(given = "Sergio", family = "Ibarra-Espinosa",
          role = c("aut", "cre"),
@@ -20,7 +20,7 @@ URL: https://github.com/noaa-gml/rtorf,
 BugReports: https://github.com/noaa-gml/rtorf/issues/
 Encoding: UTF-8
 RoxygenNote: 7.3.2
-Date: 2025-09-24
+Date: 2025-09-25
 Suggests: 
     knitr,
     rmarkdown,

NAMESPACE

@@ -21,6 +21,7 @@ export(obs_footname)
 export(obs_format)
 export(obs_format_et)
 export(obs_freq)
+export(obs_grid)
 export(obs_hysplit_ascdata)
 export(obs_hysplit_control)
 export(obs_hysplit_control_read)
@@ -42,6 +43,7 @@ export(obs_meta)
 export(obs_month_middle)
 export(obs_nc)
 export(obs_nc_get)
+export(obs_normalize_dmass)
 export(obs_out)
 export(obs_plot)
 export(obs_rbind)

NEWS.md

@@ -1,5 +1,8 @@
 NEWS
 ===========
+### rtorf 2.5.1 (Release date: 2025-09-25)
+
+add obs_normalize_dmass
 
 ### rtorf 2.5.0 (Release date: 2025-09-24)
 

R/10_obs_hysplit.R

@@ -926,3 +926,43 @@ obs_hysplit_control_read <- function(
   )
   return(l_control)
 }
+
+
+#' obs_grid
+#' Create grid for footprints
+#'
+#' @param numpix.x number of pixels in x direction
+#' @param numpix.y number of pixels in y direction
+#' @param lon.ll lower left corner of grid
+#' @param lat.ll lower left corner of grid
+#' @param lon.res resolution in degrees longitude
+#' @param lat.res resolution in degrees latitude
+#' @return A grid
+#' @note name it is not obs_hysplit_grid because we can create define grid in CONTROL
+#' @export
+#' @examples {
+#' # Do not run
+#' obs_hysplit_grid(control_file)
+#'
+#' }
+#
+obs_grid <- function(
+  numpix.x = 70,
+  numpix.y = 70,
+  lon.ll = 10,
+  lat.ll = -30,
+  lon.res = 1,
+  lat.res = 1
+) {
+  flat <- lat.ll + (lat.res) * (0:(numpix.y - 1))
+  flon <- lon.ll + (lon.res) * (0:(numpix.x - 1))
+
+  #get grid indices
+  #For horizontal grids (lower left corner of south-west gridcell: 11N,145W; resolution: 1/4 lon, 1/6 lat, 376 (x) times 324 (y))
+  lats <- seq(lat.ll, (lat.ll + (numpix.y - 1) * lat.res), lat.res)
+  lons <- seq(lon.ll, (lon.ll + (numpix.x - 1) * lon.res), lon.res)
+
+  #assumes westernmost gridcell is smallest longitude & easternmost
+  #gridcell is largest longitude; does not work over dateline
+  gitx <- floor(1 / lon.res * (part[, "lon"] - lon.ll) + 1)
+}

R/helpers_legacy.R

@@ -2773,7 +2773,7 @@ obs_make_hysplit_nc4 <- function(
 #' @title obs_traj_foot
 #' @family helpers legacy
 #' @name obs_traj_foot
-#' @description return time
+#' @description return trajectory
 #' @param ident  character value specifying the trajectory ensemble to look at
 #' @param part  object containing particle run; if NULL, then determine from ident and pathname
 #' @param pathname  path where object with particle locations is saved
@@ -2825,7 +2825,7 @@ obs_traj_foot <- function(
     #Check if object exists
     if (read.r) {
       if (existsr(ident, pathname)) {
-        print(paste("obs_traj_foot(): starting with ident=", ident, sep = "")) #found object
+        print(paste0("obs_traj_foot(): starting with ident=", ident)) #found object
         part <- getr(ident, pathname) #get it
       } else {
         #if not there, break out of function, return NA
@@ -3188,3 +3188,70 @@ obs_traj_foot <- function(
   dimnames(foot.arr) <- list(lats, lons, foottimes[1:(length(foottimes) - 1)])
   return(foot.arr)
 }
+
+
+#' @title obs_normalize_dmass
+#' @family helpers legacy
+#' @name obs_normalize_dmass
+#' @description add columns of
+#' @param part  particle data.table
+#' @return return footprint
+#' @export
+#' @examples {
+#' \dontrun{
+#' # Do not run
+#' }}
+obs_normalize_dmass <- function(
+  part
+) {
+  part <- data.table::as.data.table(part)
+
+  # check if btime is not presnet and add it
+  if (!"btime" %in% names(part)) {
+    cat("Adding btime\n")
+    time <- NULL
+    part[, btime := abs(time) / 60]
+  }
+
+  # must be in this order (might not be...)
+  # this is from lower time to bigger time by part index
+  # index from first to last repeating btime
+  data.table::setorderv(part, c("btime", "index")) #data.table
+
+  # identify particles with too strong dmass violation
+  part[
+    dmass < 1 / 1e3 |
+      dmass > 1e3,
+    unique(index)
+  ] -> ind
+
+  if (length(ind) >= length(part[, unique(index)]) / 2) {
+    stop(cat(
+      "More than 50% of particles have mass defect\n"
+    ))
+  }
+
+  # for those particles with strong dmass violation set dmass entries to NA
+  part[
+    dmass < 1 / 1e3 |
+      dmass > 1e3,
+    dmass := NA
+  ]
+
+  part[,
+    mean_dmass := mean(dmass, na.rm = T),
+    by = .(btime)
+  ]
+
+  # if mean_dmass is 0
+
+  part[mean_dmass == 0, mean_dmass := 0.0001] # legacy criteria
+
+  part[is.na(dmass), dmass := mean_dmass]
+
+  part[, "dmass"] <- part[, "dmass"] / mean.dmass
+
+  # normalize
+  part[, ndmass := dmass / mean_dmass]
+  return(part)
+}