add slimPreEvaluatorGeoTIFF

benchmarks/bohai/mesh.py

@@ -12,6 +12,7 @@ dist = slimPre.distance_shp(coast_shp, 1000, [1,2,20])
 
 def fsize(x, y, z) :
     d = R*dist.evaluate_geocentric(x,y,z)
+    h = bath.evaluate_geocentric(x,y,z)
     if d < 1500 :
         s0 = 750
     elif d < 100000 :
@@ -19,7 +20,9 @@ def fsize(x, y, z) :
         s0 = 750*(1-f) + 20000 * f
     else :
         s0 = 20000
-    return s0/R
+    h = min(25,max(5, h))
+    s1 = 1000 * h**0.5
+    return min(s0, s1)/R
 
 physicals = {1:"coast", 2:"coast", 10:"yellow_sea_north", 11:"yellow_sea_center", 12:"yellow_sea_south", 20:"yellow_river"}
 

benchmarks/bohai/prepro.py

@@ -5,73 +5,44 @@ pre_data_dir_base = 'data_'+slimPre.partition_nb()+'/'
 pre_data_dir = pre_data_dir_base+slimPre.partition_id()+'/'
 slimPre.make_directory(pre_data_dir)
 
-mesh_file_name = slimPre.fetch_ftp('slim_data/bohai/bohai_v1.msh')
-
 initial_time = '2016-01-01 00:00:00'
 final_time = '2016-01-05 00:00:00'
 
-mesh_file_name = slimPre.partition_mesh(mesh_file_name, pre_data_dir+'mesh.msh')
+mesh_file_name = slimPre.partition_mesh('bohai.msh', pre_data_dir+"mesh.msh")
 mesh = slimPre.Mesh(mesh_file_name, mesh_proj='+proj=utm +ellps=WGS84 +zone=51')
 region_global = slimPre.Region(mesh, '')
-lonlat_global = slimPre.Coordinate_system(region_global, data_proj='+proj=latlong +ellps=WGS84')
 
 ### Bathymetry Process ###
-
 print('Preprocessing Bathymetry')
-from scipy.interpolate import griddata
-bath_data_file_name = slimPre.fetch_ftp('slim_data/bohai/bohai_bath.xyz')
-num_lines = sum(1 for line in open(bath_data_file_name)) -2
-cloud = np.empty((num_lines,3))
-fin = open(bath_data_file_name, 'r')
-l = fin.readline()
-l = fin.readline()
-i = 0
-while l != '' and l[0:9] != '<end_xyz>':
-    (x, y, b) = l.split(',')
-    cloud[i,:] = [float(x), float(y), float(b)]
-    l = fin.readline()
-    i = i+1
-
-lonlat_coordinates = lonlat_global.coordinates * 180./np.pi
-bath_linear = griddata(cloud[:,0:2], cloud[:,2], (lonlat_coordinates[:,0], lonlat_coordinates[:,1]), method='linear')
-bath_nearest = griddata(cloud[:,0:2], cloud[:,2], (lonlat_coordinates[:,0], lonlat_coordinates[:,1]), method='nearest')
-bath = np.where(np.isnan(bath_linear), bath_nearest, bath_linear)
-bath = np.maximum(bath, 2.)
-
+bath_tiff = slimPre.fetch_ftp('slim_data/bohai/bohai_bathymetry.tif')
+bath = slimPre.geo_tiff(bath_tiff).evaluate_region(region_global)
 slimPre.write_file(pre_data_dir+'bathymetry.nc', region=region_global, time=None, data=[('bathymetry', bath)])
-
 slimPre.smooth_bathymetry(mesh, (pre_data_dir+'bathymetry.nc','bathymetry'), output_file_name=pre_data_dir_base+'bathymetry_smooth', coefficient=0.5)
 
 ### Coriolis Process ###
-
 print('Preprocessing Coriolis')
-lonlat_coordinates = lonlat_global.coordinates
-corio = 2 * 7.292e-5 * np.sin(lonlat_coordinates[:,1])
+lonlat = slimPre.Coordinate_system(region_global, data_proj='+proj=latlong +ellps=WGS84').coordinates
+corio = 2 * 7.292e-5 * np.sin(lonlat[:,1])
 slimPre.write_file(pre_data_dir+'coriolis.nc', region=region_global, time=None, data=[('coriolis', corio)])
 
 ### River Discharge Process ###
-
 print('Preprocessing River Discharge')
 t = [3600*j for j in range(300)]
 time = slimPre.Time(time_vector=t, initial_time=initial_time)
-
 discharge = np.ones(len(t)) * 1000
 discharge[0] = 0
-
 slimPre.write_file(pre_data_dir+'river_discharge.nc', region=None, time=time, data=[('river_discharge', discharge)])
 
 ### Open Sea Process ###
-
 print('Preprocessing Open Sea Boundary')
 time = slimPre.Time(initial_time=initial_time, final_time=final_time, time_step=900.)
+
 # Yellow Sea North
 region = slimPre.Region(mesh, physical_tags=['yellow_sea_north'])
 (ssh, ssu, ssv) = slimPre.tpxo_tide(region, time)
 flow = slimPre.residual_flow(region, 0.1*1e6, pre_data_dir_base+'bathymetry_smooth/bathymetry_smooth.idx')
-
 ssu += flow[None,:,0]
 ssv += flow[None,:,1]
-
 slimPre.write_file(pre_data_dir+'yellow_sea_north.nc', region=region, time=time, data=[('h', ssh), ('u', ssu), ('v', ssv)])
 
 # Yellow Sea Center
@@ -79,17 +50,13 @@ region = slimPre.Region(mesh, physical_tags=['yellow_sea_center'])
 (ssh, ssu, ssv) = slimPre.tpxo_tide(region, time)
 slimPre.write_file(pre_data_dir+'yellow_sea_center.nc', region=region, time=time, data=[('h', ssh), ('u', ssu), ('v', ssv)])
 
-
 # Yellow Sea South
 region = slimPre.Region(mesh, physical_tags=['yellow_sea_south'])
 (ssh, ssu, ssv) = slimPre.tpxo_tide(region, time)
 flow = slimPre.residual_flow(region, -0.1*1e6, pre_data_dir_base+'bathymetry_smooth/bathymetry_smooth.idx')
-
 ssu += flow[None,:,0]
 ssv += flow[None,:,1]
-
 slimPre.write_file(pre_data_dir+'yellow_sea_south.nc', region=region, time=time, data=[('h', ssh), ('u', ssu), ('v', ssv)])
 
-
-print( 'Preprocessing done ' )
+print('Preprocessing done')
 slimPre.exit(0)

dgpy/scripts/slimPre.py

@@ -1009,3 +1009,13 @@ class distance_shp(slim_private.pre_evaluator) :
             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))
+
+class geo_tiff(slim_private.pre_evaluator) :
+    """ extract data from a geo TIFF file """
+    def __init__(self, tiff_filename) :
+        """
+        arguments:
+        * filename : string
+            path to the geo TIFF file
+        """
+        super(geo_tiff, self).__init__(dgpy.slimPreEvaluatorGeoTIFF(tiff_filename))

modules/slimFunction/CMakeLists.txt

@@ -13,6 +13,7 @@ set(SRC
   unrefMesher.cpp
   slimPreEvaluator.cpp
   slimPreEvaluatorDistanceSHP.cpp
+  slimPreEvaluatorGeoTIFF.cpp
 )
 
 set(LIBS "")
@@ -50,6 +51,7 @@ if(ENABLE_FES)
 endif(ENABLE_FES)
 
 option(ENABLE_PROJ "Enable PROJ" ON)
+option(ENABLE_GEOTIFF "Enable GEOTIFF (requires 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)
@@ -105,6 +107,25 @@ if(ENABLE_PROJ)
     set(HAVE_FLANN ON)
     list(APPEND EXTERNAL_MODULE_DEPS libproj)
   endif (ENABLE_FLANN)
+  if (ENABLE_GEOTIFF)
+    find_package(TIFF REQUIRED)
+    ExternalProject_Add(geotiff
+      URL http://download.osgeo.org/geotiff/libgeotiff/libgeotiff-1.4.2.tar.gz
+      URL_MD5 96ab80e0d4eff7820579957245d844f8
+      PREFIX libgeotiff
+      INSTALL_COMMAND ""
+      CMAKE_ARGS -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE} -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TOOLCHAIN_FILE} -DWITH_UTILITIES=NO -DCMAKE_C_FLAGS=-fPIC -DGEOTIFF_CSV_DATA_DIR=${CMAKE_CURRENT_BINARY_DIR}/libgeotiff/src/geotiff/csv
+      BINARY_DIR libgeotiff/build
+    )
+    include_directories(${CMAKE_CURRENT_BINARY_DIR}/libgeotiff/src/geotiff)
+    include_directories(${CMAKE_CURRENT_BINARY_DIR}/libgeotiff/src/geotiff/libxtiff)
+    include_directories(${CMAKE_CURRENT_BINARY_DIR}/libgeotiff/build)
+    list(APPEND LIBS ${CMAKE_CURRENT_BINARY_DIR}/libgeotiff/build/lib/libgeotiff.a ${CMAKE_CURRENT_BINARY_DIR}/libgeotiff/build/lib/libxtiff.a)
+    include_directories(${TIFF_INCLUDE_DIRS})
+    list(APPEND LIBS ${TIFF_LIBRARY})
+    list(APPEND EXTERNAL_MODULE_DEPS geotiff)
+    set(HAVE_GEOTIFF ON)
+  endif(ENABLE_GEOTIFF)
 endif(ENABLE_PROJ)
 
 configure_file(${CMAKE_CURRENT_SOURCE_DIR}/slimFunctionConfig.h.in ${CMAKE_CURRENT_BINARY_DIR}/slimFunctionConfig.h)

modules/slimFunction/slimFunction.i

@@ -15,6 +15,7 @@
   #include "unrefMesher.h"
   #include "slimPreEvaluator.h"
   #include "slimPreEvaluatorDistanceSHP.h"
+  #include "slimPreEvaluatorGeoTIFF.h"
 %}
 
 %import(module="dgpy.dgCommon") "dgCommon.i"
@@ -36,3 +37,4 @@
 %include "slimExportNetCDF.h"
 %include "slimPreEvaluator.h"
 %include "slimPreEvaluatorDistanceSHP.h"
+%include "slimPreEvaluatorGeoTIFF.h"

modules/slimFunction/slimFunctionConfig.h.in

@@ -7,4 +7,5 @@
 #cmakedefine HAVE_PROJ
 #cmakedefine HAVE_UNREF
 #cmakedefine HAVE_FLANN
+#cmakedefine HAVE_GEOTIFF
 #endif

modules/slimFunction/slimPreEvaluatorGeoTIFF.cpp

+#include "slimFunctionConfig.h"
+#ifdef HAVE_GEOTIFF
+#include "slimPreEvaluatorGeoTIFF.h"
+#include "xtiffio.h"
+#include "geotiff.h"
+#include "geovalues.h"
+#include "geo_normalize.h" // projection and conversion to proj4
+#include "tiffio.h"
+#include "proj_api.h"
+#include "dgMessage.h"
+#include "dgConfig.h"
+
+enum {SLIM_TIFF_UINT=1, SLIM_TIFF_INT=2, SLIM_TIFF_FLOAT=3};
+enum {SLIM_TIFF_BYTE=1, SLIM_TIFF_INT16, SLIM_TIFF_INT32, SLIM_TIFF_INT64, SLIM_TIFF_UINT16, SLIM_TIFF_UINT32, SLIM_TIFF_UINT64, SLIM_TIFF_FLOAT32, SLIM_TIFF_FLOAT64};
+
+static short buildType(short base, short size) {
+  if (size == 8) return SLIM_TIFF_BYTE;
+  if(base == SLIM_TIFF_INT && size == 16) return SLIM_TIFF_INT16;
+  if(base == SLIM_TIFF_INT && size == 32) return SLIM_TIFF_INT32;
+  if(base == SLIM_TIFF_INT && size == 64) return SLIM_TIFF_INT64;
+  if(base == SLIM_TIFF_UINT && size == 16) return SLIM_TIFF_UINT16;
+  if(base == SLIM_TIFF_UINT && size == 32) return SLIM_TIFF_UINT32;
+  if(base == SLIM_TIFF_UINT && size == 64) return SLIM_TIFF_UINT64;
+  if(base == SLIM_TIFF_FLOAT && size == 32) return SLIM_TIFF_FLOAT32;
+  if(base == SLIM_TIFF_FLOAT && size == 64) return SLIM_TIFF_FLOAT64;
+  Msg::Error("Unknown type in tiff file");
+  return 0;
+}
+
+static double _asDouble(char *data, short type, int i) {
+  switch (type) {
+    case SLIM_TIFF_BYTE: return ((uint8_t*)(data))[i];
+    case SLIM_TIFF_INT16: return ((int16_t*)(data))[i];
+    case SLIM_TIFF_INT32: return ((int32_t*)(data))[i];
+    case SLIM_TIFF_INT64: return ((int64_t*)(data))[i];
+    case SLIM_TIFF_UINT16: return ((uint16_t*)(data))[i];
+    case SLIM_TIFF_UINT32: return ((uint32_t*)(data))[i];
+    case SLIM_TIFF_UINT64: return ((uint64_t*)(data))[i];
+    case SLIM_TIFF_FLOAT32: return ((float*)(data))[i];
+    case SLIM_TIFF_FLOAT64: return ((double*)(data))[i];
+  }
+  return 0.;
+}
+
+slimPreEvaluatorGeoTIFF::slimPreEvaluatorGeoTIFF(const char *filename) {
+  _data = NULL;
+  _tiff = XTIFFOpen(filename, "r");
+  if(!_tiff)
+    Msg::Error("cannot open tiff file %s", filename);
+  _gtif = GTIFNew(_tiff);
+  if(!_gtif)
+    Msg::Error("cannot open tiff file %s", filename);
+  GTIFDefn	defn;
+  if(GTIFGetDefn(_gtif, &defn))
+    _coordinateSystem = GTIFGetProj4Defn(&defn);
+  else
+    Msg::Error("cannot read coordinate system from tiff file");
+  short int nsample, sampleformat, bitspersample;
+  if(TIFFGetField(_tiff, TIFFTAG_IMAGEWIDTH, &_sizex) != 1 ||
+    TIFFGetField(_tiff, TIFFTAG_IMAGELENGTH, &_sizey) != 1 ||
+    TIFFGetField(_tiff, TIFFTAG_SAMPLESPERPIXEL, &nsample) != 1 || 
+    TIFFGetField(_tiff, TIFFTAG_BITSPERSAMPLE, &bitspersample) != 1) {
+      Msg::Error("invalid tiff file");
+  }
+  if(TIFFGetField( _tiff, TIFFTAG_SAMPLEFORMAT, &sampleformat) != 1){
+    sampleformat = 1;
+  }
+  _dataType = buildType(sampleformat, bitspersample);
+
+  short int rasterType;
+  if(GTIFKeyGet( _gtif, GTRasterTypeGeoKey, &rasterType, 0, 1) != 1) 
+    rasterType = 1;
+  if(GTIFKeyGet( _gtif, GTModelTypeGeoKey, &_modelType, 0, 1) != 1) 
+    Msg::Error("no model type defined in geotiff file");
+  if(_modelType == ModelGeographic && GTIFKeyGet( _gtif, GeogAngularUnitsGeoKey, &_angularUnit, 0, 1) == 1) 
+    if(_angularUnit !=  Angular_Degree)
+      Msg::Error("geotiff unit is not decimal degrees (Angular_Degree)");
+      
+  printf("model type : %i\n", _modelType);
+  
+  /*
+  double *tiepoint;
+  int count;
+  TIFFGetField(_tiff, TIFFTAG_GEOTIEPOINTS, &count, &tiepoint);
+  TIFFGetField(_tiff, TIFFTAG_GEOPIXELSCALE, &count, &tiepoint);
+  */
+
+  _data = new char[bitspersample/8*_sizex*_sizey];
+  for (int i = 0; i < _sizey; ++i) {
+    if(TIFFReadScanline(_tiff, _data+i*_sizex*bitspersample/8, i, 0) !=1)
+      Msg::Error("error reading tiff file %s", filename);
+  }
+}
+
+std::string slimPreEvaluatorGeoTIFF::coordinateSystem() {
+  return _coordinateSystem;
+}
+
+double slimPreEvaluatorGeoTIFF::evaluate(double x, double y, double z) {
+  if(_modelType==2){
+    x*=180/M_PI;
+    y*=180/M_PI;
+  }
+  if(!GTIFPCSToImage(_gtif, &x, &y))
+      return 0;
+  int x0 = std::min(_sizex-1,std::max(0, (int)floor(x)));
+  int x1 = std::min(_sizex-1,std::max(0, (int)ceil(x)));
+  int y0 = std::min(_sizey-1,std::max(0, (int)floor(y)));
+  int y1 = std::min(_sizey-1,std::max(0, (int)ceil(y)));
+  double a = x-x0;
+  double b = y-y0;
+  return
+     (1-b)*(1-a)*(_asDouble(_data, _dataType, y0*_sizex+x0))
+    +(1-b)*(  a)*(_asDouble(_data, _dataType, y0*_sizex+x1))
+    +(  b)*(1-a)*(_asDouble(_data, _dataType, y1*_sizex+x0))
+    +(  b)*(  a)*(_asDouble(_data, _dataType, y1*_sizex+x1));
+}
+
+slimPreEvaluatorGeoTIFF::~slimPreEvaluatorGeoTIFF() {
+  if(_gtif)
+    GTIFFree(_gtif);
+  if (_tiff)
+    XTIFFClose(_tiff);
+  if (_data)
+    delete []_data;
+}
+
+#endif

modules/slimFunction/slimPreEvaluatorGeoTIFF.h

+#ifndef SLIM_PRE_GEOTIFF_H_
+#define SLIM_PRE_GEOTIFF_H_
+#include "slimFunctionConfig.h"
+#ifdef HAVE_GEOTIFF
+#include "slimPreEvaluator.h"
+#include "unrefMesher.h"
+#undef HAVE_STRINGS_H
+#undef HAVE_STRING_H
+#undef HAVE_STDLIB_H
+typedef struct tiff TIFF;
+typedef struct gtiff GTIF;
+
+class slimPreEvaluatorGeoTIFF : public slimPreEvaluator{
+  TIFF *_tiff;
+  GTIF *_gtif;
+  std::string _coordinateSystem;
+  short _modelType;
+  short _angularUnit;
+  short _dataType;
+  int _sizex, _sizey;
+  char *_data;
+ public:
+  slimPreEvaluatorGeoTIFF(const char *filename);
+  double evaluate(double x, double y, double z);
+  std::string coordinateSystem();
+  ~slimPreEvaluatorGeoTIFF();
+};
+#endif
+#endif