OGS
SHPInterface.cpp
Go to the documentation of this file.
1 
18 #include "SHPInterface.h"
19 
21 #include "BaseLib/Logging.h"
23 #include "GeoLib/GEOObjects.h"
24 #include "GeoLib/Point.h"
25 #include "GeoLib/Polyline.h"
27 #include "MeshLib/Elements/Quad.h"
28 #include "MeshLib/Elements/Tri.h"
29 #include "MeshLib/Mesh.h"
30 #include "MeshLib/Node.h"
31 
32 namespace FileIO
33 {
34 bool SHPInterface::readSHPInfo(const std::string& filename, int& shapeType,
35  int& numberOfEntities)
36 {
37  SHPHandle hSHP = SHPOpen(filename.c_str(), "rb");
38  if (!hSHP)
39  {
40  return false;
41  }
42 
43  double padfMinBound[4];
44  double padfMaxBound[4];
45 
46  // The SHPGetInfo() function retrieves various information about shapefile
47  // as a whole. The bounds are read from the file header, and may be
48  // inaccurate if the file was improperly generated.
49  SHPGetInfo(hSHP, &numberOfEntities, &shapeType, padfMinBound, padfMaxBound);
50 
51  SHPClose(hSHP);
52  return true;
53 }
54 
55 void SHPInterface::readSHPFile(const std::string& filename, OGSType choice,
56  const std::string& listName,
57  std::string const& gmsh_path)
58 {
59  int shapeType;
60  int numberOfElements;
61  double padfMinBound[4];
62  double padfMaxBound[4];
63 
64  SHPHandle hSHP = SHPOpen(filename.c_str(), "rb");
65  SHPGetInfo(hSHP, &numberOfElements, &shapeType, padfMinBound, padfMaxBound);
66 
67  if (((shapeType - 1) % 10 == 0) && (choice == SHPInterface::OGSType::POINT))
68  {
69  readPoints(hSHP, numberOfElements, listName);
70  }
71  if (((shapeType - 1) % 10 == 0) &&
73  {
74  readStations(hSHP, numberOfElements, listName);
75  }
76  if (((shapeType - 3) % 10 == 0 || (shapeType - 5) % 10 == 0) &&
78  {
79  readPolylines(hSHP, numberOfElements, listName);
80  }
81  if (((shapeType - 3) % 10 == 0 || (shapeType - 5) % 10 == 0) &&
83  {
84  readPolygons(hSHP, numberOfElements, listName, gmsh_path);
85  }
86 }
87 
88 void SHPInterface::readPoints(const SHPHandle& hSHP, int numberOfElements,
89  std::string listName)
90 {
91  if (numberOfElements > 0)
92  {
93  std::vector<GeoLib::Point*> points;
94  SHPObject* hSHPObject;
95 
96  for (int i = 0; i < numberOfElements; i++)
97  {
98  hSHPObject = SHPReadObject(hSHP, i);
99 
100  auto* pnt =
101  new GeoLib::Point(*(hSHPObject->padfX), *(hSHPObject->padfY),
102  *(hSHPObject->padfZ));
103  points.push_back(pnt);
104  }
105 
106  _geoObjects.addPointVec(std::move(points), listName,
108  SHPDestroyObject(hSHPObject); // de-allocate SHPObject
109  }
110 }
111 
112 void SHPInterface::readStations(const SHPHandle& hSHP, int numberOfElements,
113  std::string listName)
114 {
115  if (numberOfElements > 0)
116  {
117  std::vector<GeoLib::Point*> stations;
118  stations.reserve(numberOfElements);
119  SHPObject* hSHPObject;
120 
121  for (int i = 0; i < numberOfElements; i++)
122  {
123  hSHPObject = SHPReadObject(hSHP, i);
124  GeoLib::Station* stn =
125  GeoLib::Station::createStation(std::to_string(i),
126  *(hSHPObject->padfX),
127  *(hSHPObject->padfY),
128  *(hSHPObject->padfZ));
129  stations.push_back(stn);
130  }
131 
132  _geoObjects.addStationVec(std::move(stations), listName);
133  SHPDestroyObject(hSHPObject); // de-allocate SHPObject
134  }
135 }
136 
137 void SHPInterface::readPolylines(const SHPHandle& hSHP, int numberOfElements,
138  std::string listName)
139 {
140  if (numberOfElements <= 0)
141  {
142  return;
143  }
144  std::vector<GeoLib::Point*> pnts;
145  std::vector<GeoLib::Polyline*> lines;
146 
147  std::size_t pnt_id(0);
148  // for each polyline
149  for (int i = 0; i < numberOfElements; ++i)
150  {
151  SHPObject* hSHPObject = SHPReadObject(hSHP, i);
152  int const noOfPoints = hSHPObject->nVertices;
153  int const noOfParts = hSHPObject->nParts;
154 
155  for (int p = 0; p < noOfParts; ++p)
156  {
157  int const firstPnt = *(hSHPObject->panPartStart + p);
158  int const lastPnt = (p < (noOfParts - 1))
159  ? *(hSHPObject->panPartStart + p + 1)
160  : noOfPoints;
161 
162  // for each point in that polyline
163  for (int j = firstPnt; j < lastPnt; ++j)
164  {
165  pnts.push_back(new GeoLib::Point(
166  *(hSHPObject->padfX + j), *(hSHPObject->padfY + j),
167  *(hSHPObject->padfZ + j), pnt_id));
168  pnt_id++;
169  }
170  }
171  SHPDestroyObject(hSHPObject); // de-allocate SHPObject
172  }
173 
174  _geoObjects.addPointVec(std::move(pnts), listName,
176  GeoLib::PointVec const& points(*(_geoObjects.getPointVecObj(listName)));
177  std::vector<std::size_t> const& pnt_id_map(points.getIDMap());
178 
179  pnt_id = 0;
180  for (int i = 0; i < numberOfElements; ++i)
181  {
182  SHPObject* hSHPObject = SHPReadObject(hSHP, i);
183  int const noOfPoints = hSHPObject->nVertices;
184  int const noOfParts = hSHPObject->nParts;
185 
186  for (int p = 0; p < noOfParts; ++p)
187  {
188  // output for the first part of multipart polyline
189  if (noOfParts > 1 && p == 0)
190  {
191  INFO(
192  "Polygon {:d} consists of {:d} parts (PolylineIDs "
193  "{:d}-{:d}).",
194  i, noOfParts, lines.size(), lines.size() + noOfParts - 1);
195  }
196 
197  int const firstPnt = *(hSHPObject->panPartStart + p);
198  int const lastPnt = (p < (noOfParts - 1))
199  ? *(hSHPObject->panPartStart + p + 1)
200  : noOfPoints;
201 
202  auto* line = new GeoLib::Polyline(points.getVector());
203 
204  // create polyline
205  for (int j = firstPnt; j < lastPnt; ++j)
206  {
207  line->addPoint(pnt_id_map[pnt_id]);
208  pnt_id++;
209  }
210  // add polyline to polyline vector
211  lines.push_back(line);
212  }
213  SHPDestroyObject(hSHPObject); // de-allocate SHPObject
214  }
215  _geoObjects.addPolylineVec(std::move(lines), listName,
217 }
218 
219 void SHPInterface::readPolygons(const SHPHandle& hSHP, int numberOfElements,
220  const std::string& listName,
221  std::string const& gmsh_path)
222 {
223  readPolylines(hSHP, numberOfElements, listName);
224 
225  auto const polylines = _geoObjects.getPolylineVec(listName);
226 
227  for (auto const* polyline : *polylines)
228  {
229  INFO("Creating a surface by triangulation of the polyline ...");
230  if (FileIO::createSurface(*polyline, _geoObjects, listName, gmsh_path))
231  {
232  INFO("\t done");
233  }
234  else
235  {
236  WARN(
237  "\t Creating a surface by triangulation of the polyline "
238  "failed.");
239  }
240  }
241 }
242 
243 bool SHPInterface::write2dMeshToSHP(const std::string& file_name,
244  const MeshLib::Mesh& mesh)
245 {
246  if (mesh.getDimension() != 2)
247  {
248  ERR("SHPInterface::write2dMeshToSHP(): Mesh to Shape conversion is "
249  "only working for 2D Meshes.");
250  return false;
251  }
252 
253  std::size_t const n_elements = mesh.getNumberOfElements();
254  if (n_elements < 1)
255  {
256  ERR("SHPInterface::write2dMeshToSHP(): Mesh contains no elements.");
257  return false;
258  }
259 
260  // DBF-export requires a limit of records because of the limits to the
261  // integer values. The exponent controls maximum number of elements and
262  // the maximum number of digits written in the DBF file.
263  std::size_t const max_exp = 8;
264  if (n_elements >= std::pow(10, max_exp))
265  {
266  ERR("SHPInterface::write2dMeshToSHP(): Mesh contains too many elements "
267  "for currently implemented DBF-boundaries.");
268  return false;
269  }
270 
271  SHPHandle hSHP = SHPCreate(file_name.c_str(), SHPT_POLYGON);
272  DBFHandle hDBF = DBFCreate(file_name.c_str());
273  int const elem_id_field =
274  DBFAddField(hDBF, "Elem_ID", FTInteger, max_exp, 0);
275 
276  // Writing mesh elements to shape file
277  std::size_t shp_record(0);
278  std::vector<MeshLib::Element*> const& elems = mesh.getElements();
279  for (MeshLib::Element const* const e : elems)
280  {
281  // ignore all elements except triangles and quads
282  if ((e->getGeomType() == MeshLib::MeshElemType::TRIANGLE) ||
283  (e->getGeomType() == MeshLib::MeshElemType::QUAD))
284  {
285  SHPObject* object = createShapeObject(*e, shp_record);
286  SHPWriteObject(hSHP, -1, object);
287  SHPDestroyObject(object);
288 
289  // write element ID to DBF-file
290  DBFWriteIntegerAttribute(hDBF, shp_record, elem_id_field,
291  e->getID());
292  shp_record++;
293  }
294  }
295  SHPClose(hSHP);
296 
297  // Write scalar arrays to database file
298  int const n_recs = DBFGetRecordCount(hDBF);
299  for (auto [name, property] : mesh.getProperties())
300  {
301  if (auto p = dynamic_cast<MeshLib::PropertyVector<int>*>(property))
302  {
303  int const field = DBFAddField(hDBF, name.c_str(), FTInteger, 16, 0);
304  for (int i = 0; i < n_recs; ++i)
305  {
306  std::size_t const elem_idx =
307  DBFReadIntegerAttribute(hDBF, i, elem_id_field);
308  DBFWriteIntegerAttribute(hDBF, i, field, (*p)[elem_idx]);
309  }
310  }
311  else if (auto p =
312  dynamic_cast<MeshLib::PropertyVector<double>*>(property))
313  {
314  int const field = DBFAddField(hDBF, name.c_str(), FTDouble, 33, 16);
315  for (int i = 0; i < n_recs; ++i)
316  {
317  std::size_t const elem_idx =
318  DBFReadIntegerAttribute(hDBF, i, elem_id_field);
319  DBFWriteDoubleAttribute(hDBF, i, field, (*p)[elem_idx]);
320  }
321  }
322  }
323  DBFClose(hDBF);
324  INFO("Shape export of 2D mesh '{:s}' finished.", mesh.getName());
325  return true;
326 }
327 
329  std::size_t const shp_record)
330 {
331  unsigned const nNodes(e.getNumberOfBaseNodes());
332  double* padfX = new double[nNodes + 1];
333  double* padfY = new double[nNodes + 1];
334  double* padfZ = new double[nNodes + 1];
335  for (unsigned j = 0; j < nNodes; ++j)
336  {
337  padfX[j] = (*e.getNode(j))[0];
338  padfY[j] = (*e.getNode(j))[1];
339  padfZ[j] = (*e.getNode(j))[2];
340  }
341  // Last node == first node to close the polygon
342  padfX[nNodes] = (*e.getNode(0))[0];
343  padfY[nNodes] = (*e.getNode(0))[1];
344  padfZ[nNodes] = (*e.getNode(0))[2];
345 
346  // the generated shape object now handles the memory for padfX/Y/Z
347  return SHPCreateObject(SHPT_POLYGON, shp_record, 0, nullptr, nullptr,
348  nNodes + 1, padfX, padfY, padfZ, nullptr);
349 }
350 
351 } // namespace FileIO
Definition of analytical geometry functions.
Definition of the Element class.
Definition of the GEOObjects class.
Definition of the Point class.
void INFO(char const *fmt, Args const &... args)
Definition: Logging.h:32
void ERR(char const *fmt, Args const &... args)
Definition: Logging.h:42
void WARN(char const *fmt, Args const &... args)
Definition: Logging.h:37
Definition of the Mesh class.
Definition of the Node class.
Definition of the PolyLine class.
Definition of the Quad class.
Implementation of the SHPInterface class.
Definition of the Tri class.
void readPolylines(const SHPHandle &hSHP, int numberOfElements, std::string listName)
Reads lines into a vector of Polyline objects.
static bool readSHPInfo(const std::string &filename, int &shapeType, int &numberOfEntities)
Reads the header of the shape file.
void readPoints(const SHPHandle &hSHP, int numberOfElements, std::string listName)
Reads points into a vector of Point objects.
OGSType
Connection between ESRI type system for shape files and OGS GeoLib.
Definition: SHPInterface.h:49
void readStations(const SHPHandle &hSHP, int numberOfElements, std::string listName)
Reads points into a vector of Point objects and marks them as Station.
void readSHPFile(const std::string &filename, OGSType choice, const std::string &listName, std::string const &gmsh_path)
Reads data from the shape file.
GeoLib::GEOObjects & _geoObjects
Definition: SHPInterface.h:94
static bool write2dMeshToSHP(const std::string &file_name, const MeshLib::Mesh &mesh)
static SHPObject * createShapeObject(MeshLib::Element const &e, std::size_t const shp_record)
Creates a shape object polygon out of a 2D mesh element.
void readPolygons(const SHPHandle &hSHP, int numberOfElements, const std::string &listName, std::string const &gmsh_path)
Reads lines into a vector of Polyline and Surface objects.
void addPolylineVec(std::vector< Polyline * > &&lines, std::string const &name, PolylineVec::NameIdMap &&ply_names)
Definition: GEOObjects.cpp:152
void addPointVec(std::vector< Point * > &&points, std::string &name, PointVec::NameIdMap &&pnt_id_name_map, double const eps=std::sqrt(std::numeric_limits< double >::epsilon()))
Definition: GEOObjects.cpp:47
const PointVec * getPointVecObj(const std::string &name) const
Definition: GEOObjects.cpp:85
void addStationVec(std::vector< Point * > &&stations, std::string &name)
Adds a vector of stations with the given name and colour to GEOObjects.
Definition: GEOObjects.cpp:123
const std::vector< Polyline * > * getPolylineVec(const std::string &name) const
Definition: GEOObjects.cpp:189
This class manages pointers to Points in a std::vector along with a name. It also handles the deletio...
Definition: PointVec.h:38
const std::vector< std::size_t > & getIDMap() const
Definition: PointVec.h:96
Class Polyline consists mainly of a reference to a point vector and a vector that stores the indices ...
Definition: Polyline.h:53
A Station (observation site) is basically a Point with some additional information.
Definition: Station.h:37
static Station * createStation(const std::string &line)
Creates a Station-object from information contained in a string (assuming the string has the right fo...
Definition: Station.cpp:45
std::map< std::string, std::size_t > NameIdMap
Definition: TemplateVec.h:44
std::vector< T * > const & getVector() const
Definition: TemplateVec.h:109
virtual const Node * getNode(unsigned idx) const =0
virtual unsigned getNumberOfBaseNodes() const =0
unsigned getDimension() const
Returns the dimension of the mesh (determined by the maximum dimension over all elements).
Definition: Mesh.h:71
std::vector< Element * > const & getElements() const
Get the element-vector for the mesh.
Definition: Mesh.h:98
const std::string getName() const
Get name of the mesh.
Definition: Mesh.h:92
Properties & getProperties()
Definition: Mesh.h:123
std::size_t getNumberOfElements() const
Get the number of elements.
Definition: Mesh.h:86
bool createSurface(GeoLib::Polyline const &ply, GeoLib::GEOObjects &geometries, std::string const &geometry_name, std::string const &gmsh_binary)
TemplateElement< PointRule1 > Point
Definition: Point.h:20