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 1067 : std::string &attrf(int ncid, int varId, const char *attrName,
33 : std::string &alloc)
34 : {
35 1067 : size_t len = 0;
36 1067 : nc_inq_attlen(ncid, varId, attrName, &len);
37 :
38 1067 : if (len < 1)
39 : {
40 188 : alloc.clear();
41 188 : return alloc;
42 : }
43 :
44 879 : alloc.resize(len);
45 879 : memset(&alloc[0], 0, len);
46 :
47 : // Now look through this variable for the attribute
48 879 : nc_get_att_text(ncid, varId, attrName, &alloc[0]);
49 :
50 879 : 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 218 : std::string n(property_name);
812 109 : v_ids.push_back(curr);
813 109 : v_headers.push_back(n);
814 : }
815 : }
816 : }
817 :
818 : // Exception Class Implementations
819 1 : SG_Exception_Dim_MM::SG_Exception_Dim_MM(const char *container_name,
820 : const char *field_1,
821 1 : const char *field_2)
822 : {
823 2 : std::string cn_s(container_name);
824 2 : std::string field1_s(field_1);
825 1 : std::string field2_s(field_2);
826 :
827 2 : this->err_msg = "[" + cn_s + "] One or more dimensions of " + field1_s +
828 1 : " and " + field2_s + " do not match but must match.";
829 1 : }
830 :
831 2 : SG_Exception_Existential::SG_Exception_Existential(const char *container_name,
832 2 : const char *missing_name)
833 : {
834 4 : std::string cn_s(container_name);
835 2 : std::string mn_s(missing_name);
836 :
837 4 : this->err_msg = "[" + cn_s +
838 4 : "] The property or the variable associated with " + mn_s +
839 2 : " is missing.";
840 2 : }
841 :
842 0 : SG_Exception_Dep::SG_Exception_Dep(const char *container_name, const char *arg1,
843 0 : const char *arg2)
844 : {
845 0 : std::string cn_s(container_name);
846 0 : std::string arg1_s(arg1);
847 0 : std::string arg2_s(arg2);
848 :
849 0 : this->err_msg = "[" + cn_s + "] The attribute " + arg1_s +
850 0 : " may not exist without the attribute " + arg2_s +
851 0 : " existing.";
852 0 : }
853 :
854 0 : SG_Exception_BadSum::SG_Exception_BadSum(const char *container_name,
855 0 : const char *arg1, const char *arg2)
856 : {
857 0 : std::string cn_s(container_name);
858 0 : std::string arg1_s(arg1);
859 0 : std::string arg2_s(arg2);
860 :
861 0 : this->err_msg = "[" + cn_s + "]" + " The sum of all values in " + arg1_s +
862 0 : " and " + arg2_s + " do not match.";
863 0 : }
864 :
865 0 : SG_Exception_General_Malformed ::SG_Exception_General_Malformed(const char *arg)
866 : {
867 0 : std::string arg1_s(arg);
868 :
869 : this->err_msg =
870 0 : "Corruption or malformed formatting has been detected in: " + arg1_s;
871 0 : }
872 :
873 : // to get past linker
874 7 : SG_Exception::~SG_Exception()
875 : {
876 7 : }
877 :
878 : // Helpers
879 : // following is a short hand for a clean up and exit, since goto isn't allowed
880 481 : double getCFVersion(int ncid)
881 : {
882 481 : double ver = -1.0;
883 962 : std::string attrVal;
884 :
885 : // Fetch the CF attribute
886 481 : if (attrf(ncid, NC_GLOBAL, NCDF_CONVENTIONS, attrVal) == "")
887 : {
888 56 : return ver;
889 : }
890 :
891 425 : if (sscanf(attrVal.c_str(), "CF-%lf", &ver) != 1)
892 : {
893 31 : return -1.0;
894 : }
895 :
896 394 : return ver;
897 : }
898 :
899 73 : geom_t getGeometryType(int ncid, int varid)
900 : {
901 73 : geom_t ret = UNSUPPORTED;
902 146 : std::string gt_name_s;
903 : const char *gt_name =
904 73 : attrf(ncid, varid, CF_SG_GEOMETRY_TYPE, gt_name_s).c_str();
905 :
906 73 : if (gt_name[0] == '\0')
907 : {
908 1 : return NONE;
909 : }
910 :
911 : // Points
912 72 : if (!strcmp(gt_name, CF_SG_TYPE_POINT))
913 : {
914 : // Node Count not present? Assume that it is a multipoint.
915 18 : if (nc_inq_att(ncid, varid, CF_SG_NODE_COUNT, nullptr, nullptr) ==
916 : NC_ENOTATT)
917 : {
918 12 : ret = POINT;
919 : }
920 : else
921 6 : ret = MULTIPOINT;
922 : }
923 :
924 : // Lines
925 54 : else if (!strcmp(gt_name, CF_SG_TYPE_LINE))
926 : {
927 : // Part Node Count present? Assume multiline
928 10 : if (nc_inq_att(ncid, varid, CF_SG_PART_NODE_COUNT, nullptr, nullptr) ==
929 : NC_ENOTATT)
930 : {
931 4 : ret = LINE;
932 : }
933 : else
934 6 : ret = MULTILINE;
935 : }
936 :
937 : // Polygons
938 44 : else if (!strcmp(gt_name, CF_SG_TYPE_POLY))
939 : {
940 : /* Polygons versus MultiPolygons, slightly ambiguous
941 : * Part Node Count & no Interior Ring - MultiPolygon
942 : * no Part Node Count & no Interior Ring - Polygon
943 : * Part Node Count & Interior Ring - assume that it is a MultiPolygon
944 : */
945 : int pnc_present =
946 43 : nc_inq_att(ncid, varid, CF_SG_PART_NODE_COUNT, nullptr, nullptr);
947 : int ir_present =
948 43 : nc_inq_att(ncid, varid, CF_SG_INTERIOR_RING, nullptr, nullptr);
949 :
950 43 : if (pnc_present == NC_ENOTATT && ir_present == NC_ENOTATT)
951 : {
952 10 : ret = POLYGON;
953 : }
954 : else
955 33 : ret = MULTIPOLYGON;
956 : }
957 :
958 72 : return ret;
959 : }
960 :
961 133 : void inPlaceSerialize_Point(SGeometry_Reader *ge, size_t seek_pos,
962 : std::vector<unsigned char> &buffer)
963 : {
964 133 : uint8_t order = PLATFORM_HEADER;
965 187 : uint32_t t = ge->get_axisCount() == 2 ? wkbPoint
966 54 : : ge->get_axisCount() == 3 ? wkbPoint25D
967 133 : : wkbNone;
968 :
969 133 : if (t == wkbNone)
970 0 : throw SG_Exception_BadFeature();
971 :
972 133 : add_to_buffer(buffer, order);
973 133 : add_to_buffer(buffer, t);
974 :
975 : // Now get point data;
976 133 : Point &p = (*ge)[seek_pos];
977 133 : add_to_buffer(buffer, p[0]);
978 133 : add_to_buffer(buffer, p[1]);
979 :
980 133 : if (ge->get_axisCount() >= 3)
981 : {
982 54 : add_to_buffer(buffer, p[2]);
983 : }
984 133 : }
985 :
986 45 : void inPlaceSerialize_LineString(SGeometry_Reader *ge, int node_count,
987 : size_t seek_begin,
988 : std::vector<unsigned char> &buffer)
989 : {
990 45 : uint8_t order = PLATFORM_HEADER;
991 65 : uint32_t t = ge->get_axisCount() == 2 ? wkbLineString
992 20 : : ge->get_axisCount() == 3 ? wkbLineString25D
993 45 : : wkbNone;
994 :
995 45 : if (t == wkbNone)
996 0 : throw SG_Exception_BadFeature();
997 45 : uint32_t nc = (uint32_t)node_count;
998 :
999 45 : add_to_buffer(buffer, order);
1000 45 : add_to_buffer(buffer, t);
1001 45 : add_to_buffer(buffer, nc);
1002 :
1003 : // Now serialize points
1004 170 : for (int ind = 0; ind < node_count; ind++)
1005 : {
1006 125 : Point &p = (*ge)[seek_begin + ind];
1007 125 : add_to_buffer(buffer, p[0]);
1008 125 : add_to_buffer(buffer, p[1]);
1009 :
1010 125 : if (ge->get_axisCount() >= 3)
1011 : {
1012 56 : add_to_buffer(buffer, p[2]);
1013 : }
1014 : }
1015 45 : }
1016 :
1017 432 : void inPlaceSerialize_PolygonExtOnly(SGeometry_Reader *ge, int node_count,
1018 : size_t seek_begin,
1019 : std::vector<unsigned char> &buffer)
1020 : {
1021 432 : uint8_t order = PLATFORM_HEADER;
1022 444 : uint32_t t = ge->get_axisCount() == 2 ? wkbPolygon
1023 12 : : ge->get_axisCount() == 3 ? wkbPolygon25D
1024 432 : : wkbNone;
1025 :
1026 432 : if (t == wkbNone)
1027 0 : throw SG_Exception_BadFeature();
1028 432 : uint32_t rc = 1;
1029 :
1030 432 : add_to_buffer(buffer, order);
1031 432 : add_to_buffer(buffer, t);
1032 432 : add_to_buffer(buffer, rc);
1033 432 : add_to_buffer(buffer, static_cast<uint32_t>(node_count));
1034 47383 : for (int pind = 0; pind < node_count; pind++)
1035 : {
1036 46951 : Point &pt = (*ge)[seek_begin + pind];
1037 46951 : add_to_buffer(buffer, pt[0]);
1038 46951 : add_to_buffer(buffer, pt[1]);
1039 :
1040 46951 : if (ge->get_axisCount() >= 3)
1041 : {
1042 51 : add_to_buffer(buffer, pt[2]);
1043 : }
1044 : }
1045 432 : }
1046 :
1047 166 : void inPlaceSerialize_Polygon(SGeometry_Reader *ge, std::vector<int> &pnc,
1048 : int ring_count, size_t seek_begin,
1049 : std::vector<unsigned char> &buffer)
1050 : {
1051 166 : uint8_t order = PLATFORM_HEADER;
1052 180 : uint32_t t = ge->get_axisCount() == 2 ? wkbPolygon
1053 14 : : ge->get_axisCount() == 3 ? wkbPolygon25D
1054 166 : : wkbNone;
1055 :
1056 166 : if (t == wkbNone)
1057 0 : throw SG_Exception_BadFeature();
1058 166 : uint32_t rc = static_cast<uint32_t>(ring_count);
1059 :
1060 166 : add_to_buffer(buffer, order);
1061 166 : add_to_buffer(buffer, t);
1062 166 : add_to_buffer(buffer, rc);
1063 166 : int cmoffset = 0;
1064 366 : for (int ring_c = 0; ring_c < ring_count; ring_c++)
1065 : {
1066 200 : uint32_t node_count = pnc[ring_c];
1067 200 : add_to_buffer(buffer, node_count);
1068 :
1069 200 : int pind = 0;
1070 20332 : for (pind = 0; pind < pnc[ring_c]; pind++)
1071 : {
1072 20132 : Point &pt = (*ge)[seek_begin + cmoffset + pind];
1073 20132 : add_to_buffer(buffer, pt[0]);
1074 20132 : add_to_buffer(buffer, pt[1]);
1075 :
1076 20132 : if (ge->get_axisCount() >= 3)
1077 : {
1078 78 : add_to_buffer(buffer, pt[2]);
1079 : }
1080 : }
1081 :
1082 200 : cmoffset += pind;
1083 : }
1084 166 : }
1085 :
1086 75 : int scanForGeometryContainers(int ncid, std::set<int> &r_ids)
1087 : {
1088 : int nvars;
1089 75 : if (nc_inq_nvars(ncid, &nvars) != NC_NOERR)
1090 : {
1091 0 : return -1;
1092 : }
1093 :
1094 75 : r_ids.clear();
1095 :
1096 : // For each variable check for geometry attribute
1097 : // If has geometry attribute, then check the associated variable ID
1098 :
1099 703 : for (int itr = 0; itr < nvars; itr++)
1100 : {
1101 : char c[NC_MAX_CHAR];
1102 628 : memset(c, 0, NC_MAX_CHAR);
1103 628 : if (nc_get_att_text(ncid, itr, CF_SG_GEOMETRY, c) != NC_NOERR)
1104 : {
1105 513 : continue;
1106 : }
1107 :
1108 : int varID;
1109 115 : if (nc_inq_varid(ncid, c, &varID) != NC_NOERR)
1110 : {
1111 0 : continue;
1112 : }
1113 :
1114 : // Now have variable ID. See if vector contains it, and if not
1115 : // insert
1116 115 : r_ids.insert(varID); // It's a set. No big deal sets don't allow
1117 : // duplicates anyways
1118 : }
1119 :
1120 75 : return 0;
1121 : }
1122 : } // namespace nccfdriver
|
