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