Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: netCDF read/write Driver
4 : * Purpose: GDAL bindings over netCDF library.
5 : * Author: Winor Chen <wchen329 at wisc.edu>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2019, Winor Chen <wchen329 at wisc.edu>
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 : #include <cstdio>
13 : #include <cstring>
14 : #include <set>
15 : #include <vector>
16 : #include "netcdf.h"
17 : #include "netcdfdataset.h"
18 : #include "netcdfsg.h"
19 :
20 : namespace nccfdriver
21 : {
22 :
23 : /* SGeometry_Reader
24 : * (implementations)
25 : *
26 : */
27 73 : SGeometry_Reader::SGeometry_Reader(int ncId, int geoVarId)
28 139 : : gc_varId(geoVarId), touple_order(0)
29 : {
30 :
31 : char container_name[NC_MAX_NAME + 1];
32 73 : memset(container_name, 0, NC_MAX_NAME + 1);
33 :
34 : // Get geometry container name
35 73 : if (nc_inq_varname(ncId, geoVarId, container_name) != NC_NOERR)
36 : {
37 : throw SG_Exception_Existential("new geometry container",
38 0 : "the variable of the given ID");
39 : }
40 :
41 : // Establish string version of container_name
42 73 : container_name_s = std::string(container_name);
43 :
44 : // Find geometry type
45 73 : this->type = nccfdriver::getGeometryType(ncId, geoVarId);
46 :
47 73 : if (this->type == NONE)
48 : {
49 : throw SG_Exception_Existential(
50 1 : static_cast<const char *>(container_name), CF_SG_GEOMETRY_TYPE);
51 : }
52 :
53 : // Get grid mapping variable, if it exists
54 72 : this->gm_varId = INVALID_VAR_ID;
55 118 : if (NCDFGetAttr(ncId, geoVarId, CF_GRD_MAPPING, gm_name_s) == CE_None &&
56 46 : !gm_name_s.empty())
57 : {
58 46 : const char *gm_name = gm_name_s.c_str();
59 : int gmVID;
60 46 : if (nc_inq_varid(ncId, gm_name, &gmVID) == NC_NOERR)
61 : {
62 46 : this->gm_varId = gmVID;
63 : }
64 : }
65 :
66 : // Find a list of node counts and part node count
67 144 : std::string nc_name_s;
68 144 : std::string pnc_name_s;
69 144 : std::string ir_name_s;
70 72 : int pnc_vid = INVALID_VAR_ID;
71 72 : int nc_vid = INVALID_VAR_ID;
72 72 : int ir_vid = INVALID_VAR_ID;
73 : int buf;
74 72 : size_t bound = 0;
75 72 : size_t total_node_count = 0; // used in error checks later
76 72 : size_t total_part_node_count = 0;
77 :
78 130 : if (NCDFGetAttr(ncId, geoVarId, CF_SG_NODE_COUNT, nc_name_s) == CE_None &&
79 58 : !nc_name_s.empty())
80 : {
81 58 : const char *nc_name = nc_name_s.c_str();
82 58 : nc_inq_varid(ncId, nc_name, &nc_vid);
83 651 : while (nc_get_var1_int(ncId, nc_vid, &bound, &buf) == NC_NOERR)
84 : {
85 594 : if (buf < 0)
86 : {
87 : throw SG_Exception_Value_Violation(
88 : static_cast<const char *>(container_name), CF_SG_NODE_COUNT,
89 1 : "negative");
90 : }
91 :
92 593 : this->node_counts.push_back(buf);
93 593 : total_node_count += buf;
94 593 : bound++;
95 : }
96 : }
97 :
98 71 : if (NCDFGetAttr(ncId, geoVarId, CF_SG_PART_NODE_COUNT, pnc_name_s) ==
99 109 : CE_None &&
100 38 : !pnc_name_s.empty())
101 : {
102 38 : const char *pnc_name = pnc_name_s.c_str();
103 38 : bound = 0;
104 38 : nc_inq_varid(ncId, pnc_name, &pnc_vid);
105 870 : while (nc_get_var1_int(ncId, pnc_vid, &bound, &buf) == NC_NOERR)
106 : {
107 833 : if (buf < 0)
108 : {
109 : throw SG_Exception_Value_Violation(
110 : static_cast<const char *>(container_name),
111 1 : CF_SG_PART_NODE_COUNT, "negative");
112 : }
113 :
114 832 : this->pnode_counts.push_back(buf);
115 832 : total_part_node_count += buf;
116 832 : bound++;
117 : }
118 : }
119 :
120 70 : if (NCDFGetAttr(ncId, geoVarId, CF_SG_INTERIOR_RING, ir_name_s) ==
121 82 : CE_None &&
122 12 : !ir_name_s.empty())
123 : {
124 12 : const char *ir_name = ir_name_s.c_str();
125 12 : bound = 0;
126 12 : nc_inq_varid(ncId, ir_name, &ir_vid);
127 295 : while (nc_get_var1_int(ncId, ir_vid, &bound, &buf) == NC_NOERR)
128 : {
129 284 : bool store = buf == 0 ? false : true;
130 :
131 284 : if (buf != 0 && buf != 1)
132 : {
133 : throw SG_Exception_Value_Required(
134 : static_cast<const char *>(container_name),
135 1 : CF_SG_INTERIOR_RING, "0 or 1");
136 : }
137 :
138 283 : this->int_rings.push_back(store);
139 283 : bound++;
140 : }
141 : }
142 :
143 : /* Enforcement of well formed CF files
144 : * If these are not met then the dataset is malformed and will halt further
145 : * processing of simple geometries.
146 : */
147 :
148 : // part node count exists only when node count exists
149 69 : if (pnode_counts.size() > 0 && node_counts.size() == 0)
150 : {
151 : throw SG_Exception_Dep(static_cast<const char *>(container_name),
152 0 : CF_SG_PART_NODE_COUNT, CF_SG_NODE_COUNT);
153 : }
154 :
155 : // part node counts and node counts don't match up in count
156 69 : if (pnc_vid != INVALID_VAR_ID && total_node_count != total_part_node_count)
157 : {
158 : throw SG_Exception_BadSum(static_cast<const char *>(container_name),
159 0 : CF_SG_PART_NODE_COUNT, CF_SG_PART_NODE_COUNT);
160 : }
161 :
162 : // interior rings only exist when part node counts exist
163 69 : if (int_rings.size() > 0 && pnode_counts.size() == 0)
164 : {
165 : throw SG_Exception_Dep(static_cast<const char *>(container_name),
166 0 : CF_SG_INTERIOR_RING, CF_SG_PART_NODE_COUNT);
167 : }
168 :
169 : // cardinality of part_node_counts == cardinality of interior_ring (if
170 : // interior ring > 0)
171 69 : if (int_rings.size() > 0)
172 : {
173 11 : if (int_rings.size() != pnode_counts.size())
174 : {
175 : throw SG_Exception_Dim_MM(static_cast<const char *>(container_name),
176 : CF_SG_INTERIOR_RING,
177 0 : CF_SG_PART_NODE_COUNT);
178 : }
179 : }
180 :
181 : // lines and polygons require node_counts, multipolygons are checked with
182 : // part_node_count
183 69 : if (this->type == POLYGON || this->type == LINE)
184 : {
185 14 : if (node_counts.size() < 1)
186 : {
187 : throw SG_Exception_Existential(
188 1 : static_cast<const char *>(container_name), CF_SG_NODE_COUNT);
189 : }
190 : }
191 :
192 : /* Basic Safety checks End
193 : */
194 :
195 : // Create bound list
196 68 : size_t rc = 0;
197 68 : bound_list.push_back(0); // start with 0
198 :
199 68 : if (node_counts.size() > 0)
200 : {
201 587 : for (size_t i = 0; i < node_counts.size() - 1; i++)
202 : {
203 532 : rc = rc + node_counts[i];
204 532 : bound_list.push_back(rc);
205 : }
206 : }
207 :
208 69 : std::string cart_s;
209 : // Node Coordinates
210 68 : if (NCDFGetAttr(ncId, geoVarId, CF_SG_NODE_COORDINATES, cart_s) !=
211 136 : CE_None ||
212 68 : cart_s.empty())
213 : {
214 0 : throw SG_Exception_Existential(container_name, CF_SG_NODE_COORDINATES);
215 : }
216 :
217 : // Create parts count list and an offset list for parts indexing
218 68 : if (this->node_counts.size() > 0)
219 : {
220 55 : int ind = 0;
221 55 : int parts = 0;
222 55 : int prog = 0;
223 55 : int c = 0;
224 :
225 883 : for (size_t pcnt = 0; pcnt < pnode_counts.size(); pcnt++)
226 : {
227 828 : if (prog == 0)
228 521 : pnc_bl.push_back(pcnt);
229 :
230 828 : if (int_rings.size() > 0 && !int_rings[pcnt])
231 236 : c++;
232 :
233 828 : prog = prog + pnode_counts[pcnt];
234 828 : parts++;
235 :
236 828 : if (prog == node_counts[ind])
237 : {
238 521 : ind++;
239 521 : this->parts_count.push_back(parts);
240 521 : if (int_rings.size() > 0)
241 230 : this->poly_count.push_back(c);
242 521 : c = 0;
243 521 : prog = 0;
244 521 : parts = 0;
245 : }
246 307 : else if (prog > node_counts[ind])
247 : {
248 : throw SG_Exception_BadSum(container_name, CF_SG_PART_NODE_COUNT,
249 0 : CF_SG_NODE_COUNT);
250 : }
251 : }
252 : }
253 :
254 : // (1) the tuple order for a single point
255 : // (2) the variable ids with the relevant coordinate values
256 68 : int X = INVALID_VAR_ID;
257 68 : int Y = INVALID_VAR_ID;
258 68 : int Z = INVALID_VAR_ID;
259 :
260 : char cart[NC_MAX_NAME + 1];
261 68 : memset(cart, 0, NC_MAX_NAME + 1);
262 68 : strncpy(cart, cart_s.c_str(), NC_MAX_NAME);
263 :
264 68 : char *dim = strtok(cart, " ");
265 68 : int axis_id = 0;
266 :
267 228 : while (dim != nullptr)
268 : {
269 160 : if (nc_inq_varid(ncId, dim, &axis_id) == NC_NOERR)
270 : {
271 :
272 : // Check axis signature
273 160 : std::string a_sig;
274 160 : NCDFGetAttr(ncId, axis_id, CF_AXIS, a_sig);
275 :
276 : // If valid signify axis correctly
277 160 : if (a_sig == "X")
278 : {
279 68 : X = axis_id;
280 : }
281 92 : else if (a_sig == "Y")
282 : {
283 68 : Y = axis_id;
284 : }
285 24 : else if (a_sig == "Z")
286 : {
287 24 : Z = axis_id;
288 : }
289 : else
290 : {
291 : throw SG_Exception_Dep(container_name,
292 0 : "A node_coordinates variable", CF_AXIS);
293 : }
294 :
295 160 : this->touple_order++;
296 : }
297 : else
298 : {
299 0 : throw SG_Exception_Existential(container_name, dim);
300 : }
301 :
302 160 : dim = strtok(nullptr, " ");
303 : }
304 :
305 : // Write axis in X, Y, Z order
306 :
307 68 : if (X != INVALID_VAR_ID)
308 68 : this->nodec_varIds.push_back(X);
309 : else
310 : {
311 : throw SG_Exception_Existential(container_name,
312 0 : "node_coordinates: X-axis");
313 : }
314 68 : if (Y != INVALID_VAR_ID)
315 68 : this->nodec_varIds.push_back(Y);
316 : else
317 : {
318 : throw SG_Exception_Existential(container_name,
319 0 : "node_coordinates: Y-axis");
320 : }
321 68 : if (Z != INVALID_VAR_ID)
322 24 : this->nodec_varIds.push_back(Z);
323 :
324 : /* Final Checks for node coordinates
325 : * (1) Each axis has one and only one dimension, dim length of each node
326 : * coordinates are all equal (2) total node count == dim length of each node
327 : * coordinates (if node_counts not empty) (3) there are at least two node
328 : * coordinate variable ids
329 : */
330 :
331 68 : int all_dim = INVALID_VAR_ID;
332 68 : bool dim_set = false;
333 : //(1) one dimension check, each node_coordinates have same dimension
334 227 : for (size_t nvitr = 0; nvitr < nodec_varIds.size(); nvitr++)
335 : {
336 160 : int dimC = 0;
337 160 : nc_inq_varndims(ncId, nodec_varIds[nvitr], &dimC);
338 :
339 160 : if (dimC != 1)
340 : {
341 0 : throw SG_Exception_Not1D();
342 : }
343 :
344 : // check that all have same dimension
345 : int inter_dim[1];
346 160 : if (nc_inq_vardimid(ncId, nodec_varIds[nvitr], inter_dim) != NC_NOERR)
347 : {
348 : throw SG_Exception_Existential(
349 0 : container_name, "one or more node_coordinate dimensions");
350 : }
351 :
352 160 : if (!dim_set)
353 : {
354 68 : all_dim = inter_dim[0];
355 68 : dim_set = true;
356 : }
357 :
358 : else
359 : {
360 92 : if (inter_dim[0] != all_dim)
361 : throw SG_Exception_Dim_MM(
362 : container_name, "X, Y",
363 1 : "in general all node coordinate axes");
364 : }
365 : }
366 :
367 : // (2) check equality one
368 67 : if (node_counts.size() > 0)
369 : {
370 54 : size_t diml = 0;
371 54 : nc_inq_dimlen(ncId, all_dim, &diml);
372 :
373 54 : if (diml != total_node_count)
374 : throw SG_Exception_BadSum(container_name, "node_count",
375 0 : "node coordinate dimension length");
376 : }
377 :
378 : // (3) check tuple order
379 67 : if (this->touple_order < 2)
380 : {
381 : throw SG_Exception_Existential(
382 : container_name,
383 0 : "insufficient node coordinates must have at least two axis");
384 : }
385 :
386 : /* Investigate for instance dimension
387 : * The procedure is as follows
388 : *
389 : * (1) if there's node_count, use the dimension used to index node count
390 : * (2) otherwise it's point (singleton) data, in this case use the node
391 : * coordinate dimension
392 : */
393 67 : size_t instance_dim_len = 0;
394 :
395 67 : if (node_counts.size() >= 1)
396 : {
397 54 : int nc_dims = 0;
398 54 : nc_inq_varndims(ncId, nc_vid, &nc_dims);
399 :
400 54 : if (nc_dims != 1)
401 0 : throw SG_Exception_Not1D();
402 :
403 : int nc_dim_id[1];
404 :
405 54 : if (nc_inq_vardimid(ncId, nc_vid, nc_dim_id) != NC_NOERR)
406 : {
407 : throw SG_Exception_Existential(container_name,
408 0 : "node_count dimension");
409 : }
410 :
411 54 : this->inst_dimId = nc_dim_id[0];
412 : }
413 :
414 : else
415 : {
416 13 : this->inst_dimId = all_dim;
417 : }
418 :
419 67 : nc_inq_dimlen(ncId, this->inst_dimId, &instance_dim_len);
420 :
421 67 : if (instance_dim_len == 0)
422 0 : throw SG_Exception_EmptyDim();
423 :
424 : // Set values accordingly
425 67 : this->inst_dimLen = instance_dim_len;
426 67 : this->pt_buffer = std::make_unique<Point>(this->touple_order);
427 67 : this->gc_varId = geoVarId;
428 67 : this->ncid = ncId;
429 67 : }
430 :
431 67341 : Point &SGeometry_Reader::operator[](size_t index)
432 : {
433 202262 : for (int order = 0; order < touple_order; order++)
434 : {
435 134921 : Point &pt = *(this->pt_buffer);
436 : double data;
437 134921 : size_t real_ind = index;
438 :
439 : // Read a single coord
440 : int err =
441 134921 : nc_get_var1_double(ncid, nodec_varIds[order], &real_ind, &data);
442 :
443 134921 : if (err != NC_NOERR)
444 : {
445 0 : throw SG_Exception_BadPoint();
446 : }
447 :
448 134921 : pt[order] = data;
449 : }
450 :
451 67341 : return *(this->pt_buffer);
452 : }
453 :
454 529 : size_t SGeometry_Reader::get_geometry_count()
455 : {
456 529 : if (type == POINT)
457 : {
458 47 : if (this->nodec_varIds.size() < 1)
459 0 : return 0;
460 :
461 : // If more than one dim, then error. Otherwise inquire its length and
462 : // return that
463 : int dims;
464 47 : if (nc_inq_varndims(this->ncid, nodec_varIds[0], &dims) != NC_NOERR)
465 0 : return 0;
466 47 : if (dims != 1)
467 0 : return 0;
468 :
469 : // Find which dimension is used for x
470 : int index;
471 47 : if (nc_inq_vardimid(this->ncid, nodec_varIds[0], &index) != NC_NOERR)
472 : {
473 0 : return 0;
474 : }
475 :
476 : // Finally find the length
477 : size_t len;
478 47 : if (nc_inq_dimlen(this->ncid, index, &len) != NC_NOERR)
479 : {
480 0 : return 0;
481 : }
482 47 : return len;
483 : }
484 :
485 : else
486 482 : return this->node_counts.size();
487 : }
488 :
489 1200 : static void add_to_buffer(std::vector<unsigned char> &buffer, uint8_t v)
490 : {
491 1200 : buffer.push_back(v);
492 1200 : }
493 :
494 2899 : static void add_to_buffer(std::vector<unsigned char> &buffer, uint32_t v)
495 : {
496 2899 : const size_t old_size = buffer.size();
497 2899 : buffer.resize(old_size + sizeof(v));
498 2899 : memcpy(&buffer[old_size], &v, sizeof(v));
499 2899 : }
500 :
501 : static void add_to_buffer(std::vector<unsigned char> &, int) = delete;
502 :
503 134921 : static void add_to_buffer(std::vector<unsigned char> &buffer, double v)
504 : {
505 134921 : const size_t old_size = buffer.size();
506 134921 : buffer.resize(old_size + sizeof(v));
507 134921 : memcpy(&buffer[old_size], &v, sizeof(v));
508 134921 : }
509 :
510 : /* serializeToWKB
511 : * Takes the geometry in SGeometry at a given index and converts it into WKB
512 : * format. Converting SGeometry into WKB automatically allocates the required
513 : * buffer space and returns a buffer that MUST be free'd
514 : */
515 492 : std::vector<unsigned char> SGeometry_Reader::serializeToWKB(size_t featureInd)
516 : {
517 492 : std::vector<unsigned char> ret;
518 492 : int nc = 0;
519 492 : size_t sb = 0;
520 :
521 : // Points don't have node_count entry... only inspect and set node_counts if
522 : // not a point
523 492 : if (this->getGeometryType() != POINT)
524 : {
525 453 : nc = node_counts[featureInd];
526 453 : sb = bound_list[featureInd];
527 : }
528 :
529 : // Serialization occurs differently depending on geometry
530 : // The memory requirements also differ between geometries
531 492 : switch (this->getGeometryType())
532 : {
533 39 : case POINT:
534 39 : inPlaceSerialize_Point(this, featureInd, ret);
535 39 : break;
536 :
537 16 : case LINE:
538 16 : inPlaceSerialize_LineString(this, nc, sb, ret);
539 16 : break;
540 :
541 13 : case POLYGON:
542 : // A polygon has:
543 : // 1 byte header
544 : // 4 byte Type
545 : // 4 byte ring count (1 (exterior)) [assume only exterior rings,
546 : // otherwise multipolygon] For each ring: 4 byte point count, 8 byte
547 : // double x point count x # dimension (This is equivalent to
548 : // requesting enough for all points and a point count header for
549 : // each point) (Or 8 byte double x node count x # dimension + 4 byte
550 : // point count x part_count)
551 :
552 : // if interior ring, then assume that it will be a multipolygon
553 : // (maybe future work?)
554 13 : inPlaceSerialize_PolygonExtOnly(this, nc, sb, ret);
555 13 : break;
556 :
557 25 : case MULTIPOINT:
558 : {
559 25 : uint8_t header = PLATFORM_HEADER;
560 32 : uint32_t t = this->touple_order == 2 ? wkbMultiPoint
561 7 : : this->touple_order == 3 ? wkbMultiPoint25D
562 : : wkbNone;
563 :
564 25 : if (t == wkbNone)
565 0 : throw SG_Exception_BadFeature();
566 :
567 : // Add metadata
568 25 : add_to_buffer(ret, header);
569 25 : add_to_buffer(ret, t);
570 25 : add_to_buffer(ret, static_cast<uint32_t>(nc));
571 :
572 : // Add points
573 119 : for (int pts = 0; pts < nc; pts++)
574 : {
575 94 : inPlaceSerialize_Point(this, static_cast<size_t>(sb + pts),
576 : ret);
577 : }
578 : }
579 :
580 25 : break;
581 :
582 18 : case MULTILINE:
583 : {
584 18 : uint8_t header = PLATFORM_HEADER;
585 25 : uint32_t t = this->touple_order == 2 ? wkbMultiLineString
586 7 : : this->touple_order == 3 ? wkbMultiLineString25D
587 : : wkbNone;
588 :
589 18 : if (t == wkbNone)
590 0 : throw SG_Exception_BadFeature();
591 18 : int32_t lc = parts_count[featureInd];
592 18 : size_t seek_begin = sb;
593 : size_t pc_begin =
594 18 : pnc_bl[featureInd]; // initialize with first part count, list
595 : // of part counts is contiguous
596 36 : std::vector<int> pnc;
597 :
598 : // Build sub vector for part_node_counts
599 : // + Calculate wkbSize
600 47 : for (int itr = 0; itr < lc; itr++)
601 : {
602 29 : pnc.push_back(pnode_counts[pc_begin + itr]);
603 : }
604 :
605 18 : size_t cur_point = seek_begin;
606 18 : size_t pcount = pnc.size();
607 :
608 : // Begin Writing
609 18 : add_to_buffer(ret, header);
610 18 : add_to_buffer(ret, t);
611 18 : add_to_buffer(ret, static_cast<uint32_t>(pcount));
612 :
613 47 : for (size_t itr = 0; itr < pcount; itr++)
614 : {
615 29 : inPlaceSerialize_LineString(this, pnc[itr], cur_point, ret);
616 29 : cur_point = pnc[itr] + cur_point;
617 : }
618 : }
619 :
620 18 : break;
621 :
622 381 : case MULTIPOLYGON:
623 : {
624 381 : uint8_t header = PLATFORM_HEADER;
625 395 : uint32_t t = this->touple_order == 2 ? wkbMultiPolygon
626 14 : : this->touple_order == 3 ? wkbMultiPolygon25D
627 : : wkbNone;
628 :
629 381 : if (t == wkbNone)
630 0 : throw SG_Exception_BadFeature();
631 381 : bool noInteriors = this->int_rings.size() == 0 ? true : false;
632 381 : int32_t rc = parts_count[featureInd];
633 381 : size_t seek_begin = sb;
634 : size_t pc_begin =
635 381 : pnc_bl[featureInd]; // initialize with first part count, list
636 : // of part counts is contiguous
637 762 : std::vector<int> pnc;
638 :
639 : // Build sub vector for part_node_counts
640 1000 : for (int itr = 0; itr < rc; itr++)
641 : {
642 619 : pnc.push_back(pnode_counts[pc_begin + itr]);
643 : }
644 :
645 : // Create Multipolygon headers
646 381 : add_to_buffer(ret, header);
647 381 : add_to_buffer(ret, t);
648 :
649 381 : if (noInteriors)
650 : {
651 221 : size_t cur_point = seek_begin;
652 221 : size_t pcount = pnc.size();
653 221 : add_to_buffer(ret, static_cast<uint32_t>(pcount));
654 :
655 640 : for (size_t itr = 0; itr < pcount; itr++)
656 : {
657 419 : inPlaceSerialize_PolygonExtOnly(this, pnc[itr], cur_point,
658 : ret);
659 419 : cur_point = pnc[itr] + cur_point;
660 : }
661 : }
662 :
663 : else
664 : {
665 160 : int32_t polys = poly_count[featureInd];
666 160 : add_to_buffer(ret, static_cast<uint32_t>(polys));
667 :
668 160 : size_t base = pnc_bl[featureInd]; // beginning of parts_count
669 : // for this multigeometry
670 160 : size_t seek = seek_begin; // beginning of node range for this
671 : // multigeometry
672 160 : size_t ir_base = base + 1;
673 :
674 : // has interior rings,
675 326 : for (int32_t itr = 0; itr < polys; itr++)
676 : {
677 166 : int rc_m = 1;
678 :
679 : // count how many parts belong to each Polygon
680 200 : while (ir_base < int_rings.size() && int_rings[ir_base])
681 : {
682 34 : rc_m++;
683 34 : ir_base++;
684 : }
685 :
686 166 : if (rc_m == 1)
687 142 : ir_base++; // single ring case
688 :
689 166 : std::vector<int> poly_parts;
690 :
691 : // Make part node count sub vector
692 366 : for (int itr_2 = 0; itr_2 < rc_m; itr_2++)
693 : {
694 200 : poly_parts.push_back(pnode_counts[base + itr_2]);
695 : }
696 :
697 166 : inPlaceSerialize_Polygon(this, poly_parts, rc_m, seek, ret);
698 :
699 : // Update seek position
700 366 : for (size_t itr_3 = 0; itr_3 < poly_parts.size(); itr_3++)
701 : {
702 200 : seek += poly_parts[itr_3];
703 : }
704 :
705 166 : base += poly_parts.size();
706 : // cppcheck-suppress redundantAssignment
707 166 : ir_base = base + 1;
708 : }
709 : }
710 : }
711 381 : break;
712 :
713 0 : default:
714 :
715 0 : throw SG_Exception_BadFeature();
716 : }
717 :
718 492 : return ret;
719 : }
720 :
721 67 : void SGeometry_PropertyScanner::open(int container_id)
722 : {
723 : // First check for container_id, if variable doesn't exist error out
724 67 : if (nc_inq_var(this->nc, container_id, nullptr, nullptr, nullptr, nullptr,
725 67 : nullptr) != NC_NOERR)
726 : {
727 0 : return; // change to exception
728 : }
729 :
730 : // Now exists, see what variables refer to this one
731 : // First get name of this container
732 : char contname[NC_MAX_NAME + 1];
733 67 : memset(contname, 0, NC_MAX_NAME + 1);
734 67 : if (nc_inq_varname(this->nc, container_id, contname) != NC_NOERR)
735 : {
736 0 : return;
737 : }
738 :
739 : // Then scan throughout the netcdfDataset if those variables
740 : // geometry_container attribute matches the container
741 67 : int varCount = 0;
742 67 : if (nc_inq_nvars(this->nc, &varCount) != NC_NOERR)
743 : {
744 0 : return;
745 : }
746 :
747 134 : std::string buf;
748 615 : for (int curr = 0; curr < varCount; curr++)
749 : {
750 669 : if (NCDFGetAttr(this->nc, curr, CF_SG_GEOMETRY, buf) != CE_None ||
751 121 : buf.empty())
752 : {
753 427 : continue;
754 : }
755 :
756 : // If matches, then establish a reference by placing this variable's
757 : // data in both vectors
758 121 : if (!strcmp(contname, buf.c_str()))
759 : {
760 218 : std::string property_name;
761 :
762 : // look for special OGR original name field
763 109 : if (NCDFGetAttr(this->nc, curr, OGR_SG_ORIGINAL_LAYERNAME,
764 109 : property_name) != CE_None)
765 : {
766 : // if it doesn't exist, then just use the variable name
767 109 : property_name.resize(NC_MAX_NAME + 1, 0);
768 :
769 109 : if (nc_inq_varname(this->nc, curr, property_name.data()) !=
770 : NC_NOERR)
771 : {
772 0 : throw SG_Exception_General_Malformed(contname);
773 : }
774 : }
775 :
776 109 : v_ids.push_back(curr);
777 109 : v_headers.push_back(property_name);
778 : }
779 : }
780 : }
781 :
782 : // Exception Class Implementations
783 1 : SG_Exception_Dim_MM::SG_Exception_Dim_MM(const char *container_name,
784 : const char *field_1,
785 1 : const char *field_2)
786 : {
787 2 : std::string cn_s(container_name);
788 2 : std::string field1_s(field_1);
789 1 : std::string field2_s(field_2);
790 :
791 2 : this->err_msg = "[" + cn_s + "] One or more dimensions of " + field1_s +
792 1 : " and " + field2_s + " do not match but must match.";
793 1 : }
794 :
795 2 : SG_Exception_Existential::SG_Exception_Existential(const char *container_name,
796 2 : const char *missing_name)
797 : {
798 4 : std::string cn_s(container_name);
799 2 : std::string mn_s(missing_name);
800 :
801 4 : this->err_msg = "[" + cn_s +
802 4 : "] The property or the variable associated with " + mn_s +
803 2 : " is missing.";
804 2 : }
805 :
806 0 : SG_Exception_Dep::SG_Exception_Dep(const char *container_name, const char *arg1,
807 0 : const char *arg2)
808 : {
809 0 : std::string cn_s(container_name);
810 0 : std::string arg1_s(arg1);
811 0 : std::string arg2_s(arg2);
812 :
813 0 : this->err_msg = "[" + cn_s + "] The attribute " + arg1_s +
814 0 : " may not exist without the attribute " + arg2_s +
815 0 : " existing.";
816 0 : }
817 :
818 0 : SG_Exception_BadSum::SG_Exception_BadSum(const char *container_name,
819 0 : const char *arg1, const char *arg2)
820 : {
821 0 : std::string cn_s(container_name);
822 0 : std::string arg1_s(arg1);
823 0 : std::string arg2_s(arg2);
824 :
825 0 : this->err_msg = "[" + cn_s + "]" + " The sum of all values in " + arg1_s +
826 0 : " and " + arg2_s + " do not match.";
827 0 : }
828 :
829 0 : SG_Exception_General_Malformed ::SG_Exception_General_Malformed(const char *arg)
830 : {
831 0 : std::string arg1_s(arg);
832 :
833 : this->err_msg =
834 0 : "Corruption or malformed formatting has been detected in: " + arg1_s;
835 0 : }
836 :
837 : // to get past linker
838 7 : SG_Exception::~SG_Exception()
839 : {
840 7 : }
841 :
842 : // Helpers
843 : // following is a short hand for a clean up and exit, since goto isn't allowed
844 515 : double getCFVersion(int ncid)
845 : {
846 515 : double ver = -1.0;
847 1030 : std::string attrVal;
848 :
849 : // Fetch the CF attribute
850 971 : if (NCDFGetAttr(ncid, NC_GLOBAL, NCDF_CONVENTIONS, attrVal) != CE_None ||
851 456 : attrVal.empty())
852 : {
853 59 : return ver;
854 : }
855 :
856 456 : if (sscanf(attrVal.c_str(), "CF-%lf", &ver) != 1)
857 : {
858 33 : return -1.0;
859 : }
860 :
861 423 : return ver;
862 : }
863 :
864 73 : geom_t getGeometryType(int ncid, int varid)
865 : {
866 73 : geom_t ret = UNSUPPORTED;
867 146 : std::string gt_name_s;
868 :
869 145 : if (NCDFGetAttr(ncid, varid, CF_SG_GEOMETRY_TYPE, gt_name_s) != CE_None ||
870 72 : gt_name_s.empty())
871 : {
872 1 : return NONE;
873 : }
874 :
875 : // Points
876 72 : if (!strcmp(gt_name_s.c_str(), CF_SG_TYPE_POINT))
877 : {
878 : // Node Count not present? Assume that it is a multipoint.
879 18 : if (nc_inq_att(ncid, varid, CF_SG_NODE_COUNT, nullptr, nullptr) ==
880 : NC_ENOTATT)
881 : {
882 12 : ret = POINT;
883 : }
884 : else
885 6 : ret = MULTIPOINT;
886 : }
887 :
888 : // Lines
889 54 : else if (!strcmp(gt_name_s.c_str(), CF_SG_TYPE_LINE))
890 : {
891 : // Part Node Count present? Assume multiline
892 10 : if (nc_inq_att(ncid, varid, CF_SG_PART_NODE_COUNT, nullptr, nullptr) ==
893 : NC_ENOTATT)
894 : {
895 4 : ret = LINE;
896 : }
897 : else
898 6 : ret = MULTILINE;
899 : }
900 :
901 : // Polygons
902 44 : else if (!strcmp(gt_name_s.c_str(), CF_SG_TYPE_POLY))
903 : {
904 : /* Polygons versus MultiPolygons, slightly ambiguous
905 : * Part Node Count & no Interior Ring - MultiPolygon
906 : * no Part Node Count & no Interior Ring - Polygon
907 : * Part Node Count & Interior Ring - assume that it is a MultiPolygon
908 : */
909 : int pnc_present =
910 43 : nc_inq_att(ncid, varid, CF_SG_PART_NODE_COUNT, nullptr, nullptr);
911 : int ir_present =
912 43 : nc_inq_att(ncid, varid, CF_SG_INTERIOR_RING, nullptr, nullptr);
913 :
914 43 : if (pnc_present == NC_ENOTATT && ir_present == NC_ENOTATT)
915 : {
916 10 : ret = POLYGON;
917 : }
918 : else
919 33 : ret = MULTIPOLYGON;
920 : }
921 :
922 72 : return ret;
923 : }
924 :
925 133 : void inPlaceSerialize_Point(SGeometry_Reader *ge, size_t seek_pos,
926 : std::vector<unsigned char> &buffer)
927 : {
928 133 : uint8_t order = PLATFORM_HEADER;
929 187 : uint32_t t = ge->get_axisCount() == 2 ? wkbPoint
930 54 : : ge->get_axisCount() == 3 ? wkbPoint25D
931 133 : : wkbNone;
932 :
933 133 : if (t == wkbNone)
934 0 : throw SG_Exception_BadFeature();
935 :
936 133 : add_to_buffer(buffer, order);
937 133 : add_to_buffer(buffer, t);
938 :
939 : // Now get point data;
940 133 : Point &p = (*ge)[seek_pos];
941 133 : add_to_buffer(buffer, p[0]);
942 133 : add_to_buffer(buffer, p[1]);
943 :
944 133 : if (ge->get_axisCount() >= 3)
945 : {
946 54 : add_to_buffer(buffer, p[2]);
947 : }
948 133 : }
949 :
950 45 : void inPlaceSerialize_LineString(SGeometry_Reader *ge, int node_count,
951 : size_t seek_begin,
952 : std::vector<unsigned char> &buffer)
953 : {
954 45 : uint8_t order = PLATFORM_HEADER;
955 65 : uint32_t t = ge->get_axisCount() == 2 ? wkbLineString
956 20 : : ge->get_axisCount() == 3 ? wkbLineString25D
957 45 : : wkbNone;
958 :
959 45 : if (t == wkbNone)
960 0 : throw SG_Exception_BadFeature();
961 45 : uint32_t nc = (uint32_t)node_count;
962 :
963 45 : add_to_buffer(buffer, order);
964 45 : add_to_buffer(buffer, t);
965 45 : add_to_buffer(buffer, nc);
966 :
967 : // Now serialize points
968 170 : for (int ind = 0; ind < node_count; ind++)
969 : {
970 125 : Point &p = (*ge)[seek_begin + ind];
971 125 : add_to_buffer(buffer, p[0]);
972 125 : add_to_buffer(buffer, p[1]);
973 :
974 125 : if (ge->get_axisCount() >= 3)
975 : {
976 56 : add_to_buffer(buffer, p[2]);
977 : }
978 : }
979 45 : }
980 :
981 432 : void inPlaceSerialize_PolygonExtOnly(SGeometry_Reader *ge, int node_count,
982 : size_t seek_begin,
983 : std::vector<unsigned char> &buffer)
984 : {
985 432 : uint8_t order = PLATFORM_HEADER;
986 444 : uint32_t t = ge->get_axisCount() == 2 ? wkbPolygon
987 12 : : ge->get_axisCount() == 3 ? wkbPolygon25D
988 432 : : wkbNone;
989 :
990 432 : if (t == wkbNone)
991 0 : throw SG_Exception_BadFeature();
992 432 : uint32_t rc = 1;
993 :
994 432 : add_to_buffer(buffer, order);
995 432 : add_to_buffer(buffer, t);
996 432 : add_to_buffer(buffer, rc);
997 432 : add_to_buffer(buffer, static_cast<uint32_t>(node_count));
998 47383 : for (int pind = 0; pind < node_count; pind++)
999 : {
1000 46951 : Point &pt = (*ge)[seek_begin + pind];
1001 46951 : add_to_buffer(buffer, pt[0]);
1002 46951 : add_to_buffer(buffer, pt[1]);
1003 :
1004 46951 : if (ge->get_axisCount() >= 3)
1005 : {
1006 51 : add_to_buffer(buffer, pt[2]);
1007 : }
1008 : }
1009 432 : }
1010 :
1011 166 : void inPlaceSerialize_Polygon(SGeometry_Reader *ge, std::vector<int> &pnc,
1012 : int ring_count, size_t seek_begin,
1013 : std::vector<unsigned char> &buffer)
1014 : {
1015 166 : uint8_t order = PLATFORM_HEADER;
1016 180 : uint32_t t = ge->get_axisCount() == 2 ? wkbPolygon
1017 14 : : ge->get_axisCount() == 3 ? wkbPolygon25D
1018 166 : : wkbNone;
1019 :
1020 166 : if (t == wkbNone)
1021 0 : throw SG_Exception_BadFeature();
1022 166 : uint32_t rc = static_cast<uint32_t>(ring_count);
1023 :
1024 166 : add_to_buffer(buffer, order);
1025 166 : add_to_buffer(buffer, t);
1026 166 : add_to_buffer(buffer, rc);
1027 166 : int cmoffset = 0;
1028 366 : for (int ring_c = 0; ring_c < ring_count; ring_c++)
1029 : {
1030 200 : uint32_t node_count = pnc[ring_c];
1031 200 : add_to_buffer(buffer, node_count);
1032 :
1033 200 : int pind = 0;
1034 20332 : for (pind = 0; pind < pnc[ring_c]; pind++)
1035 : {
1036 20132 : Point &pt = (*ge)[seek_begin + cmoffset + pind];
1037 20132 : add_to_buffer(buffer, pt[0]);
1038 20132 : add_to_buffer(buffer, pt[1]);
1039 :
1040 20132 : if (ge->get_axisCount() >= 3)
1041 : {
1042 78 : add_to_buffer(buffer, pt[2]);
1043 : }
1044 : }
1045 :
1046 200 : cmoffset += pind;
1047 : }
1048 166 : }
1049 :
1050 75 : int scanForGeometryContainers(int ncid, std::set<int> &r_ids)
1051 : {
1052 : int nvars;
1053 75 : if (nc_inq_nvars(ncid, &nvars) != NC_NOERR)
1054 : {
1055 0 : return -1;
1056 : }
1057 :
1058 75 : r_ids.clear();
1059 :
1060 : // For each variable check for a geometry attribute
1061 : // If it has geometry attribute, then check the associated variable ID
1062 75 : std::string buf;
1063 703 : for (int itr = 0; itr < nvars; itr++)
1064 : {
1065 628 : if (NCDFGetAttr(ncid, itr, CF_SG_GEOMETRY, buf) != CE_None)
1066 : {
1067 513 : continue;
1068 : }
1069 :
1070 : int varID;
1071 115 : if (nc_inq_varid(ncid, buf.c_str(), &varID) != NC_NOERR)
1072 : {
1073 0 : continue;
1074 : }
1075 :
1076 : // Now have variable ID. See if vector contains it, and if not
1077 : // insert
1078 115 : r_ids.insert(varID); // It's a set. No big deal sets don't allow
1079 : // duplicates anyways
1080 : }
1081 :
1082 75 : return 0;
1083 : }
1084 : } // namespace nccfdriver
|
