add slimPreEvaluatorDistanceSHP

benchmarks/bohai/mesh.py

@@ -11,7 +11,15 @@ bath = slimPre.geo_tiff(bath_tiff)
 dist = slimPre.distance_shp(coast_shp, 1000, [1,2,20])
 
 def fsize(x, y, z) :
-    return 10000/R
+    d = R*dist.evaluate_geocentric(x,y,z)
+    if d < 1500 :
+        s0 = 750
+    elif d < 100000 :
+        f = (d-1500)/(100000-1500)
+        s0 = 750*(1-f) + 20000 * f
+    else :
+        s0 = 20000
+    return s0/R
 
 physicals = {1:"coast", 2:"coast", 10:"yellow_sea_north", 11:"yellow_sea_center", 12:"yellow_sea_south", 20:"yellow_river"}
 

dgpy/scripts/slimPre.py

@@ -995,3 +995,17 @@ def mesh_shp(shp_inputs, lon, lat, mesh_size, mesh_size_min, output, projection
         to assign physicals names to the boundary elements of the mesh, a field (of type integer64) called "entity" should be present in the shape files. Then this parameter can be used to map the values of this field to physical names
     """
     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=[]) :
+        """
+        arguments:
+        * shp_filename : string
+            path to the shp file. The shape file can be composed of points, lines or polygons. The projection used in the shape file should be longitude-latitude (in degrees).
+        * dmax : float
+            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'
+        """
+        super(distance_shp, self).__init__(dgpy.slimPreEvaluatorDistanceSHP(shp_filename, dmax, tags))

dgpy/scripts/slim_private.py

@@ -259,3 +259,47 @@ class OpenBnd_3d:
         self.temperature = temperature
         self.flux = flux
         self.transport = transport
+
+class pre_evaluator :
+    """ base class for objects that can be evaluated one point at a time (in mesh generation)
+    or on a whole region (in preprocessing)
+    """
+    def __init__(self, cobject) :
+        """ constructor
+        arguments:
+        * cobject : dgpy.slimPreEvaluator
+        """
+        self._cobject = cobject
+    def coordinate_system(self) :
+        """ return the description of the projection in which the coordinates should be provided """
+        return self._cobject.coordinateSystem()
+    def evaluate_matrix(self, xyz) :
+        """ evaluate at several points
+        arguments :
+        * xyz : bi-dimensional numpy array
+           one point on each line (in the coordinate system specified by coordinate_system())
+        """
+        return self._cobject.evaluateMatrix(xyz)
+    def evaluate_geocentric(self, x, y, z) :
+        """ evaluate at a single point
+        arguments :
+        * x,y,z : float
+            the point in a cartesian coordinate system with its origin at the center of a sphere of radius 1
+        """
+        return self._cobject.evaluateGeocentric(x, y, z)
+    def evaluate(self, x, y, z) :
+        """ evaluate at a single point
+        arguments :
+        * x,y,z : float
+            the point in the projection specified by coordinate_system()
+        """
+        return self._cobject.evaluate(x, y, z)
+    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)

modules/slimFunction/CMakeLists.txt

@@ -11,6 +11,8 @@ set(SRC
   slimFES.cpp
   slimExportNetCDF.cpp
   unrefMesher.cpp
+  slimPreEvaluator.cpp
+  slimPreEvaluatorDistanceSHP.cpp
 )
 
 set(LIBS "")
@@ -49,6 +51,7 @@ endif(ENABLE_FES)
 
 option(ENABLE_PROJ "Enable PROJ" ON)
 option(ENABLE_UNREF "Enable UNREF mesh generator (requires ENABLE_PROJ)" ON)
+option(ENABLE_FLANN "Enable flann to compute distance from coastlines (requires ENABLE_PROJ and ENABLE_UNREF)" ON)
 if(ENABLE_PROJ)
   if(WIN32)
   set (PROJ_LIB ${CMAKE_CURRENT_BINARY_DIR}/libproj/build/lib/libproj_4_9.a)
@@ -89,6 +92,19 @@ if(ENABLE_PROJ)
     list(APPEND LIBS ${CMAKE_CURRENT_BINARY_DIR}/unref/build/libunref.a)
     include_directories(${CMAKE_CURRENT_BINARY_DIR}/unref/src/unref)
   endif(ENABLE_UNREF)
+  if (ENABLE_FLANN)
+    ExternalProject_Add(flann
+      URL http://www.cs.ubc.ca/research/flann/uploads/FLANN/flann-1.8.4-src.zip
+      URL_MD5 a0ecd46be2ee11a68d2a7d9c6b4ce701
+      PREFIX flann
+      CONFIGURE_COMMAND ""
+      BUILD_COMMAND ""
+      INSTALL_COMMAND ""
+    )
+    include_directories(${CMAKE_CURRENT_BINARY_DIR}/flann/src/flann/src/cpp)
+    set(HAVE_FLANN ON)
+    list(APPEND EXTERNAL_MODULE_DEPS libproj)
+  endif (ENABLE_FLANN)
 endif(ENABLE_PROJ)
 
 configure_file(${CMAKE_CURRENT_SOURCE_DIR}/slimFunctionConfig.h.in ${CMAKE_CURRENT_BINARY_DIR}/slimFunctionConfig.h)

modules/slimFunction/slimFunction.i

@@ -8,12 +8,13 @@
   #include "slimTpxo.h"
   #include "slimGebco.h"
   #include "slimTemporalSerie.h"
-  #include "slimFields.h"
   #include "slimLonLat.h"
   #include "slimIterateHelper.h"
   #include "slimFES.h"
   #include "slimExportNetCDF.h"
   #include "unrefMesher.h"
+  #include "slimPreEvaluator.h"
+  #include "slimPreEvaluatorDistanceSHP.h"
 %}
 
 %import(module="dgpy.dgCommon") "dgCommon.i"
@@ -33,3 +34,5 @@
 %include "slimIterateHelper.h"
 %include "slimFES.h"
 %include "slimExportNetCDF.h"
+%include "slimPreEvaluator.h"
+%include "slimPreEvaluatorDistanceSHP.h"

modules/slimFunction/slimFunctionConfig.h.in

@@ -6,4 +6,5 @@
 #cmakedefine FES_DATA
 #cmakedefine HAVE_PROJ
 #cmakedefine HAVE_UNREF
+#cmakedefine HAVE_FLANN
 #endif

modules/slimFunction/slimPreEvaluator.cpp

+#include "slimPreEvaluator.h"
+#include "proj_api.h"
+
+fullMatrix<double> slimPreEvaluator::evaluateMatrix(const fullMatrix<double> &xyz)
+{
+  fullMatrix<double> r(xyz.size1(), 1);
+  for (int i = 0; i < xyz.size1(); ++i)
+    r(i, 0) = evaluate(xyz(i,0), xyz(i,1), xyz(i,2));
+  return r;
+}
+
+slimPreEvaluator::slimPreEvaluator() {
+  pjLocal = pjGlobal = NULL;
+}
+
+slimPreEvaluator::~slimPreEvaluator() {
+  if(pjLocal)
+    pj_free(pjLocal);
+  if(pjGlobal)
+    pj_free(pjGlobal);
+}
+
+double slimPreEvaluator::evaluateGeocentric(double x, double y, double z){
+  if (pjLocal == NULL) {
+    pjLocal = pj_init_plus(coordinateSystem().c_str());
+    pjGlobal = pj_init_plus("+proj=geocent +ellps=sphere +a=1 +b=1");
+    if (!pjLocal || !pjGlobal) Msg::Fatal("invalid coordinate system : %s", coordinateSystem().c_str());
+  }
+  pj_transform(pjGlobal,pjLocal, 1, 1, &x, &y, &z);
+  return evaluate(x,y,z);
+}

modules/slimFunction/slimPreEvaluator.h

+#ifndef SLIM_PRE_EVALUATOR_H
+#define SLIM_PRE_EVALUATOR_H
+#include "slimFunctionConfig.h"
+#ifdef HAVE_PROJ
+#include "proj_api.h"
+#include "fullMatrix.h"
+// we follow the name convention of the python interface for this class
+class slimPreEvaluator { 
+  projPJ pjLocal, pjGlobal;
+  public :
+  virtual fullMatrix<double> evaluateMatrix(const fullMatrix<double> &xyz);
+  slimPreEvaluator();
+  virtual ~slimPreEvaluator();
+  virtual double evaluateGeocentric(double x, double y, double z);
+  virtual std::string coordinateSystem() = 0;
+  virtual double evaluate(double x, double y, double z) = 0;
+};
+#endif
+#endif

modules/slimFunction/slimPreEvaluatorDistanceSHP.cpp

+#include "slimFunctionConfig.h"
+#if defined(HAVE_UNREF) and defined(HAVE_FLANN)
+#include "slimPreEvaluatorDistanceSHP.h"
+#include <cmath> // required for flann
+#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);
+  int n = 0;
+  for (const auto &v : vertices){
+    if ((!tags.empty()) && std::find(tags.begin(), tags.end(), v->tag()) == tags.end()) continue;
+    n++;
+  }
+  _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();
+    (*_points)[n][1] = v->y();
+    (*_points)[n][2] = v->z();
+    n += 1;
+  }
+  _index = new flann::Index<flann::L2<float> >(*_points, flann::KDTreeSingleIndexParams());
+  _index->buildIndex();
+}
+
+slimPreEvaluatorDistanceSHP::~slimPreEvaluatorDistanceSHP() {
+  if (_index) {
+    delete _index;
+  }
+  if (_points) {
+    delete _points;
+    delete _points->ptr();
+  }
+}
+
+double slimPreEvaluatorDistanceSHP::evaluate(double x, double y, double z) {
+  float query[3] = {(float)x,(float)y,(float)z};
+  int id;
+  float dist;
+  flann::Matrix<float> queryM(query,1,3);
+  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));
+}
+
+#endif

modules/slimFunction/slimPreEvaluatorDistanceSHP.h

+#ifndef SLIM_PRE_EVLUATOR_DISTANCE_SHP_H_
+#define SLIM_PRE_EVLUATOR_DISTANCE_SHP_H_
+#include "slimFunctionConfig.h"
+#ifdef HAVE_PROJ
+#ifdef HAVE_UNREF
+#ifdef HAVE_FLANN
+#include "slimPreEvaluator.h"
+#ifndef SWIG
+#include <cmath> // required for flann
+#include "flann/flann.h"
+#endif
+class slimPreEvaluatorDistanceSHP:public slimPreEvaluator{
+#ifndef SWIG
+  flann::Index<flann::L2<float> > *_index;
+  flann::Matrix<float> *_points;
+#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();
+  double evaluate(double x, double y, double z) ;
+};
+#endif
+#endif
+#endif
+#endif