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 71 : SGeometry_Reader::SGeometry_Reader(int ncId, int geoVarId)
28 137 : : gc_varId(geoVarId), touple_order(0)
29 : {
30 :
31 : char container_name[NC_MAX_NAME + 1];
32 71 : memset(container_name, 0, NC_MAX_NAME + 1);
33 :
34 : // Get geometry container name
35 71 : 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 71 : container_name_s = std::string(container_name);
43 :
44 : // Find geometry type
45 71 : this->type = nccfdriver::getGeometryType(ncId, geoVarId);
46 :
47 71 : 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 70 : this->gm_varId = INVALID_VAR_ID;
55 116 : 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 140 : std::string nc_name_s;
68 140 : std::string pnc_name_s;
69 140 : std::string ir_name_s;
70 70 : int pnc_vid = INVALID_VAR_ID;
71 70 : int nc_vid = INVALID_VAR_ID;
72 70 : int ir_vid = INVALID_VAR_ID;
73 : int buf;
74 70 : size_t bound = 0;
75 70 : size_t total_node_count = 0; // used in error checks later
76 70 : size_t total_part_node_count = 0;
77 :
78 128 : 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 69 : if (NCDFGetAttr(ncId, geoVarId, CF_SG_PART_NODE_COUNT, pnc_name_s) ==
99 108 : CE_None &&
100 39 : !pnc_name_s.empty())
101 : {
102 39 : const char *pnc_name = pnc_name_s.c_str();
103 39 : bound = 0;
104 39 : nc_inq_varid(ncId, pnc_name, &pnc_vid);
105 881 : while (nc_get_var1_int(ncId, pnc_vid, &bound, &buf) == NC_NOERR)
106 : {
107 843 : 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 842 : this->pnode_counts.push_back(buf);
115 842 : total_part_node_count += buf;
116 842 : bound++;
117 : }
118 : }
119 :
120 68 : if (NCDFGetAttr(ncId, geoVarId, CF_SG_INTERIOR_RING, ir_name_s) ==
121 80 : 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 67 : 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 67 : 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 67 : 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 67 : 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 67 : if (this->type == POLYGON || this->type == LINE)
184 : {
185 13 : 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 66 : size_t rc = 0;
197 66 : bound_list.push_back(0); // start with 0
198 :
199 66 : 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 67 : std::string cart_s;
209 : // Node Coordinates
210 66 : if (NCDFGetAttr(ncId, geoVarId, CF_SG_NODE_COORDINATES, cart_s) !=
211 132 : CE_None ||
212 66 : 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 66 : 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 893 : for (size_t pcnt = 0; pcnt < pnode_counts.size(); pcnt++)
226 : {
227 838 : if (prog == 0)
228 531 : pnc_bl.push_back(pcnt);
229 :
230 838 : if (int_rings.size() > 0 && !int_rings[pcnt])
231 236 : c++;
232 :
233 838 : prog = prog + pnode_counts[pcnt];
234 838 : parts++;
235 :
236 838 : if (prog == node_counts[ind])
237 : {
238 531 : ind++;
239 531 : this->parts_count.push_back(parts);
240 531 : if (int_rings.size() > 0)
241 230 : this->poly_count.push_back(c);
242 531 : c = 0;
243 531 : prog = 0;
244 531 : 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 66 : int X = INVALID_VAR_ID;
257 66 : int Y = INVALID_VAR_ID;
258 66 : int Z = INVALID_VAR_ID;
259 :
260 : char cart[NC_MAX_NAME + 1];
261 66 : memset(cart, 0, NC_MAX_NAME + 1);
262 66 : strncpy(cart, cart_s.c_str(), NC_MAX_NAME);
263 :
264 66 : char *dim = strtok(cart, " ");
265 66 : int axis_id = 0;
266 :
267 220 : while (dim != nullptr)
268 : {
269 154 : if (nc_inq_varid(ncId, dim, &axis_id) == NC_NOERR)
270 : {
271 :
272 : // Check axis signature
273 154 : std::string a_sig;
274 154 : CPL_IGNORE_RET_VAL(NCDFGetAttr(ncId, axis_id, CF_AXIS, a_sig));
275 :
276 : // If valid signify axis correctly
277 154 : if (a_sig == "X")
278 : {
279 66 : X = axis_id;
280 : }
281 88 : else if (a_sig == "Y")
282 : {
283 66 : Y = axis_id;
284 : }
285 22 : else if (a_sig == "Z")
286 : {
287 22 : 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 154 : this->touple_order++;
296 : }
297 : else
298 : {
299 0 : throw SG_Exception_Existential(container_name, dim);
300 : }
301 :
302 154 : dim = strtok(nullptr, " ");
303 : }
304 :
305 : // Write axis in X, Y, Z order
306 :
307 66 : if (X != INVALID_VAR_ID)
308 66 : this->nodec_varIds.push_back(X);
309 : else
310 : {
311 : throw SG_Exception_Existential(container_name,
312 0 : "node_coordinates: X-axis");
313 : }
314 66 : if (Y != INVALID_VAR_ID)
315 66 : this->nodec_varIds.push_back(Y);
316 : else
317 : {
318 : throw SG_Exception_Existential(container_name,
319 0 : "node_coordinates: Y-axis");
320 : }
321 66 : if (Z != INVALID_VAR_ID)
322 22 : 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 66 : int all_dim = INVALID_VAR_ID;
332 66 : bool dim_set = false;
333 : //(1) one dimension check, each node_coordinates have same dimension
334 219 : for (size_t nvitr = 0; nvitr < nodec_varIds.size(); nvitr++)
335 : {
336 154 : int dimC = 0;
337 154 : nc_inq_varndims(ncId, nodec_varIds[nvitr], &dimC);
338 :
339 154 : 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 154 : 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 154 : if (!dim_set)
353 : {
354 66 : all_dim = inter_dim[0];
355 66 : dim_set = true;
356 : }
357 :
358 : else
359 : {
360 88 : 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 65 : 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 65 : 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 65 : size_t instance_dim_len = 0;
394 :
395 65 : 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 11 : this->inst_dimId = all_dim;
417 : }
418 :
419 65 : nc_inq_dimlen(ncId, this->inst_dimId, &instance_dim_len);
420 :
421 65 : if (instance_dim_len == 0)
422 0 : throw SG_Exception_EmptyDim();
423 :
424 : // Set values accordingly
425 65 : this->inst_dimLen = instance_dim_len;
426 65 : this->pt_buffer = std::make_unique<Point>(this->touple_order);
427 65 : this->gc_varId = geoVarId;
428 65 : this->ncid = ncId;
429 65 : }
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 65 : void SGeometry_PropertyScanner::open(int container_id)
722 : {
723 : // First check for container_id, if variable doesn't exist error out
724 65 : if (nc_inq_var(this->nc, container_id, nullptr, nullptr, nullptr, nullptr,
725 65 : 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 65 : memset(contname, 0, NC_MAX_NAME + 1);
734 65 : 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 65 : int varCount = 0;
742 65 : if (nc_inq_nvars(this->nc, &varCount) != NC_NOERR)
743 : {
744 0 : return;
745 : }
746 :
747 130 : std::string buf;
748 598 : for (int curr = 0; curr < varCount; curr++)
749 : {
750 646 : if (NCDFGetAttr(this->nc, curr, CF_SG_GEOMETRY, buf) != CE_None ||
751 113 : buf.empty())
752 : {
753 420 : continue;
754 : }
755 :
756 : // If matches, then establish a reference by placing this variable's
757 : // data in both vectors
758 113 : if (!strcmp(contname, buf.c_str()))
759 : {
760 202 : std::string property_name;
761 :
762 : // look for special OGR original name field
763 101 : if (NCDFGetAttr(this->nc, curr, OGR_SG_ORIGINAL_LAYERNAME,
764 101 : property_name) != CE_None)
765 : {
766 : // if it doesn't exist, then just use the variable name
767 101 : property_name.resize(NC_MAX_NAME + 1, 0);
768 :
769 101 : 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 101 : v_ids.push_back(curr);
777 101 : v_headers.push_back(std::move(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 552 : bool getCFVersion(int ncid, int &major, int &minor)
845 : {
846 1104 : std::string attrVal;
847 552 : major = 0;
848 552 : minor = 0;
849 :
850 : // Fetch the CF attribute
851 1044 : if (NCDFGetAttr(ncid, NC_GLOBAL, NCDF_CONVENTIONS, attrVal) != CE_None ||
852 492 : attrVal.empty())
853 : {
854 60 : return false;
855 : }
856 :
857 492 : if (STARTS_WITH(attrVal.c_str(), "CF-"))
858 : {
859 : const CPLStringList aosTokens(
860 459 : CSLTokenizeString2(attrVal.c_str() + strlen("CF-"), ".", 0));
861 459 : if (aosTokens.size() == 2)
862 : {
863 457 : major = atoi(aosTokens[0]);
864 457 : minor = atoi(aosTokens[1]);
865 457 : return true;
866 : }
867 : }
868 :
869 35 : return false;
870 : }
871 :
872 71 : geom_t getGeometryType(int ncid, int varid)
873 : {
874 71 : geom_t ret = UNSUPPORTED;
875 142 : std::string gt_name_s;
876 :
877 141 : if (NCDFGetAttr(ncid, varid, CF_SG_GEOMETRY_TYPE, gt_name_s) != CE_None ||
878 70 : gt_name_s.empty())
879 : {
880 1 : return NONE;
881 : }
882 :
883 : // Points
884 70 : if (!strcmp(gt_name_s.c_str(), CF_SG_TYPE_POINT))
885 : {
886 : // Node Count not present? Assume that it is a multipoint.
887 16 : if (nc_inq_att(ncid, varid, CF_SG_NODE_COUNT, nullptr, nullptr) ==
888 : NC_ENOTATT)
889 : {
890 10 : ret = POINT;
891 : }
892 : else
893 6 : ret = MULTIPOINT;
894 : }
895 :
896 : // Lines
897 54 : else if (!strcmp(gt_name_s.c_str(), CF_SG_TYPE_LINE))
898 : {
899 : // Part Node Count present? Assume multiline
900 10 : if (nc_inq_att(ncid, varid, CF_SG_PART_NODE_COUNT, nullptr, nullptr) ==
901 : NC_ENOTATT)
902 : {
903 4 : ret = LINE;
904 : }
905 : else
906 6 : ret = MULTILINE;
907 : }
908 :
909 : // Polygons
910 44 : else if (!strcmp(gt_name_s.c_str(), CF_SG_TYPE_POLY))
911 : {
912 : /* Polygons versus MultiPolygons, slightly ambiguous
913 : * Part Node Count & no Interior Ring - MultiPolygon
914 : * no Part Node Count & no Interior Ring - Polygon
915 : * Part Node Count & Interior Ring - assume that it is a MultiPolygon
916 : */
917 : int pnc_present =
918 43 : nc_inq_att(ncid, varid, CF_SG_PART_NODE_COUNT, nullptr, nullptr);
919 : int ir_present =
920 43 : nc_inq_att(ncid, varid, CF_SG_INTERIOR_RING, nullptr, nullptr);
921 :
922 43 : if (pnc_present == NC_ENOTATT && ir_present == NC_ENOTATT)
923 : {
924 9 : ret = POLYGON;
925 : }
926 : else
927 34 : ret = MULTIPOLYGON;
928 : }
929 :
930 70 : return ret;
931 : }
932 :
933 133 : void inPlaceSerialize_Point(SGeometry_Reader *ge, size_t seek_pos,
934 : std::vector<unsigned char> &buffer)
935 : {
936 133 : uint8_t order = PLATFORM_HEADER;
937 187 : uint32_t t = ge->get_axisCount() == 2 ? wkbPoint
938 54 : : ge->get_axisCount() == 3 ? wkbPoint25D
939 133 : : wkbNone;
940 :
941 133 : if (t == wkbNone)
942 0 : throw SG_Exception_BadFeature();
943 :
944 133 : add_to_buffer(buffer, order);
945 133 : add_to_buffer(buffer, t);
946 :
947 : // Now get point data;
948 133 : Point &p = (*ge)[seek_pos];
949 133 : add_to_buffer(buffer, p[0]);
950 133 : add_to_buffer(buffer, p[1]);
951 :
952 133 : if (ge->get_axisCount() >= 3)
953 : {
954 54 : add_to_buffer(buffer, p[2]);
955 : }
956 133 : }
957 :
958 45 : void inPlaceSerialize_LineString(SGeometry_Reader *ge, int node_count,
959 : size_t seek_begin,
960 : std::vector<unsigned char> &buffer)
961 : {
962 45 : uint8_t order = PLATFORM_HEADER;
963 65 : uint32_t t = ge->get_axisCount() == 2 ? wkbLineString
964 20 : : ge->get_axisCount() == 3 ? wkbLineString25D
965 45 : : wkbNone;
966 :
967 45 : if (t == wkbNone)
968 0 : throw SG_Exception_BadFeature();
969 45 : uint32_t nc = (uint32_t)node_count;
970 :
971 45 : add_to_buffer(buffer, order);
972 45 : add_to_buffer(buffer, t);
973 45 : add_to_buffer(buffer, nc);
974 :
975 : // Now serialize points
976 170 : for (int ind = 0; ind < node_count; ind++)
977 : {
978 125 : Point &p = (*ge)[seek_begin + ind];
979 125 : add_to_buffer(buffer, p[0]);
980 125 : add_to_buffer(buffer, p[1]);
981 :
982 125 : if (ge->get_axisCount() >= 3)
983 : {
984 56 : add_to_buffer(buffer, p[2]);
985 : }
986 : }
987 45 : }
988 :
989 432 : void inPlaceSerialize_PolygonExtOnly(SGeometry_Reader *ge, int node_count,
990 : size_t seek_begin,
991 : std::vector<unsigned char> &buffer)
992 : {
993 432 : uint8_t order = PLATFORM_HEADER;
994 444 : uint32_t t = ge->get_axisCount() == 2 ? wkbPolygon
995 12 : : ge->get_axisCount() == 3 ? wkbPolygon25D
996 432 : : wkbNone;
997 :
998 432 : if (t == wkbNone)
999 0 : throw SG_Exception_BadFeature();
1000 432 : uint32_t rc = 1;
1001 :
1002 432 : add_to_buffer(buffer, order);
1003 432 : add_to_buffer(buffer, t);
1004 432 : add_to_buffer(buffer, rc);
1005 432 : add_to_buffer(buffer, static_cast<uint32_t>(node_count));
1006 47383 : for (int pind = 0; pind < node_count; pind++)
1007 : {
1008 46951 : Point &pt = (*ge)[seek_begin + pind];
1009 46951 : add_to_buffer(buffer, pt[0]);
1010 46951 : add_to_buffer(buffer, pt[1]);
1011 :
1012 46951 : if (ge->get_axisCount() >= 3)
1013 : {
1014 51 : add_to_buffer(buffer, pt[2]);
1015 : }
1016 : }
1017 432 : }
1018 :
1019 166 : void inPlaceSerialize_Polygon(SGeometry_Reader *ge, std::vector<int> &pnc,
1020 : int ring_count, size_t seek_begin,
1021 : std::vector<unsigned char> &buffer)
1022 : {
1023 166 : uint8_t order = PLATFORM_HEADER;
1024 180 : uint32_t t = ge->get_axisCount() == 2 ? wkbPolygon
1025 14 : : ge->get_axisCount() == 3 ? wkbPolygon25D
1026 166 : : wkbNone;
1027 :
1028 166 : if (t == wkbNone)
1029 0 : throw SG_Exception_BadFeature();
1030 166 : uint32_t rc = static_cast<uint32_t>(ring_count);
1031 :
1032 166 : add_to_buffer(buffer, order);
1033 166 : add_to_buffer(buffer, t);
1034 166 : add_to_buffer(buffer, rc);
1035 166 : int cmoffset = 0;
1036 366 : for (int ring_c = 0; ring_c < ring_count; ring_c++)
1037 : {
1038 200 : uint32_t node_count = pnc[ring_c];
1039 200 : add_to_buffer(buffer, node_count);
1040 :
1041 200 : int pind = 0;
1042 20332 : for (pind = 0; pind < pnc[ring_c]; pind++)
1043 : {
1044 20132 : Point &pt = (*ge)[seek_begin + cmoffset + pind];
1045 20132 : add_to_buffer(buffer, pt[0]);
1046 20132 : add_to_buffer(buffer, pt[1]);
1047 :
1048 20132 : if (ge->get_axisCount() >= 3)
1049 : {
1050 78 : add_to_buffer(buffer, pt[2]);
1051 : }
1052 : }
1053 :
1054 200 : cmoffset += pind;
1055 : }
1056 166 : }
1057 :
1058 73 : int scanForGeometryContainers(int ncid, std::set<int> &r_ids)
1059 : {
1060 : int nvars;
1061 73 : if (nc_inq_nvars(ncid, &nvars) != NC_NOERR)
1062 : {
1063 0 : return -1;
1064 : }
1065 :
1066 73 : r_ids.clear();
1067 :
1068 : // For each variable check for a geometry attribute
1069 : // If it has geometry attribute, then check the associated variable ID
1070 73 : std::string buf;
1071 686 : for (int itr = 0; itr < nvars; itr++)
1072 : {
1073 613 : if (NCDFGetAttr(ncid, itr, CF_SG_GEOMETRY, buf) != CE_None)
1074 : {
1075 506 : continue;
1076 : }
1077 :
1078 : int varID;
1079 107 : if (nc_inq_varid(ncid, buf.c_str(), &varID) != NC_NOERR)
1080 : {
1081 0 : continue;
1082 : }
1083 :
1084 : // Now have variable ID. See if vector contains it, and if not
1085 : // insert
1086 107 : r_ids.insert(varID); // It's a set. No big deal sets don't allow
1087 : // duplicates anyways
1088 : }
1089 :
1090 73 : return 0;
1091 : }
1092 : } // namespace nccfdriver
|
