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