fix Euclidian distance calculation

DESCRIPTION

@@ -1,8 +1,8 @@
 Type: Package
 Package: openair
 Title: Tools for the Analysis of Air Pollution Data
-Version: 2.18-1
-Date: 2023-10-20
+Version: 2.18-2
+Date: 2023-11-07
 Authors@R: c(
     person("David", "Carslaw", , "[email protected]", role = c("aut", "cre")),
     person("Jack", "Davison", , "[email protected]", role = "aut"),

NEWS.md

@@ -3,6 +3,7 @@
 ## Bug fixes
 
 - fix date formatting issue in `aqStats`.
+- fix wrong formula for Euclidian distances in `trajCluster` that did not transform coordinates before distance matrix was calculated. Thanks to Dan Jaffe.
 
 # openair 2.18-0
 

R/trajCluster.R

@@ -101,7 +101,7 @@ trajCluster <- function(traj, method = "Euclid", n.cluster = 5,
     mutate(traj_len = length(date))
 
   if (length(unique(traj$traj_len)) > 1) {
-    warning("Trajectory lengths differ, using most common length.")
+
     ux <- unique(traj$traj_len)
     nmax <- ux[which.max(tabulate(match(traj$traj_len, ux)))]
     traj <- ungroup(traj) %>%

src/cluster.cpp

@@ -1,54 +1,56 @@
-#include <Rcpp.h>
-using namespace Rcpp;
-
-// function to calculate Eucledian distance for trajectories
-RcppExport SEXP distEuclid(SEXP x, SEXP y, SEXP z){
-  NumericMatrix Am(x); // trajectory matrix of lats
-  NumericMatrix Bm(y); // trajectory matrix of lons
-  NumericMatrix Zm(z); // results matrix
-  double sum = 0.0;
-  int ncolumns = Am.ncol();
-  int nrows = Am.nrow();
-  for (int i = 0; i < ncolumns; i++) {
-    for (int j = 0; j < ncolumns; j++) {
-      sum = 0;
-      if (j >= i) { // only need lower triangle
-        for (int k = 0; k < nrows; k++) {
-          sum += pow(Am(k , i) - Am(k , j), 2)   + pow(Bm(k , i) - Bm(k , j), 2);
-        }
-        Zm(j, i) = sqrt(sum);
-      }
-    }
-  }
-  return Zm ;
-}
-
-
-// distance matrix based on similarity of trajectory angles
-RcppExport SEXP distAngle(SEXP x, SEXP y, SEXP z){
-  NumericMatrix Lonm(x); // matrix of longitudes
-  NumericMatrix Latm(y); // matrix of latitudes
-  NumericMatrix Zm(z);
-  double A, B, C; // intermediate trajectory variables
-  double X0, Y0; // origin of lon, lat
-  double sum = 0.0;
-  int ncolumns = Lonm.ncol();
-  int nrows = Lonm.nrow();
-  X0 = (double) Lonm(nrows - 1, 0);
-  Y0 = (double) Latm(nrows - 1, 0);
-  for (int i = 0; i < ncolumns; i++) {
-    for (int j = 0; j < ncolumns; j++) {
-      sum = 0;
-      if (j > i) { // only need lower triangle
-        for (int k = 0; k < nrows - 1; k++) { // do not include origin
-          A = pow((Lonm(k, i) - X0), 2) + pow((Latm(k, i) - Y0), 2);
-          B = pow((Lonm(k, j) - X0), 2) + pow((Latm(k, j) - Y0), 2);
-          C = pow((Lonm(k, j) - Lonm(k, i)), 2) + pow((Latm(k, j) - Latm(k, i)), 2);
-          sum += acos((0.5 * (A + B - C) / pow(A * B, 0.5)));
-        }
-        Zm(j, i) = sum / (nrows - 1) ;
-      }
-    }
-  }
-  return Zm;
-}
+#include <Rcpp.h>
+using namespace Rcpp;
+
+// function to calculate Eucledian distance for trajectories
+RcppExport SEXP distEuclid(SEXP x, SEXP y, SEXP z){
+  NumericMatrix Am(x); // trajectory matrix of lats
+  NumericMatrix Bm(y); // trajectory matrix of lons
+  NumericMatrix Zm(z); // results matrix
+  double sum = 0.0;
+  int ncolumns = Am.ncol();
+  int nrows = Am.nrow();
+  for (int i = 0; i < ncolumns; i++) {
+    for (int j = 0; j < ncolumns; j++) {
+      sum = 0;
+      if (j >= i) { // only need lower triangle
+        for (int k = 0; k < nrows; k++) {
+
+          sum += 110 * sqrt(pow(Am(k , i) - Am(k , j), 2) + 
+            pow((Bm(k , i) - Bm(k , j)) * cos(Am(k , i) * 3.14159 / 180), 2));
+        }
+        Zm(j, i) = sum ; 
+      }
+    }
+  }
+  return Zm ;
+}
+
+
+// distance matrix based on similarity of trajectory angles
+RcppExport SEXP distAngle(SEXP x, SEXP y, SEXP z){
+  NumericMatrix Lonm(x); // matrix of longitudes
+  NumericMatrix Latm(y); // matrix of latitudes
+  NumericMatrix Zm(z);
+  double A, B, C; // intermediate trajectory variables
+  double X0, Y0; // origin of lon, lat
+  double sum = 0.0;
+  int ncolumns = Lonm.ncol();
+  int nrows = Lonm.nrow();
+  X0 = (double) Lonm(nrows - 1, 0);
+  Y0 = (double) Latm(nrows - 1, 0);
+  for (int i = 0; i < ncolumns; i++) {
+    for (int j = 0; j < ncolumns; j++) {
+      sum = 0;
+      if (j > i) { // only need lower triangle
+        for (int k = 0; k < nrows - 1; k++) { // do not include origin
+          A = pow((Lonm(k, i) - X0), 2) + pow((Latm(k, i) - Y0), 2);
+          B = pow((Lonm(k, j) - X0), 2) + pow((Latm(k, j) - Y0), 2);
+          C = pow((Lonm(k, j) - Lonm(k, i)), 2) + pow((Latm(k, j) - Latm(k, i)), 2);
+          sum += acos((0.5 * (A + B - C) / pow(A * B, 0.5)));
+        }
+        Zm(j, i) = sum / (nrows - 1) ;
+      }
+    }
+  }
+  return Zm;
+}