distance computed the same way for msh and shp

works with unprojected meshes

dgpy/scripts/slimPre.py

@@ -868,7 +868,7 @@ def get_data_from_netcdf(nc_file_name, variable_name):
     val = f.variables[variable_name][:]
     return val
 
-def get_distance(region, physical_tags):
+def get_distance(region, physical_tags, dmax=1000):
     """Return the minimum distance from each node on the region to the physical tags. The index of the 1-D array returned is the nodal index, same as Region.coordinates 
     
     keyword arguments:
@@ -877,10 +877,12 @@ def get_distance(region, physical_tags):
         slimPre Region object of the nodes on which compute the distance
     * physical_tags
         vector containing the physical tags from which the minimum distance will be computed
+    * dmax : float (default: 1000)
+        points will be inserted on lines every dmax interval
     """
-    dist = distance_msh(region, 1000, physical_tags)
-    proj = Coordinate_system(region, data_proj="+proj=geocent +ellps=sphere +a=1 +b=1")
-    return dist.evaluate_matrix(proj.coordinates)
+    _regionTags = Region(region._mesh, physical_tags)
+    dist = distance_msh(_regionTags, dmax)
+    return dist.evaluate_region(region)
 
 def extrude(mesh_file_name, bath_file_name, nb_layers=None, z_layers=None, layers_function=None, mesh_file_name_out='', factor_show=0, addLayerDepthAtBot=1e10, periodicity=None):
 
@@ -1073,11 +1075,15 @@ class distance_msh(slim_private.pre_evaluator) :
         """
         if region._type is not "interface":
             dgpy.Msg.Fatal('distance_msh is only implemented for region of interfaces')
-        super(distance_msh, self).__init__(dgpy.slimPreEvaluatorDistance(region._groups, region._mesh._mesh_proj, dist))
+        if region._mesh._mesh_proj is None:
+            _mesh_proj = ""
+        else:
+            _mesh_proj = region._mesh._mesh_proj
+        super(distance_msh, self).__init__(dgpy.slimPreEvaluatorDistance(region._groups, _mesh_proj, dist))
 
 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 an ellipsoid."""
-    def __init__(self, shp_filename,  dmax, tags=[], ellps="+ellps=WGSG4") :
+    def __init__(self, shp_filename,  dmax, tags=[], ellps="+ellps=WGS84") :
         """
         arguments:
         * shp_filename : string

dgpy/scripts/slim_private.py

@@ -283,9 +283,12 @@ class pre_evaluator :
     def evaluate_region(self, region) :
         """ evaluate at all the points of a given region
         """
-        pjMesh = dgpy.proj.pj_init_plus(region._mesh._mesh_proj)
-        pjData = dgpy.proj.pj_init_plus(self.coordinate_system())
-        coordinates = dgpy.proj.pjTransform(pjMesh, pjData, region.coordinates)
-        dgpy.proj.pj_free(pjMesh)
-        dgpy.proj.pj_free(pjData)
-        return self.evaluate_matrix(coordinates)
+        if region._mesh._mesh_proj is None:
+            return self.evaluate_matrix(region.coordinates)
+        else:
+            pjMesh = dgpy.proj.pj_init_plus(region._mesh._mesh_proj)
+            pjData = dgpy.proj.pj_init_plus(self.coordinate_system())
+            coordinates = dgpy.proj.pjTransform(pjMesh, pjData, region.coordinates)
+            dgpy.proj.pj_free(pjMesh)
+            dgpy.proj.pj_free(pjData)
+            return self.evaluate_matrix(coordinates)

modules/slimFunction/slimPreEvaluatorDistance.cpp

@@ -10,6 +10,7 @@ namespace flann {namespace serialization {BASIC_TYPE_SERIALIZER(unsigned long lo
 #include "unref.h"
 slimPreEvaluatorDistance::slimPreEvaluatorDistance(std::string shp, double dmax, std::vector<int> tags, std::string ellipsoid) {
   _ellipsoid = ellipsoid;
+  _isProj = true; 
   auto vertices = unref::readSHP(shp.c_str(), dmax, ellipsoid);
   int n = 0;
   for (const auto &v : vertices){
@@ -30,11 +31,27 @@ slimPreEvaluatorDistance::slimPreEvaluatorDistance(std::string shp, double dmax,
 }
 
 slimPreEvaluatorDistance::slimPreEvaluatorDistance(std::vector<dgGroupOfFaces*> groups, std::string meshProj, double dmax) {
+  if (strcmp(meshProj.c_str(), "") == 0)
+    _isProj = false;
+  else
+    _isProj = true;
   fullMatrix<double> coord;
   std::vector<std::vector<double> > tmpcoord;
-
-  projPJ pjLocal = pj_init_plus(meshProj.c_str());
-  projPJ pjGlobal = pj_init_plus("+proj=geocent +ellps=sphere +a=1 +b=1");
+  double major_axis, eccentricity_squared;
+  projPJ pjLocal, pjGlobal;
+  if (_isProj){
+    pjLocal = pj_init_plus(meshProj.c_str());
+    pj_get_spheroid_defn(pjLocal, &major_axis, &eccentricity_squared);
+    std::string majAx = std::to_string(major_axis);
+    std::string eccSq = std::to_string(eccentricity_squared);
+    _ellipsoid = "+ellps=sphere +a=";
+    _ellipsoid.append(majAx);
+    _ellipsoid.append(" +es=");
+    _ellipsoid.append(eccSq);
+    std::string pjG = "+proj=geocent ";
+    pjG.append(_ellipsoid);
+    pjGlobal = pj_init_plus(pjG.c_str());
+  }
   int n = 0;
   for (size_t iFaceGroup = 0; iFaceGroup < groups.size(); ++iFaceGroup){
     dgGroupOfFaces *faceGroup= groups[iFaceGroup];
@@ -53,8 +70,10 @@ slimPreEvaluatorDistance::slimPreEvaluatorDistance(std::vector<dgGroupOfFaces*>
       double x1 = coord(cl[1], 0);
       double y1 = coord(cl[1], 1);
       double z1 = coord(cl[1], 2);
-      pj_transform(pjLocal,pjGlobal, 1, 1, &x0, &y0, &z0);
-      pj_transform(pjLocal,pjGlobal, 1, 1, &x1, &y1, &z1);
+      if (_isProj){
+        pj_transform(pjLocal,pjGlobal, 1, 1, &x0, &y0, &z0);
+        pj_transform(pjLocal,pjGlobal, 1, 1, &x1, &y1, &z1);
+      }
       double dist = sqrt((x0-x1)*(x0-x1)+(y0-y1)*(y0-y1)+(z0-z1)*(z0-z1));
       int nPoints = ceil(dist/dmax) + 2; 
       for (int i = 0; i < nPoints; ++i){
@@ -67,9 +86,10 @@ slimPreEvaluatorDistance::slimPreEvaluatorDistance(std::vector<dgGroupOfFaces*>
       }
     }
   }
-  pj_free(pjLocal);
-  pj_free(pjGlobal);
-
+  if (_isProj){
+    pj_free(pjLocal);
+    pj_free(pjGlobal);
+  }
   _points = new flann::Matrix<float>(new float[3*n],n,3);
   for (int i = 0; i < n; ++i)
     for (int j = 0; j < 3; ++j)
@@ -96,8 +116,12 @@ double slimPreEvaluatorDistance::evaluate(double x, double y, double z) {
   flann::Matrix<int> idM(&id,1,1);
   flann::Matrix<float> distM(&dist,1,1);
   _index->knnSearch(queryM, idM, distM, 1, flann::SearchParams(0));
-  double R = sqrt(x*x + y*y + z*z);
-  return 2*R*asin(sqrt(dist)/(2*R));
+  if (_isProj) {
+    double R = sqrt(x*x + y*y + z*z);
+    return 2*R*asin(sqrt(dist)/(2*R));
+  }
+  else
+    return sqrt(dist);
 }
 
 #endif

modules/slimFunction/slimPreEvaluatorDistance.h

@@ -15,6 +15,7 @@ class slimPreEvaluatorDistance:public slimPreEvaluator{
   flann::Index<flann::L2<float> > *_index;
   flann::Matrix<float> *_points;
   std::string _ellipsoid;
+  bool _isProj;
 #endif
  public:
   std::string coordinateSystem() {return "+proj=geocent " + _ellipsoid;}

modules/slimFunction/tests/distanceToTag/test.py

 from dgpy import *
 import numpy as np
-from scipy import spatial
-from dgpy.scripts.distance import distance
 import slimPre
 
-mesh = slimPre.Mesh("square.msh")
+mesh_file_name = "square.msh"
+mesh = slimPre.Mesh(mesh_file_name)
 groups = mesh._groups
-sol = distance(mesh, ["paste"], nPoints=2)
+region = slimPre.Region(mesh)
+dist = slimPre.get_distance(region, ["paste"])
+slimPre.write_file('dist.nc', region=region, time=None, data=[('dist', dist)])
+slimPre.netcdf_to_msh(mesh_file_name, 'dist.nc', 'dist', 'dist')
+sol = dgDofContainer(groups, 1)
+sol.importIdx('dist/dist.idx')
 
 x = [-1, -0.5, 0, 0.5, 1]
 y = [-1, 0, 1]