slimPreEvaluatorDistanceSHP : compute distance on an ellipsoid

benchmarks/bohai/mesh.py

@@ -11,7 +11,7 @@ bath = slimPre.geo_tiff(bath_tiff)
 dist = slimPre.distance_shp(coast_shp, 1000, [1,2,20])
 
 def fsize(x, y, z) :
-    d = R*dist.evaluate_geocentric(x,y,z)
+    d = dist.evaluate_geocentric(x,y,z)
     h = bath.evaluate_geocentric(x,y,z)
     if d < 1500 :
         s0 = 750

dgpy/scripts/slimPre.py

@@ -1013,8 +1013,8 @@ def mesh_shp(shp_inputs, lon, lat, mesh_size, mesh_size_min, output, projection
     dgpy.meshShp(shp_inputs, lon, lat, mesh_size, mesh_size_min, output, projection, boundary_names)
 
 class distance_shp(slim_private.pre_evaluator) :
-    """ compute the distance from a set of points given in an esri shape file (shp). Distances are evaluated on a sphere of radius 1."""
-    def __init__(self, shp_filename,  dmax, tags=[]) :
+    """ compute the distance from a set of points given in an esri shape file (shp). Distances are evaluated on an ellipsoid."""
+    def __init__(self, shp_filename,  dmax, tags=[], ellps="+ellps=WGSG4") :
         """
         arguments:
         * shp_filename : string
@@ -1023,6 +1023,8 @@ class distance_shp(slim_private.pre_evaluator) :
             if the shp file contains lines (or polygons), points will be inserted on those lines every dmax interval
         * tags : [int]
             if there is 'entity' field of integer associated with the shp features, this argument can be used to keep only the features with given values of 'entity'
+        * ellps : string (default: "+ellps=WGS84")
+            libproj string describing the ellipsoid
         """
         super(distance_shp, self).__init__(dgpy.slimPreEvaluatorDistanceSHP(shp_filename, dmax, tags))
 

modules/slimFunction/slimPreEvaluatorDistanceSHP.cpp

@@ -5,8 +5,9 @@
 #include "flann/flann.hpp"
 namespace flann {namespace serialization {BASIC_TYPE_SERIALIZER(unsigned long long);}} // hack to build build on mingw_w64 x86_64
 #include "unref.h"
-slimPreEvaluatorDistanceSHP::slimPreEvaluatorDistanceSHP(std::string shp, double dmax, std::vector<int> tags) {
-  auto vertices = unref::readSHP(shp.c_str(), dmax);
+slimPreEvaluatorDistanceSHP::slimPreEvaluatorDistanceSHP(std::string shp, double dmax, std::vector<int> tags, std::string ellipsoid) {
+  _ellipsoid = ellipsoid;
+  auto vertices = unref::readSHP(shp.c_str(), dmax, ellipsoid);
   int n = 0;
   for (const auto &v : vertices){
     if ((!tags.empty()) && std::find(tags.begin(), tags.end(), v->tag()) == tags.end()) continue;
@@ -14,6 +15,7 @@ slimPreEvaluatorDistanceSHP::slimPreEvaluatorDistanceSHP(std::string shp, double
   }
   _points = new flann::Matrix<float>(new float[3*n],n,3);
   n = 0;
+
   for (const auto &v : vertices) {
     if ((!tags.empty()) && std::find(tags.begin(), tags.end(), v->tag()) == tags.end()) continue;
     (*_points)[n][0] = v->x();

modules/slimFunction/slimPreEvaluatorDistanceSHP.h

@@ -13,10 +13,11 @@ class slimPreEvaluatorDistanceSHP:public slimPreEvaluator{
 #ifndef SWIG
   flann::Index<flann::L2<float> > *_index;
   flann::Matrix<float> *_points;
+  std::string _ellipsoid;
 #endif
  public:
-  slimPreEvaluatorDistanceSHP(std::string shp, double dmax, std::vector<int> tags = std::vector<int>());
-  std::string coordinateSystem() {return "+proj=geocent +ellps=sphere +a=1 +b=1";}
+  slimPreEvaluatorDistanceSHP(std::string shp, double dmax, std::vector<int> tags = std::vector<int>(), std::string ellipsoid="+ellps=WGS84");
+  std::string coordinateSystem() {return "+proj=geocent " + _ellipsoid;}
   ~slimPreEvaluatorDistanceSHP();
   double evaluate(double x, double y, double z) ;
 };

modules/slimFunction/unrefMesher.cpp

@@ -38,7 +38,16 @@ void meshShp(std::vector<std::string> shp_inputs, double lon, double lat, PyObje
   physicals[1] = boundary_names;
   std::vector<unref::TaggedVertex*> vertices;
   for (auto & filename : shp_inputs){
-    auto vs = unref::readSHP(filename.c_str(), lcmin);
+    auto vs = unref::readSHP(filename.c_str(), lcmin, "+ellps=sphere +a=1 +b=1");
+    for (auto v : vs) {
+      v->x() += rand() * 1.e-12/RAND_MAX;
+      v->y() += rand() * 1.e-12/RAND_MAX;
+      v->z() += rand() * 1.e-12/RAND_MAX;
+      double r = sqrt(v->x()*v->x() + v->y()*v->y() + v->z()*v->z());
+      v->x() /= r;
+      v->y() /= r;
+      v->z() /= r;
+    }
     vertices.insert(vertices.end(), vs.begin(), vs.end());
   }
   auto f = std::bind(&evalPythonCallback, callback, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3);