Skip to content

Commit af8d503

Browse files
AlePaluAlePalu
committed
dcel cell support, dynamic edge insertion (also known as split_cell), shapefile reader allow to loop over single polygons of a multipolygon
1 parent 8cc95fb commit af8d503

File tree

4 files changed

+179
-127
lines changed

4 files changed

+179
-127
lines changed

‎fdaPDE/geoframe/shp.h

Lines changed: 93 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -14,8 +14,8 @@
1414
// You should have received a copy of the GNU General Public License
1515
// along with this program. If not, see <http://www.gnu.org/licenses/>.
1616

17-
#ifndef __SHAPEFILE_READER_H__
18-
#define __SHAPEFILE_READER_H__
17+
#ifndef __SHP_H__
18+
#define __SHP_H__
1919

2020
#include <array>
2121
#include <charconv>
@@ -25,7 +25,7 @@
2525
#include <string>
2626
#include <vector>
2727

28-
#include "../assert.h"
28+
#include "../utils/assert.h"
2929

3030
namespace fdapde {
3131

@@ -61,17 +61,19 @@ class shp_reader {
6161
int shape_type;
6262
};
6363
struct sf_point_t {
64+
int shape_type;
6465
int record_number;
6566
double x, y, z, m; // point coordinates (z and m optional)
6667

6768
sf_point_t() = default;
6869
sf_point_t(double x_, double y_) : x(x_), y(y_) { }
69-
sf_point_t(int record_number_, shp_reader& file) : record_number(record_number_) {
70+
sf_point_t(int record_number_, shp_reader& file) :
71+
shape_type(file.shape_type()), record_number(record_number_) {
7072
// point specific content
7173
x = read<double>(file.buff_, file.head_, LittleEndian);
7274
y = read<double>(file.buff_, file.head_, LittleEndian);
73-
if (file.shape_type() == ShapeType::PointM) { m = read<double>(file.buff_, file.head_, LittleEndian); }
74-
if (file.shape_type() == ShapeType::PointZ) {
75+
if (file.shape_type() == shape_t::PointM) { m = read<double>(file.buff_, file.head_, LittleEndian); }
76+
if (file.shape_type() == shape_t::PointZ) {
7577
z = read<double>(file.buff_, file.head_, LittleEndian);
7678
m = read<double>(file.buff_, file.head_, LittleEndian);
7779
}
@@ -86,14 +88,16 @@ class shp_reader {
8688
for (int i = 0; i < n_points; ++i) points[i] = read<double>(file.buff_, file.head_, LittleEndian);
8789
}
8890
public:
91+
int shape_type;
8992
int record_number;
9093
std::array<double, 4> bbox; // x_min, y_min, x_max, y_max
9194
std::array<double, 2> m_range, z_range; // m_min, m_max, z_min, z_max
9295
int n_points; // number of points in the record
9396
std::vector<double> x, y, z, m; // points coordinates (z and m optional)
9497

9598
sf_multipoint_t() = default;
96-
sf_multipoint_t(int record_number_, shp_reader& file) : record_number(record_number_) {
99+
sf_multipoint_t(int record_number_, shp_reader& file) :
100+
shape_type(file.shape_type()), record_number(record_number_) {
97101
// multipoint specific content
98102
for (int i = 0; i < 4; ++i) { bbox[i] = read<double>(file.buff_, file.head_, LittleEndian); }
99103
n_points = read<std::int32_t>(file.buff_, file.head_, LittleEndian);
@@ -103,8 +107,8 @@ class shp_reader {
103107
x[i] = read<double>(file.buff_, file.head_, LittleEndian);
104108
y[i] = read<double>(file.buff_, file.head_, LittleEndian);
105109
}
106-
if (file.shape_type() == ShapeType::MultiPointM) { read_zm_block(n_points, m_range, m, file); }
107-
if (file.shape_type() == ShapeType::MultiPointZ) {
110+
if (file.shape_type() == shape_t::MultiPointM) { read_zm_block(n_points, m_range, m, file); }
111+
if (file.shape_type() == shape_t::MultiPointZ) {
108112
read_zm_block(n_points, z_range, z, file);
109113
read_zm_block(n_points, m_range, m, file);
110114
}
@@ -117,15 +121,16 @@ class shp_reader {
117121
range[0] = read<double>(file.buff_, file.head_, LittleEndian);
118122
range[1] = read<double>(file.buff_, file.head_, LittleEndian);
119123
for (int i = 0; i < n_points; ++i) {
120-
if (file.shape_type() == ShapeType::PolygonZ || file.shape_type() == ShapeType::PolyLineZ) {
124+
if (file.shape_type() == shape_t::PolygonZ || file.shape_type() == shape_t::PolyLineZ) {
121125
points[i].z = read<double>(file.buff_, file.head_, LittleEndian);
122126
}
123-
if (file.shape_type() == ShapeType::PolygonM || file.shape_type() == ShapeType::PolyLineM) {
127+
if (file.shape_type() == shape_t::PolygonM || file.shape_type() == shape_t::PolyLineM) {
124128
points[i].m = read<double>(file.buff_, file.head_, LittleEndian);
125129
}
126130
}
127131
}
128132
public:
133+
int shape_type;
129134
int record_number;
130135
std::array<double, 4> bbox; // x_min, y_min, x_max, y_max
131136
std::array<double, 2> m_range, z_range; // m_min, m_max, z_min, z_max
@@ -135,7 +140,8 @@ class shp_reader {
135140
std::vector<sf_point_t> points; // points coordinates (z and m optional)
136141

137142
sf_polygon_t() = default;
138-
sf_polygon_t(int record_number_, shp_reader& file) : record_number(record_number_) {
143+
sf_polygon_t(int record_number_, shp_reader& file) :
144+
shape_type(file.shape_type()), record_number(record_number_) {
139145
// polygon specific content
140146
for (int i = 0; i < 4; ++i) { bbox[i] = read<double>(file.buff_, file.head_, LittleEndian); }
141147
n_rings = read<std::int32_t>(file.buff_, file.head_, LittleEndian);
@@ -159,10 +165,10 @@ class shp_reader {
159165
double y = read<double>(file.buff_, file.head_, LittleEndian);
160166
points.emplace_back(x, y);
161167
}
162-
if (file.shape_type() == ShapeType::PolygonM || file.shape_type() == ShapeType::PolyLineM) {
168+
if (file.shape_type() == shape_t::PolygonM || file.shape_type() == shape_t::PolyLineM) {
163169
read_zm_block(n_points, m_range, points, file);
164170
}
165-
if (file.shape_type() == ShapeType::PolygonZ || file.shape_type() == ShapeType::PolyLineZ) {
171+
if (file.shape_type() == shape_t::PolygonZ || file.shape_type() == shape_t::PolyLineZ) {
166172
read_zm_block(n_points, z_range, points, file);
167173
read_zm_block(n_points, m_range, points, file);
168174
}
@@ -187,6 +193,24 @@ class shp_reader {
187193
friend bool operator!=(const ring_iterator& it1, const ring_iterator& it2) {
188194
return it1.index_ != it2.index_;
189195
}
196+
int n_nodes() const { return polygon_->ring_end[index_] - polygon_->ring_begin[index_]; }
197+
// RowMajor expansion of polygon coordinates
198+
Eigen::Matrix<double, Dynamic, Dynamic> nodes() {
199+
int n_rows = n_nodes();
200+
int n_cols =
201+
2 +
202+
(polygon_->shape_type == shape_t::PointM ? 1 : (polygon_->shape_type == shape_t::PointZ ? 2 : 0));
203+
Eigen::Matrix<double, Dynamic, Dynamic> coords(n_rows, n_cols);
204+
int row = 0;
205+
for (auto jt = begin(), ht = end(); jt != ht; ++jt) {
206+
coords(row, 0) = jt->x;
207+
coords(row, 1) = jt->y;
208+
if (polygon_->shape_type == shape_t::PointM) { coords(row, 2) = jt->m; }
209+
if (polygon_->shape_type == shape_t::PointZ) { coords(row, 3) = jt->z; }
210+
row++;
211+
}
212+
return coords;
213+
}
190214
};
191215
ring_iterator begin() const { return ring_iterator(this, 0); }
192216
ring_iterator end() const { return ring_iterator(this, n_rings); }
@@ -219,7 +243,7 @@ class shp_reader {
219243
public:
220244
// supported shapefile format
221245
// all the non-null shapes in a shapefile are required to be of the same shape type (cit. Shapefile standard)
222-
enum ShapeType {
246+
enum shape_t {
223247
Point = 1, PolyLine = 3, Polygon = 5, MultiPoint = 8,
224248
PointZ = 11, PolyLineZ = 13, PolygonZ = 15, MultiPointZ = 18,
225249
PointM = 21, PolyLineM = 23, PolygonM = 25, MultiPointM = 28
@@ -245,13 +269,14 @@ class shp_reader {
245269
buff_ = new char[header_.file_length - header_.size];
246270
file.read(buff_, header_.file_length - header_.size);
247271
const int& st = header_.shape_type;
248-
if (st == ShapeType::Point || st == ShapeType::PointZ || st == ShapeType::PointM) read_into(points_);
249-
if (st == ShapeType::PolyLine || st == ShapeType::PolyLineZ || st == ShapeType::PolyLineM)
272+
if (st == shape_t::Point || st == shape_t::PointZ || st == shape_t::PointM) { read_into(points_); }
273+
if (st == shape_t::PolyLine || st == shape_t::PolyLineZ || st == shape_t::PolyLineM) {
250274
read_into(polylines_);
251-
if (st == ShapeType::Polygon || st == ShapeType::PolygonZ || st == ShapeType::PolygonM)
252-
read_into(polygons_);
253-
if (st == ShapeType::MultiPoint || st == ShapeType::MultiPointZ || st == ShapeType::MultiPointM)
275+
}
276+
if (st == shape_t::Polygon || st == shape_t::PolygonZ || st == shape_t::PolygonM) { read_into(polygons_); }
277+
if (st == shape_t::MultiPoint || st == shape_t::MultiPointZ || st == shape_t::MultiPointM) {
254278
read_into(multipoints_);
279+
}
255280
file.close();
256281
delete[] buff_;
257282
} else {
@@ -263,10 +288,32 @@ class shp_reader {
263288
std::array<double, 4> bbox() const { return {header_.bbox[0], header_.bbox[1], header_.bbox[2], header_.bbox[3]}; }
264289
int n_records() const { return n_records_; }
265290
const sf_header_t& header() const { return header_; }
266-
const std::vector<sf_point_t>& points() const { return points_; }
267-
const std::vector<sf_polyline_t>& polylines() const { return polylines_; }
268-
const std::vector<sf_polygon_t>& polygons() const { return polygons_; }
269-
const std::vector<sf_multipoint_t>& multipoints() const { return multipoints_; }
291+
int n_points() const { return points_.size(); }
292+
const sf_point_t& point(int index) const {
293+
fdapde_assert(
294+
shape_type() == shape_t::Point || shape_type() == shape_t::PointZ || shape_type() == shape_t::PointM);
295+
return points_[index];
296+
}
297+
int n_polylines() const { return polylines_.size(); }
298+
const sf_polyline_t& polyline(int index) const {
299+
fdapde_assert(
300+
shape_type() == shape_t::PolyLine || shape_type() == shape_t::PolyLineZ ||
301+
shape_type() == shape_t::PolyLineM);
302+
return polylines_[index];
303+
}
304+
int n_polygons() const { return polygons_.size(); }
305+
const sf_polygon_t& polygon(int index) const {
306+
fdapde_assert(
307+
shape_type() == shape_t::Polygon || shape_type() == shape_t::PolygonZ || shape_type() == shape_t::PolygonM);
308+
return polygons_[index];
309+
}
310+
int n_multipoints() const { return multipoints_.size(); }
311+
const sf_multipoint_t& multipoint(int index) const {
312+
fdapde_assert(
313+
shape_type() == shape_t::MultiPoint || shape_type() == shape_t::MultiPointZ ||
314+
shape_type() == shape_t::MultiPointM);
315+
return multipoints_[index];
316+
}
270317
};
271318

272319
// .dbf reader. file specification dBase level 5
@@ -356,7 +403,7 @@ class dbf_reader {
356403
return values;
357404
}
358405
}
359-
std::vector<std::pair<std::string, char>> field_descriptors() const { // vector of fields names and associated type
406+
std::vector<std::pair<std::string, char>> field_descriptors() const { // vector of (fields name, type) pairs
360407
std::vector<std::pair<std::string, char>> result;
361408
for (const auto& field : fields_) { result.emplace_back(field.name, field.type); }
362409
return result;
@@ -367,40 +414,36 @@ class ShapeFile {
367414
private:
368415
shp_reader shp_;
369416
dbf_reader dbf_;
370-
std::string gcs = "UNDEFINED"; // geographic coordinate system (GCS)
371-
std::string folder_name_;
417+
std::string gcs_ = "UNDEFINED"; // geographic coordinate system (GCS)
418+
std::string filename_;
372419
public:
373-
ShapeFile(std::string folder_name) : folder_name_(folder_name) {
374-
std::string file_name;
375-
std::size_t pos_start = 0, pos_end = 0;
376-
while ((pos_end = folder_name.find('/', pos_start)) != std::string::npos) { pos_start = pos_end + 1; }
377-
file_name = folder_name.substr(pos_start, pos_end - pos_start);
420+
ShapeFile(std::string filename) : filename_(filename) {
378421
// load geometric features and associated data
379-
shp_ = shp_reader(folder_name + "/" + file_name + ".shp");
380-
dbf_ = dbf_reader(folder_name + "/" + file_name + ".dbf");
422+
shp_ = shp_reader(filename_ + ".shp");
423+
dbf_ = dbf_reader(filename_ + ".dbf");
381424
// retrieve GCS informations from .prj file
382-
std::ifstream prj_file;
383-
prj_file.open(folder_name + "/" + file_name + ".prj");
384-
if (prj_file) {
425+
std::ifstream prj;
426+
prj.open(filename_ + ".prj");
427+
if (prj) {
385428
std::string line;
386-
getline(prj_file, line);
429+
getline(prj, line);
387430
std::size_t i = line.find("GEOGCS", 0);
388431
i += std::string("GEOGCS[\"").size();
389432
std::size_t j = line.find("\"", i);
390-
gcs = line.substr(i, j - i);
433+
gcs_ = line.substr(i, j - i);
391434
}
392435
}
393436
// getters
394437
const shp_reader& shp() const { return shp_; }
395438
template <typename T> std::vector<T> get_as(std::string colname) const { return dbf_.get_as<T>(colname); }
396439
std::vector<std::pair<std::string, char>> field_descriptors() const { return dbf_.field_descriptors(); }
397440
friend std::ostream& operator<<(std::ostream& os, const ShapeFile& sf) {
398-
os << "file: " << sf.folder_name_ << std::endl;
441+
os << "file: " << sf.filename_ << std::endl;
399442
std::string shape_type;
400-
if(sf.shp().shape_type() == shp_reader::ShapeType::Point) shape_type = "POINT";
401-
if(sf.shp().shape_type() == shp_reader::ShapeType::PolyLine) shape_type = "POLYLINE";
402-
if(sf.shp().shape_type() == shp_reader::ShapeType::Polygon) shape_type = "POLYGON";
403-
if(sf.shp().shape_type() == shp_reader::ShapeType::MultiPoint) shape_type = "MULTIPOINT";
443+
if(sf.shp().shape_type() == shp_reader::shape_t::Point) shape_type = "POINT";
444+
if(sf.shp().shape_type() == shp_reader::shape_t::PolyLine) shape_type = "POLYLINE";
445+
if(sf.shp().shape_type() == shp_reader::shape_t::Polygon) shape_type = "POLYGON";
446+
if(sf.shp().shape_type() == shp_reader::shape_t::MultiPoint) shape_type = "MULTIPOINT";
404447
os << "shape_type: " << shape_type << std::endl;
405448
os << "file size: " << sf.shp().header().file_length * 2 << " Bytes" << std::endl;
406449
os << "number of records: " << sf.shp().n_records() << std::endl;
@@ -415,6 +458,11 @@ class ShapeFile {
415458
}
416459
};
417460

461+
ShapeFile read_shp(const std::string& filename) {
462+
ShapeFile shp(filename);
463+
return shp;
464+
}
465+
418466
} // namespace fdapde
419467

420-
#endif // __SHAPEFILE_READER_H__
468+
#endif // __SHP_H__

0 commit comments

Comments
 (0)