diff --git a/src/distVector.cpp b/src/distVector.cpp
index 1e016a210..05ed4fcce 100644
--- a/src/distVector.cpp
+++ b/src/distVector.cpp
@@ -19,11 +19,10 @@
 #include "spatVector.h"
 #include "distance.h"
 #include "geosphere.h"
+#include "vecmath.h"
 
 
-
-
-std::vector<double> SpatVector::nearestDistLonLat(std::vector<double>&x, std::vector<double>&y, std::string unit, std::string method) {
+std::vector<double> SpatVector::nearestDistLonLat(std::vector<double> x, std::vector<double> y, std::string unit, std::string method) {
 
 // for use with rasterize 
 	std::vector<double> d;
@@ -89,7 +88,7 @@ std::vector<double> SpatVector::nearestDistLonLat(std::vector<double>&x, std::ve
 					}
 					size_t nseg = vx.size() - 1;
 					for (size_t i=0; i<np; i++) {
-						if ((d[i] != 0) && (insect[i] == 0)) {
+						if (d[i] != 0) { // && (inside[i] == 0)) {
 							for (size_t j=0; j<nseg; j++) {
 								d[i] = std::min(d[i],
 									d2seg(x[i], y[i], vx[j], vy[j], vx[j+1], vy[j+1], r));
@@ -174,7 +173,11 @@ std::vector<double> SpatVector::distance(SpatVector x, bool pairwise, std::strin
 	std::string gtype = type();
 	std::string xtype = x.type();
 
-	if ((gtype != "points") || (xtype != "points")) {
+	if ((gtype == "points") && (xtype == "points")) {
+		std::vector<std::vector<double>> p = coordinates();
+		std::vector<std::vector<double>> px = x.coordinates();
+		return pointdistance(p[0], p[1], px[0], px[1], pairwise, m, lonlat, method);
+	} else if ((gtype == "points") || (xtype == "points")) {
 		if (lonlat) {
 			// not ok for multi-points
 			if (gtype == "points") {
@@ -185,19 +188,31 @@ std::vector<double> SpatVector::distance(SpatVector x, bool pairwise, std::strin
 				return nearestDistLonLat(xy[0], xy[1], unit, method);					
 			}
 		} else {
-			std::string distfun="";
-			d = geos_distance(x, pairwise, distfun);
-			if (m != 1) {
-				for (double &i : d) i *= m;
+			return geos_distance(x, pairwise, "", m);
+		}
+	} else {
+		if (lonlat) {
+			size_t n = size() * x.size();
+			d.reserve(n);
+
+			for (size_t i=0; i<size(); i++) {
+				SpatVector tmp1 = subset_rows({(int)i});
+				std::vector<std::vector<double>> xy1 = tmp1.coordinates();
+				for (size_t j=0; j<x.size(); j++) {
+					SpatVector tmp2 = x.subset_rows( {(int)j} );
+					std::vector<double> d1 = tmp2.nearestDistLonLat(xy1[0], xy1[1], unit, method);
+
+					std::vector<std::vector<double>> xy2 = tmp2.coordinates();
+					std::vector<double> d2 = tmp1.nearestDistLonLat(xy2[0], xy2[1], unit, method);
+					
+					d.push_back(std::min(vmin(d1, false), vmin(d2, false)));
+				}
 			}
-			return d;
+		} else {
+			d = geos_distance(x, pairwise, "", m);
 		}
 	}
-
-	std::vector<std::vector<double>> p = coordinates();
-	std::vector<std::vector<double>> px = x.coordinates();
-
-	return pointdistance(p[0], p[1], px[0], px[1], pairwise, m, lonlat, method);
+	return d;
 }
 
 
@@ -218,12 +233,7 @@ std::vector<double> SpatVector::distance(bool sequential, std::string unit, cons
 	}
 	std::string gtype = type();
 	if (gtype != "points") {
-		std::string distfun="";
-		d = geos_distance(sequential, distfun);
-		if (m != 1) {
-			for (double &i : d) i *= m;
-		}
-		return d;
+		return geos_distance(sequential, "", m);
 	} else {
 		if (sequential) {
 			std::vector<std::vector<double>> p = coordinates();
diff --git a/src/geos_methods.cpp b/src/geos_methods.cpp
index 0d1bb983a..0728609a2 100644
--- a/src/geos_methods.cpp
+++ b/src/geos_methods.cpp
@@ -2349,7 +2349,7 @@ bool get_dist_fun(dist_fn &f, std::string s) {
 }
 
 
-std::vector<double> SpatVector::geos_distance(SpatVector v, bool parallel, std::string fun) {
+std::vector<double> SpatVector::geos_distance(SpatVector v, bool parallel, std::string fun, double m) {
 
 	std::vector<double> out;
 
@@ -2415,10 +2415,13 @@ std::vector<double> SpatVector::geos_distance(SpatVector v, bool parallel, std::
 		}
 	}
 	geos_finish(hGEOSCtxt);
+	if (m != 1) {
+		for (double &d : out) d *= m;
+	}
 	return out;
  }
 
-std::vector<double> SpatVector::geos_distance(bool sequential, std::string fun) {
+std::vector<double> SpatVector::geos_distance(bool sequential, std::string fun, double m) {
 
 	std::vector<double> out;
 	dist_fn distfun;
@@ -2457,6 +2460,9 @@ std::vector<double> SpatVector::geos_distance(bool sequential, std::string fun)
 		out.push_back(0);
 	}
 	geos_finish(hGEOSCtxt);
+	if (m != 1) {
+		for (double &d : out) d *= m;
+	}
 	return out;
  }
 
@@ -2927,6 +2933,8 @@ SpatVector SpatVector::gaps() {
 
 
 
+// also use GEOSPreparedNearestPoints_r()
+
 SpatVector SpatVector::nearest_point(SpatVector v, bool parallel, const std::string method) {
 	SpatVector out;
 
diff --git a/src/spatVector.h b/src/spatVector.h
index c8038a3f3..7afcfa36f 100644
--- a/src/spatVector.h
+++ b/src/spatVector.h
@@ -203,8 +203,7 @@ class SpatVector {
 //		std::vector<double> pointdistance_seq(const std::vector<double>& px, const std::vector<double>& py, double m, bool lonlat);
 
 
-		std::vector<double> nearestDistLonLat(std::vector<double> &x, std::vector<double> &y, std::string unit, std::string method);
-
+		std::vector<double> nearestDistLonLat(std::vector<double> x, std::vector<double> y, std::string unit, std::string method);
 		std::vector<std::vector<size_t>> knearest(size_t k);
 
 		size_t size();
@@ -384,8 +383,8 @@ class SpatVector {
 		std::vector<unsigned> equals_exact(SpatVector v, double tol);
 		std::vector<unsigned> equals_exact(bool symmetrical, double tol);
 
-		std::vector<double> geos_distance(SpatVector v, bool parallel, std::string fun);
-		std::vector<double> geos_distance(bool sequential, std::string fun);
+		std::vector<double> geos_distance(SpatVector v, bool parallel, std::string fun, double m);
+		std::vector<double> geos_distance(bool sequential, std::string fun, double m);
 
 		SpatVector nearest_point(SpatVector v, bool parallel, const std::string method);
 		SpatVector nearest_point(const std::string method);