Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: WCS Client Driver
4 : * Purpose: Implementation of Dataset class for WCS 2.0.
5 : * Author: Ari Jolma <ari dot jolma at gmail dot com>
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2006, Frank Warmerdam
9 : * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
10 : * Copyright (c) 2017, Ari Jolma
11 : * Copyright (c) 2017, Finnish Environment Institute
12 : *
13 : * SPDX-License-Identifier: MIT
14 : ****************************************************************************/
15 :
16 : #include "cpl_string.h"
17 : #include "cpl_minixml.h"
18 : #include "cpl_http.h"
19 : #include "cpl_conv.h"
20 : #include "gmlutils.h"
21 : #include "gdal_frmts.h"
22 : #include "gdal_pam.h"
23 : #include "ogr_spatialref.h"
24 : #include "gmlcoverage.h"
25 :
26 : #include <algorithm>
27 :
28 : #include "wcsdataset.h"
29 : #include "wcsutils.h"
30 :
31 : using namespace WCSUtils;
32 :
33 : /************************************************************************/
34 : /* CoverageSubtype() */
35 : /* */
36 : /************************************************************************/
37 :
38 21 : static std::string CoverageSubtype(CPLXMLNode *coverage)
39 : {
40 : std::string subtype =
41 21 : CPLGetXMLValue(coverage, "ServiceParameters.CoverageSubtype", "");
42 21 : size_t pos = subtype.find("Coverage");
43 21 : if (pos != std::string::npos)
44 : {
45 21 : subtype.erase(pos);
46 : }
47 21 : return subtype;
48 : }
49 :
50 : /************************************************************************/
51 : /* GetGridNode() */
52 : /* */
53 : /************************************************************************/
54 :
55 21 : static CPLXMLNode *GetGridNode(CPLXMLNode *coverage, const std::string &subtype)
56 : {
57 21 : CPLXMLNode *grid = nullptr;
58 : // Construct the name of the node that we look under domainSet.
59 : // For now we can handle RectifiedGrid and ReferenceableGridByVectors.
60 : // Note that if this is called at GetCoverage stage, the grid should not be
61 : // NULL.
62 21 : std::string path = "domainSet";
63 21 : if (subtype == "RectifiedGrid")
64 : {
65 18 : grid = CPLGetXMLNode(coverage, (path + "." + subtype).c_str());
66 : }
67 3 : else if (subtype == "ReferenceableGrid")
68 : {
69 3 : grid = CPLGetXMLNode(coverage,
70 6 : (path + "." + subtype + "ByVectors").c_str());
71 : }
72 21 : if (!grid)
73 : {
74 0 : CPLError(CE_Failure, CPLE_AppDefined,
75 : "Can't handle coverages of type '%s'.", subtype.c_str());
76 : }
77 42 : return grid;
78 : }
79 :
80 : /************************************************************************/
81 : /* ParseParameters() */
82 : /* */
83 : /************************************************************************/
84 :
85 42 : static void ParseParameters(CPLXMLNode *service,
86 : std::vector<std::string> &dimensions,
87 : std::string &range,
88 : std::vector<std::vector<std::string>> &others)
89 : {
90 : std::vector<std::string> parameters =
91 84 : Split(CPLGetXMLValue(service, "Parameters", ""), "&");
92 84 : for (unsigned int i = 0; i < parameters.size(); ++i)
93 : {
94 42 : std::vector<std::string> kv = Split(parameters[i].c_str(), "=");
95 42 : if (kv.size() < 2)
96 : {
97 0 : continue;
98 : }
99 42 : kv[0] = CPLString(kv[0]).toupper();
100 42 : if (kv[0] == "RANGESUBSET")
101 : {
102 0 : range = kv[1];
103 : }
104 42 : else if (kv[0] == "SUBSET")
105 : {
106 0 : dimensions = Split(kv[1].c_str(), ";");
107 : }
108 : else
109 : {
110 126 : others.push_back(std::vector<std::string>{kv[0], kv[1]});
111 : }
112 : }
113 : // fallback to service values, if any
114 42 : if (range == "")
115 : {
116 42 : range = CPLGetXMLValue(service, "RangeSubset", "");
117 : }
118 42 : if (dimensions.size() == 0)
119 : {
120 42 : dimensions = Split(CPLGetXMLValue(service, "Subset", ""), ";");
121 : }
122 42 : }
123 :
124 : /************************************************************************/
125 : /* GetExtent() */
126 : /* */
127 : /************************************************************************/
128 :
129 21 : std::vector<double> WCSDataset201::GetExtent(int nXOff, int nYOff, int nXSize,
130 : int nYSize,
131 : CPL_UNUSED int nBufXSize,
132 : CPL_UNUSED int nBufYSize)
133 : {
134 21 : std::vector<double> extent;
135 : // WCS 2.0 extents are the outer edges of outer pixels.
136 21 : extent.push_back(adfGeoTransform[0] + (nXOff)*adfGeoTransform[1]);
137 0 : extent.push_back(adfGeoTransform[3] +
138 21 : (nYOff + nYSize) * adfGeoTransform[5]);
139 0 : extent.push_back(adfGeoTransform[0] +
140 21 : (nXOff + nXSize) * adfGeoTransform[1]);
141 21 : extent.push_back(adfGeoTransform[3] + (nYOff)*adfGeoTransform[5]);
142 21 : return extent;
143 : }
144 :
145 : /************************************************************************/
146 : /* GetCoverageRequest() */
147 : /* */
148 : /************************************************************************/
149 :
150 : std::string
151 21 : WCSDataset201::GetCoverageRequest(bool scaled, int nBufXSize, int nBufYSize,
152 : const std::vector<double> &extent,
153 : const std::string & /*osBandList*/)
154 : {
155 21 : std::string request = CPLGetXMLValue(psService, "ServiceURL", "");
156 21 : request = CPLURLAddKVP(request.c_str(), "SERVICE", "WCS");
157 21 : request += "&REQUEST=GetCoverage";
158 : request +=
159 21 : "&VERSION=" + std::string(CPLGetXMLValue(psService, "Version", ""));
160 42 : request += "&COVERAGEID=" +
161 63 : URLEncode(CPLGetXMLValue(psService, "CoverageName", ""));
162 :
163 : // note: native_crs is not really supported
164 21 : if (!native_crs)
165 : {
166 0 : std::string crs = URLEncode(CPLGetXMLValue(psService, "SRS", ""));
167 0 : request += "&OUTPUTCRS=" + crs;
168 0 : request += "&SUBSETTINGCRS=" + crs;
169 : }
170 :
171 : std::vector<std::string> domain =
172 42 : Split(CPLGetXMLValue(psService, "Domain", ""), ",");
173 21 : if (domain.size() < 2)
174 : {
175 : // eek!
176 0 : domain.push_back("E");
177 0 : domain.push_back("N");
178 : }
179 21 : const char *x = domain[0].c_str();
180 21 : const char *y = domain[1].c_str();
181 21 : if (CPLGetXMLBoolean(psService, "SubsetAxisSwap"))
182 : {
183 6 : const char *tmp = x;
184 6 : x = y;
185 6 : y = tmp;
186 : }
187 :
188 : std::vector<std::string> low =
189 42 : Split(CPLGetXMLValue(psService, "Low", ""), ",");
190 : std::vector<std::string> high =
191 42 : Split(CPLGetXMLValue(psService, "High", ""), ",");
192 42 : std::string a = CPLString().Printf("%.17g", extent[0]);
193 21 : if (low.size() > 1 && CPLAtof(low[0].c_str()) > extent[0])
194 : {
195 0 : a = low[0];
196 : }
197 42 : std::string b = CPLString().Printf("%.17g", extent[2]);
198 21 : if (high.size() > 1 && CPLAtof(high[0].c_str()) < extent[2])
199 : {
200 1 : b = high[0];
201 : }
202 : /*
203 : std::string a = CPLString().Printf(
204 : "%.17g", MAX(adfGeoTransform[0], extent[0]));
205 : std::string b = CPLString().Printf(
206 : "%.17g", MIN(adfGeoTransform[0] + nRasterXSize * adfGeoTransform[1],
207 : extent[2]));
208 : */
209 :
210 : // 09-147 KVP Protocol: subset keys must be unique
211 : // GeoServer: seems to require plain SUBSET for x and y
212 :
213 : request +=
214 21 : CPLString().Printf("&SUBSET=%s%%28%s,%s%%29", x, a.c_str(), b.c_str());
215 :
216 21 : a = CPLString().Printf("%.17g", extent[1]);
217 21 : if (low.size() > 1 && CPLAtof(low[1].c_str()) > extent[1])
218 : {
219 0 : a = low[1];
220 : }
221 21 : b = CPLString().Printf("%.17g", extent[3]);
222 21 : if (high.size() > 1 && CPLAtof(high[1].c_str()) < extent[3])
223 : {
224 0 : b = high[1];
225 : }
226 : /*
227 : a = CPLString().Printf(
228 : "%.17g", MAX(adfGeoTransform[3] + nRasterYSize * adfGeoTransform[5],
229 : extent[1])); b = CPLString().Printf(
230 : "%.17g", MIN(adfGeoTransform[3], extent[3]));
231 : */
232 :
233 : request +=
234 21 : CPLString().Printf("&SUBSET=%s%%28%s,%s%%29", y, a.c_str(), b.c_str());
235 :
236 : // Dimension and range parameters:
237 42 : std::vector<std::string> dimensions;
238 42 : std::string range;
239 42 : std::vector<std::vector<std::string>> others;
240 21 : ParseParameters(psService, dimensions, range, others);
241 :
242 : // set subsets for axis other than x/y
243 24 : for (unsigned int i = 0; i < dimensions.size(); ++i)
244 : {
245 3 : size_t pos = dimensions[i].find("(");
246 3 : std::string dim = dimensions[i].substr(0, pos);
247 3 : if (IndexOf(dim, domain) != -1)
248 : {
249 0 : continue;
250 : }
251 : std::vector<std::string> params =
252 6 : Split(FromParenthesis(dimensions[i]).c_str(), ",");
253 : request +=
254 3 : "&SUBSET" + CPLString().Printf("%i", i) + "=" + dim + "%28"; // (
255 6 : for (unsigned int j = 0; j < params.size(); ++j)
256 : {
257 : // todo: %22 (") should be used only for non-numbers
258 3 : request += "%22" + params[j] + "%22";
259 : }
260 3 : request += "%29"; // )
261 : }
262 :
263 21 : if (scaled)
264 : {
265 14 : CPLString tmp;
266 : // scaling is expressed in grid axes
267 7 : if (CPLGetXMLBoolean(psService, "UseScaleFactor"))
268 : {
269 1 : double fx = fabs((extent[2] - extent[0]) / adfGeoTransform[1] /
270 1 : ((double)nBufXSize + 0.5));
271 1 : double fy = fabs((extent[3] - extent[1]) / adfGeoTransform[5] /
272 1 : ((double)nBufYSize + 0.5));
273 1 : tmp.Printf("&SCALEFACTOR=%.15g", MIN(fx, fy));
274 : }
275 : else
276 : {
277 : std::vector<std::string> grid_axes =
278 12 : Split(CPLGetXMLValue(psService, "GridAxes", ""), ",");
279 6 : if (grid_axes.size() < 2)
280 : {
281 : // eek!
282 0 : grid_axes.push_back("E");
283 0 : grid_axes.push_back("N");
284 : }
285 : tmp.Printf("&SCALESIZE=%s%%28%i%%29,%s%%28%i%%29",
286 12 : grid_axes[0].c_str(), nBufXSize, grid_axes[1].c_str(),
287 12 : nBufYSize);
288 : }
289 7 : request += tmp;
290 : }
291 :
292 21 : if (range != "" && range != "*")
293 : {
294 0 : request += "&RANGESUBSET=" + range;
295 : }
296 :
297 : // other parameters may come from
298 : // 1) URL (others)
299 : // 2) Service file
300 21 : const char *keys[] = {WCS_URL_PARAMETERS};
301 231 : for (unsigned int i = 0; i < CPL_ARRAYSIZE(keys); i++)
302 : {
303 420 : std::string value;
304 210 : int ix = IndexOf(CPLString(keys[i]).toupper(), others);
305 210 : if (ix >= 0)
306 : {
307 0 : value = others[ix][1];
308 : }
309 : else
310 : {
311 210 : value = CPLGetXMLValue(psService, keys[i], "");
312 : }
313 210 : if (value != "")
314 : {
315 21 : request = CPLURLAddKVP(request.c_str(), keys[i], value.c_str());
316 : }
317 : }
318 : // add extra parameters
319 42 : std::string extra = CPLGetXMLValue(psService, "Parameters", "");
320 21 : if (extra != "")
321 : {
322 42 : std::vector<std::string> pairs = Split(extra.c_str(), "&");
323 42 : for (unsigned int i = 0; i < pairs.size(); ++i)
324 : {
325 21 : std::vector<std::string> pair = Split(pairs[i].c_str(), "=");
326 : request =
327 21 : CPLURLAddKVP(request.c_str(), pair[0].c_str(), pair[1].c_str());
328 : }
329 : }
330 : std::vector<std::string> pairs =
331 42 : Split(CPLGetXMLValue(psService, "GetCoverageExtra", ""), "&");
332 42 : for (unsigned int i = 0; i < pairs.size(); ++i)
333 : {
334 42 : std::vector<std::string> pair = Split(pairs[i].c_str(), "=");
335 21 : if (pair.size() > 1)
336 : {
337 : request =
338 21 : CPLURLAddKVP(request.c_str(), pair[0].c_str(), pair[1].c_str());
339 : }
340 : }
341 :
342 21 : CPLDebug("WCS", "Requesting %s", request.c_str());
343 42 : return request;
344 : }
345 :
346 : /************************************************************************/
347 : /* DescribeCoverageRequest() */
348 : /* */
349 : /************************************************************************/
350 :
351 7 : std::string WCSDataset201::DescribeCoverageRequest()
352 : {
353 7 : std::string request = CPLGetXMLValue(psService, "ServiceURL", "");
354 7 : request = CPLURLAddKVP(request.c_str(), "SERVICE", "WCS");
355 7 : request = CPLURLAddKVP(request.c_str(), "REQUEST", "DescribeCoverage");
356 14 : request = CPLURLAddKVP(request.c_str(), "VERSION",
357 14 : CPLGetXMLValue(psService, "Version", "2.0.1"));
358 14 : request = CPLURLAddKVP(request.c_str(), "COVERAGEID",
359 14 : CPLGetXMLValue(psService, "CoverageName", ""));
360 14 : std::string extra = CPLGetXMLValue(psService, "Parameters", "");
361 7 : if (extra != "")
362 : {
363 14 : std::vector<std::string> pairs = Split(extra.c_str(), "&");
364 14 : for (unsigned int i = 0; i < pairs.size(); ++i)
365 : {
366 7 : std::vector<std::string> pair = Split(pairs[i].c_str(), "=");
367 : request =
368 7 : CPLURLAddKVP(request.c_str(), pair[0].c_str(), pair[1].c_str());
369 : }
370 : }
371 7 : extra = CPLGetXMLValue(psService, "DescribeCoverageExtra", "");
372 7 : if (extra != "")
373 : {
374 0 : std::vector<std::string> pairs = Split(extra.c_str(), "&");
375 0 : for (unsigned int i = 0; i < pairs.size(); ++i)
376 : {
377 0 : std::vector<std::string> pair = Split(pairs[i].c_str(), "=");
378 : request =
379 0 : CPLURLAddKVP(request.c_str(), pair[0].c_str(), pair[1].c_str());
380 : }
381 : }
382 7 : CPLDebug("WCS", "Requesting %s", request.c_str());
383 14 : return request;
384 : }
385 :
386 : /************************************************************************/
387 : /* GridOffsets() */
388 : /* */
389 : /************************************************************************/
390 :
391 21 : bool WCSDataset201::GridOffsets(CPLXMLNode *grid, const std::string &subtype,
392 : bool swap_grid_axis,
393 : std::vector<double> &origin,
394 : std::vector<std::vector<double>> &offset,
395 : std::vector<std::string> axes, char ***metadata)
396 : {
397 : // todo: use domain_index
398 :
399 : // origin position, center of cell
400 21 : CPLXMLNode *point = CPLGetXMLNode(grid, "origin.Point.pos");
401 42 : origin = Flist(
402 63 : Split(CPLGetXMLValue(point, nullptr, ""), " ", axis_order_swap), 0, 2);
403 :
404 : // offsets = coefficients of affine transformation from cell coords to
405 : // CRS coords, (1,2) and (4,5)
406 :
407 21 : if (subtype == "RectifiedGrid")
408 : {
409 :
410 : // for rectified grid the geo transform is from origin and offsetVectors
411 18 : int i = 0;
412 126 : for (CPLXMLNode *node = grid->psChild; node != nullptr;
413 108 : node = node->psNext)
414 : {
415 126 : if (node->eType != CXT_Element ||
416 90 : !EQUAL(node->pszValue, "offsetVector"))
417 : {
418 90 : continue;
419 : }
420 36 : offset.push_back(Flist(
421 72 : Split(CPLGetXMLValue(node, nullptr, ""), " ", axis_order_swap),
422 : 0, 2));
423 36 : if (i == 1)
424 : {
425 18 : break;
426 : }
427 18 : i++;
428 : }
429 18 : if (offset.size() < 2)
430 : {
431 : // error or not?
432 0 : offset.push_back(std::vector<double>{1, 0}); // x
433 0 : offset.push_back(std::vector<double>{0, 1}); // y
434 : }
435 : // if axis_order_swap
436 : // the offset order should be swapped
437 : // Rasdaman does it
438 : // MapServer and GeoServer not
439 18 : if (swap_grid_axis)
440 : {
441 3 : std::swap(offset[0], offset[1]);
442 : }
443 : }
444 : else
445 : { // if (coverage_type == "ReferenceableGrid"(ByVector)) {
446 :
447 : // for vector referenceable grid the geo transform is from
448 : // offsetVector, coefficients, gridAxesSpanned, sequenceRule
449 : // in generalGridAxis.GeneralGridAxis
450 33 : for (CPLXMLNode *node = grid->psChild; node != nullptr;
451 30 : node = node->psNext)
452 : {
453 30 : CPLXMLNode *axis = CPLGetXMLNode(node, "GeneralGridAxis");
454 30 : if (!axis)
455 : {
456 21 : continue;
457 : }
458 9 : std::string spanned = CPLGetXMLValue(axis, "gridAxesSpanned", "");
459 9 : int index = IndexOf(spanned, axes);
460 9 : if (index == -1)
461 : {
462 0 : CPLError(CE_Failure, CPLE_AppDefined,
463 : "This is not a rectilinear grid(?).");
464 0 : return false;
465 : }
466 9 : std::string coeffs = CPLGetXMLValue(axis, "coefficients", "");
467 9 : if (coeffs != "")
468 : {
469 6 : *metadata = CSLSetNameValue(
470 : *metadata,
471 6 : CPLString().Printf("DIMENSION_%i_COEFFS", index).c_str(),
472 : coeffs.c_str());
473 : }
474 : std::string order =
475 9 : CPLGetXMLValue(axis, "sequenceRule.axisOrder", "");
476 9 : std::string rule = CPLGetXMLValue(axis, "sequenceRule", "");
477 9 : if (!(order == "+1" && rule == "Linear"))
478 : {
479 0 : CPLError(CE_Failure, CPLE_AppDefined,
480 : "Grids with sequence rule '%s' and axis order '%s' "
481 : "are not supported.",
482 : rule.c_str(), order.c_str());
483 0 : return false;
484 : }
485 9 : CPLXMLNode *offset_node = CPLGetXMLNode(axis, "offsetVector");
486 9 : if (offset_node)
487 : {
488 9 : offset.push_back(
489 18 : Flist(Split(CPLGetXMLValue(offset_node, nullptr, ""), " ",
490 9 : axis_order_swap),
491 : 0, 2));
492 : }
493 : else
494 : {
495 0 : CPLError(CE_Failure, CPLE_AppDefined,
496 : "Missing offset vector in grid axis.");
497 0 : return false;
498 : }
499 : }
500 : // todo: make sure offset order is the same as the axes order but see
501 : // above
502 : }
503 21 : if (origin.size() < 2 || offset.size() < 2)
504 : {
505 0 : CPLError(CE_Failure, CPLE_AppDefined,
506 : "Could not parse origin or offset vectors from grid.");
507 0 : return false;
508 : }
509 21 : return true;
510 : }
511 :
512 : /************************************************************************/
513 : /* GetSubdataset() */
514 : /* */
515 : /************************************************************************/
516 :
517 0 : std::string WCSDataset201::GetSubdataset(const std::string &coverage)
518 : {
519 0 : char **metadata = GDALPamDataset::GetMetadata("SUBDATASETS");
520 0 : std::string subdataset;
521 0 : if (metadata != nullptr)
522 : {
523 0 : for (int i = 0; metadata[i] != nullptr; ++i)
524 : {
525 : char *key;
526 0 : std::string url = CPLParseNameValue(metadata[i], &key);
527 0 : if (key != nullptr && strstr(key, "SUBDATASET_") &&
528 0 : strstr(key, "_NAME"))
529 : {
530 0 : if (coverage == CPLURLGetValue(url.c_str(), "coverageId"))
531 : {
532 0 : subdataset = key;
533 0 : subdataset.erase(subdataset.find("_NAME"), 5);
534 0 : CPLFree(key);
535 0 : break;
536 : }
537 : }
538 0 : CPLFree(key);
539 : }
540 : }
541 0 : return subdataset;
542 : }
543 :
544 : /************************************************************************/
545 : /* SetFormat() */
546 : /* */
547 : /************************************************************************/
548 :
549 21 : bool WCSDataset201::SetFormat(CPLXMLNode *coverage)
550 : {
551 : // set the Format value in service,
552 : // unless it is set by the user
553 42 : std::string format = CPLGetXMLValue(psService, "Format", "");
554 :
555 : // todo: check the value against list of supported formats?
556 21 : if (format != "")
557 : {
558 14 : return true;
559 : }
560 :
561 : /* We will prefer anything that sounds like TIFF, otherwise */
562 : /* falling back to the first supported format. Should we */
563 : /* consider preferring the nativeFormat if available? */
564 :
565 7 : char **metadata = GDALPamDataset::GetMetadata(nullptr);
566 : const char *value =
567 7 : CSLFetchNameValue(metadata, "WCS_GLOBAL#formatSupported");
568 7 : if (value == nullptr)
569 : {
570 0 : format = CPLGetXMLValue(coverage, "ServiceParameters.nativeFormat", "");
571 : }
572 : else
573 : {
574 14 : std::vector<std::string> format_list = Split(value, ",");
575 32 : for (unsigned j = 0; j < format_list.size(); ++j)
576 : {
577 32 : if (CPLString(format_list[j]).ifind("tiff") != std::string::npos)
578 : {
579 7 : format = format_list[j];
580 7 : break;
581 : }
582 : }
583 7 : if (format == "" && format_list.size() > 0)
584 : {
585 0 : format = format_list[0];
586 : }
587 : }
588 7 : if (format != "")
589 : {
590 7 : CPLSetXMLValue(psService, "Format", format.c_str());
591 7 : bServiceDirty = true;
592 7 : return true;
593 : }
594 : else
595 : {
596 0 : return false;
597 : }
598 : }
599 :
600 : /************************************************************************/
601 : /* ParseGridFunction() */
602 : /* */
603 : /************************************************************************/
604 :
605 21 : bool WCSDataset201::ParseGridFunction(CPLXMLNode *coverage,
606 : std::vector<int> &axisOrder)
607 : {
608 : CPLXMLNode *function =
609 21 : CPLGetXMLNode(coverage, "coverageFunction.GridFunction");
610 21 : if (function)
611 : {
612 15 : std::string path = "sequenceRule";
613 15 : std::string sequenceRule = CPLGetXMLValue(function, path.c_str(), "");
614 15 : path += ".axisOrder";
615 : axisOrder =
616 15 : Ilist(Split(CPLGetXMLValue(function, path.c_str(), ""), " "));
617 : // for now require simple
618 15 : if (sequenceRule != "Linear")
619 : {
620 0 : CPLError(CE_Failure, CPLE_AppDefined,
621 : "Can't handle '%s' coverages.", sequenceRule.c_str());
622 0 : return false;
623 : }
624 : }
625 21 : return true;
626 : }
627 :
628 : /************************************************************************/
629 : /* ParseRange() */
630 : /* */
631 : /************************************************************************/
632 :
633 21 : int WCSDataset201::ParseRange(CPLXMLNode *coverage,
634 : const std::string &range_subset, char ***metadata)
635 : {
636 21 : int fields = 0;
637 : // Default is to include all (types permitting?)
638 : // Can also be controlled with Range parameter
639 :
640 : // The contents of a rangeType is a swe:DataRecord
641 21 : const char *path = "rangeType.DataRecord";
642 21 : CPLXMLNode *record = CPLGetXMLNode(coverage, path);
643 21 : if (!record)
644 : {
645 0 : CPLError(CE_Failure, CPLE_AppDefined,
646 : "Attributes are not defined in a DataRecord, giving up.");
647 0 : return 0;
648 : }
649 :
650 : // mapserver does not like field names, it wants indexes
651 : // so we should be able to give those
652 :
653 : // if Range is set remove those not in it
654 42 : std::vector<std::string> range = Split(range_subset.c_str(), ",");
655 : // todo: add check for range subsetting profile existence in server metadata
656 : // here
657 21 : unsigned int range_index = 0; // index for reading from range
658 21 : bool in_band_range = false;
659 :
660 21 : unsigned int field_index = 1;
661 21 : std::vector<std::string> nodata_array;
662 :
663 63 : for (CPLXMLNode *field = record->psChild; field != nullptr;
664 42 : field = field->psNext)
665 : {
666 42 : if (field->eType != CXT_Element || !EQUAL(field->pszValue, "field"))
667 : {
668 0 : continue;
669 : }
670 42 : std::string fname = CPLGetXMLValue(field, "name", "");
671 42 : bool include = true;
672 :
673 42 : if (range.size() > 0)
674 : {
675 0 : include = false;
676 0 : if (range_index < range.size())
677 : {
678 0 : std::string current_range = range[range_index];
679 0 : std::string fname_test;
680 :
681 0 : if (atoi(current_range.c_str()) != 0)
682 : {
683 0 : fname_test = CPLString().Printf("%i", field_index);
684 : }
685 : else
686 : {
687 0 : fname_test = fname;
688 : }
689 :
690 0 : if (current_range == "*")
691 : {
692 0 : include = true;
693 : }
694 0 : else if (current_range == fname_test)
695 : {
696 0 : include = true;
697 0 : range_index += 1;
698 : }
699 0 : else if (current_range.find(fname_test + ":") !=
700 : std::string::npos)
701 : {
702 0 : include = true;
703 0 : in_band_range = true;
704 : }
705 0 : else if (current_range.find(":" + fname_test) !=
706 : std::string::npos)
707 : {
708 0 : include = true;
709 0 : in_band_range = false;
710 0 : range_index += 1;
711 : }
712 0 : else if (in_band_range)
713 : {
714 0 : include = true;
715 : }
716 : }
717 : }
718 :
719 42 : if (include)
720 : {
721 : const std::string key =
722 84 : CPLString().Printf("FIELD_%i_", field_index);
723 42 : *metadata = CSLSetNameValue(*metadata, (key + "NAME").c_str(),
724 : fname.c_str());
725 :
726 : std::string nodata =
727 84 : CPLGetXMLValue(field, "Quantity.nilValues.NilValue", "");
728 42 : if (nodata != "")
729 : {
730 0 : *metadata = CSLSetNameValue(*metadata, (key + "NODATA").c_str(),
731 : nodata.c_str());
732 : }
733 :
734 : std::string descr =
735 84 : CPLGetXMLValue(field, "Quantity.description", "");
736 42 : if (descr != "")
737 : {
738 18 : *metadata = CSLSetNameValue(*metadata, (key + "DESCR").c_str(),
739 : descr.c_str());
740 : }
741 :
742 42 : path = "Quantity.constraint.AllowedValues.interval";
743 42 : std::string interval = CPLGetXMLValue(field, path, "");
744 42 : if (interval != "")
745 : {
746 54 : *metadata = CSLSetNameValue(
747 54 : *metadata, (key + "INTERVAL").c_str(), interval.c_str());
748 : }
749 :
750 42 : nodata_array.push_back(std::move(nodata));
751 42 : fields += 1;
752 : }
753 :
754 42 : field_index += 1;
755 : }
756 :
757 21 : if (fields == 0)
758 : {
759 0 : CPLError(CE_Failure, CPLE_AppDefined,
760 : "No data fields found (bad Range?).");
761 : }
762 : else
763 : {
764 : // todo: default to the first one?
765 21 : bServiceDirty = CPLUpdateXML(psService, "NoDataValue",
766 59 : Join(nodata_array, ",").c_str()) ||
767 17 : bServiceDirty;
768 : }
769 :
770 21 : return fields;
771 : }
772 :
773 : /************************************************************************/
774 : /* ExtractGridInfo() */
775 : /* */
776 : /* Collect info about grid from describe coverage for WCS 2.0. */
777 : /* */
778 : /************************************************************************/
779 :
780 21 : bool WCSDataset201::ExtractGridInfo()
781 : {
782 : // this is for checking what's in service and for filling in empty slots in
783 : // it if the service file can be considered ready for use, this could be
784 : // skipped
785 :
786 21 : CPLXMLNode *coverage = CPLGetXMLNode(psService, "CoverageDescription");
787 :
788 21 : if (coverage == nullptr)
789 : {
790 0 : CPLError(CE_Failure, CPLE_AppDefined,
791 : "CoverageDescription missing from service.");
792 0 : return false;
793 : }
794 :
795 42 : std::string subtype = CoverageSubtype(coverage);
796 :
797 : // get CRS from boundedBy.Envelope and set the native flag to true
798 : // below we may set the CRS again but that won't be native (however, non
799 : // native CRS is not yet supported) also axis order swap is set
800 42 : std::string path = "boundedBy.Envelope";
801 21 : CPLXMLNode *envelope = CPLGetXMLNode(coverage, path.c_str());
802 21 : if (envelope == nullptr)
803 : {
804 0 : path = "boundedBy.EnvelopeWithTimePeriod";
805 0 : envelope = CPLGetXMLNode(coverage, path.c_str());
806 0 : if (envelope == nullptr)
807 : {
808 0 : CPLError(CE_Failure, CPLE_AppDefined, "Missing boundedBy.Envelope");
809 0 : return false;
810 : }
811 : }
812 42 : std::vector<std::string> bbox = ParseBoundingBox(envelope);
813 21 : if (!SetCRS(ParseCRS(envelope), true) || bbox.size() < 2)
814 : {
815 0 : return false;
816 : }
817 :
818 : // has the user set the domain?
819 : std::vector<std::string> domain =
820 42 : Split(CPLGetXMLValue(psService, "Domain", ""), ",");
821 :
822 : // names of axes
823 : std::vector<std::string> axes =
824 0 : Split(CPLGetXMLValue(coverage, (path + ".axisLabels").c_str(), ""), " ",
825 42 : axis_order_swap);
826 : std::vector<std::string> uoms =
827 0 : Split(CPLGetXMLValue(coverage, (path + ".uomLabels").c_str(), ""), " ",
828 42 : axis_order_swap);
829 :
830 21 : if (axes.size() < 2)
831 : {
832 0 : CPLError(CE_Failure, CPLE_AppDefined,
833 : "The coverage has less than 2 dimensions or no axisLabels.");
834 0 : return false;
835 : }
836 :
837 42 : std::vector<int> domain_indexes = IndexOf(domain, axes);
838 21 : if (Contains(domain_indexes, -1))
839 : {
840 0 : CPLError(CE_Failure, CPLE_AppDefined,
841 : "Axis in given domain does not exist in coverage.");
842 0 : return false;
843 : }
844 21 : if (domain_indexes.size() == 0)
845 : { // default is the first two
846 7 : domain_indexes.push_back(0);
847 7 : domain_indexes.push_back(1);
848 : }
849 21 : if (domain.size() == 0)
850 : {
851 7 : domain.push_back(axes[0]);
852 7 : domain.push_back(axes[1]);
853 7 : CPLSetXMLValue(psService, "Domain", Join(domain, ",").c_str());
854 7 : bServiceDirty = true;
855 : }
856 :
857 : // GridFunction (is optional)
858 : // We support only linear grid functions.
859 : // axisOrder determines how data is arranged in the grid <order><axis
860 : // number> specifically: +2 +1 => swap grid envelope and the order of the
861 : // offsets
862 42 : std::vector<int> axisOrder;
863 21 : if (!ParseGridFunction(coverage, axisOrder))
864 : {
865 0 : return false;
866 : }
867 :
868 21 : const char *md_domain = "";
869 63 : char **metadata = CSLDuplicate(
870 21 : GetMetadata(md_domain)); // coverage metadata to be added/updated
871 :
872 21 : metadata = CSLSetNameValue(metadata, "DOMAIN", Join(domain, ",").c_str());
873 :
874 : // add coverage metadata: GeoServer TimeDomain
875 :
876 : CPLXMLNode *timedomain =
877 21 : CPLGetXMLNode(coverage, "metadata.Extension.TimeDomain");
878 21 : if (timedomain)
879 : {
880 0 : std::vector<std::string> timePositions;
881 : // "//timePosition"
882 0 : for (CPLXMLNode *node = timedomain->psChild; node != nullptr;
883 0 : node = node->psNext)
884 : {
885 0 : if (node->eType != CXT_Element ||
886 0 : strcmp(node->pszValue, "TimeInstant") != 0)
887 : {
888 0 : continue;
889 : }
890 0 : for (CPLXMLNode *node2 = node->psChild; node2 != nullptr;
891 0 : node2 = node2->psNext)
892 : {
893 0 : if (node2->eType != CXT_Element ||
894 0 : strcmp(node2->pszValue, "timePosition") != 0)
895 : {
896 0 : continue;
897 : }
898 0 : timePositions.push_back(CPLGetXMLValue(node2, "", ""));
899 : }
900 : }
901 0 : metadata = CSLSetNameValue(metadata, "TimeDomain",
902 0 : Join(timePositions, ",").c_str());
903 : }
904 :
905 : // dimension metadata
906 :
907 : std::vector<std::string> slow =
908 42 : Split(bbox[0].c_str(), " ", axis_order_swap);
909 : std::vector<std::string> shigh =
910 42 : Split(bbox[1].c_str(), " ", axis_order_swap);
911 35 : bServiceDirty = CPLUpdateXML(psService, "Low", Join(slow, ",").c_str()) ||
912 14 : bServiceDirty;
913 35 : bServiceDirty = CPLUpdateXML(psService, "High", Join(shigh, ",").c_str()) ||
914 14 : bServiceDirty;
915 21 : if (slow.size() < 2 || shigh.size() < 2)
916 : {
917 0 : CPLError(CE_Failure, CPLE_AppDefined,
918 : "The coverage has less than 2 dimensions.");
919 0 : CSLDestroy(metadata);
920 0 : return false;
921 : }
922 : // todo: if our x,y domain is not the first two? use domain_indexes?
923 42 : std::vector<double> low = Flist(slow, 0, 2);
924 42 : std::vector<double> high = Flist(shigh, 0, 2);
925 42 : std::vector<double> env;
926 21 : env.insert(env.end(), low.begin(), low.begin() + 2);
927 21 : env.insert(env.end(), high.begin(), high.begin() + 2);
928 :
929 66 : for (unsigned int i = 0; i < axes.size(); ++i)
930 : {
931 90 : const std::string key = CPLString().Printf("DIMENSION_%i_", i);
932 45 : metadata =
933 45 : CSLSetNameValue(metadata, (key + "AXIS").c_str(), axes[i].c_str());
934 45 : if (i < uoms.size())
935 : {
936 39 : metadata = CSLSetNameValue(metadata, (key + "UOM").c_str(),
937 39 : uoms[i].c_str());
938 : }
939 45 : if (i < 2)
940 : {
941 42 : metadata = CSLSetNameValue(
942 84 : metadata, (key + "INTERVAL").c_str(),
943 84 : CPLString().Printf("%.15g,%.15g", low[i], high[i]));
944 : }
945 3 : else if (i < slow.size() && i < shigh.size())
946 : {
947 3 : metadata = CSLSetNameValue(
948 6 : metadata, (key + "INTERVAL").c_str(),
949 6 : CPLString().Printf("%s,%s", slow[i].c_str(), shigh[i].c_str()));
950 : }
951 0 : else if (i < bbox.size())
952 : {
953 0 : metadata = CSLSetNameValue(metadata, (key + "INTERVAL").c_str(),
954 0 : bbox[i].c_str());
955 : }
956 : }
957 :
958 : // domainSet
959 : // requirement 23: the srsName here _shall_ be the same as in boundedBy
960 : // => we ignore it
961 : // the CRS of this dataset is from boundedBy (unless it is overridden)
962 : // this is the size of this dataset
963 : // this gives the geotransform of this dataset (unless there is CRS
964 : // override)
965 :
966 21 : CPLXMLNode *grid = GetGridNode(coverage, subtype);
967 21 : if (!grid)
968 : {
969 0 : CSLDestroy(metadata);
970 0 : return false;
971 : }
972 :
973 : //
974 21 : bool swap_grid_axis = false;
975 30 : if (axisOrder.size() >= 2 && axisOrder[domain_indexes[0]] == 2 &&
976 9 : axisOrder[domain_indexes[1]] == 1)
977 : {
978 9 : swap_grid_axis = !CPLGetXMLBoolean(psService, "NoGridAxisSwap");
979 : }
980 21 : path = "limits.GridEnvelope";
981 : std::vector<std::vector<int>> size =
982 42 : ParseGridEnvelope(CPLGetXMLNode(grid, path.c_str()), swap_grid_axis);
983 42 : std::vector<int> grid_size;
984 21 : if (size.size() < 2)
985 : {
986 0 : CPLError(CE_Failure, CPLE_AppDefined, "Can't parse the grid envelope.");
987 0 : CSLDestroy(metadata);
988 0 : return false;
989 : }
990 :
991 21 : grid_size.push_back(size[1][domain_indexes[0]] -
992 21 : size[0][domain_indexes[0]] + 1);
993 21 : grid_size.push_back(size[1][domain_indexes[1]] -
994 21 : size[0][domain_indexes[1]] + 1);
995 :
996 21 : path = "axisLabels";
997 : bool swap_grid_axis_labels =
998 21 : swap_grid_axis || CPLGetXMLBoolean(psService, "GridAxisLabelSwap");
999 : std::vector<std::string> grid_axes = Split(
1000 42 : CPLGetXMLValue(grid, path.c_str(), ""), " ", swap_grid_axis_labels);
1001 : // autocorrect MapServer thing
1002 24 : if (grid_axes.size() >= 2 && grid_axes[0] == "lat" &&
1003 3 : grid_axes[1] == "long")
1004 : {
1005 3 : grid_axes[0] = "long";
1006 3 : grid_axes[1] = "lat";
1007 : }
1008 21 : bServiceDirty =
1009 35 : CPLUpdateXML(psService, "GridAxes", Join(grid_axes, ",").c_str()) ||
1010 14 : bServiceDirty;
1011 :
1012 42 : std::vector<double> origin;
1013 42 : std::vector<std::vector<double>> offsets;
1014 21 : if (!GridOffsets(grid, subtype, swap_grid_axis, origin, offsets, axes,
1015 : &metadata))
1016 : {
1017 0 : CSLDestroy(metadata);
1018 0 : return false;
1019 : }
1020 :
1021 21 : SetGeometry(grid_size, origin, offsets);
1022 :
1023 : // subsetting and dimension to bands
1024 42 : std::vector<std::string> dimensions;
1025 42 : std::string range;
1026 42 : std::vector<std::vector<std::string>> others;
1027 21 : ParseParameters(psService, dimensions, range, others);
1028 :
1029 : // it is ok to have trimming or even slicing for x/y, it just affects our
1030 : // bounding box but that is a todo item todo: BoundGeometry(domain_trim) if
1031 : // domain_trim.size() > 0
1032 42 : std::vector<std::vector<double>> domain_trim;
1033 :
1034 : // are all dimensions that are not x/y domain sliced?
1035 : // if not, bands can't be defined, see below
1036 21 : bool dimensions_are_ok = true;
1037 66 : for (unsigned int i = 0; i < axes.size(); ++i)
1038 : {
1039 45 : std::vector<std::string> params;
1040 51 : for (unsigned int j = 0; j < dimensions.size(); ++j)
1041 : {
1042 9 : if (dimensions[j].find(axes[i] + "(") != std::string::npos)
1043 : {
1044 3 : params = Split(FromParenthesis(dimensions[j]).c_str(), ",");
1045 3 : break;
1046 : }
1047 : }
1048 45 : int domain_index = IndexOf(axes[i], domain);
1049 45 : if (domain_index != -1)
1050 : {
1051 42 : domain_trim.push_back(Flist(params, 0, 2));
1052 42 : continue;
1053 : }
1054 : // size == 1 => sliced
1055 3 : if (params.size() != 1)
1056 : {
1057 0 : dimensions_are_ok = false;
1058 : }
1059 : }
1060 : // todo: add metadata: note: no bands, you need to subset to get data
1061 :
1062 : // check for CRS override
1063 42 : std::string crs = CPLGetXMLValue(psService, "SRS", "");
1064 21 : if (crs != "" && crs != osCRS)
1065 : {
1066 0 : if (!SetCRS(crs, false))
1067 : {
1068 0 : CSLDestroy(metadata);
1069 0 : return false;
1070 : }
1071 : // todo: support CRS override, it requires warping the grid to the new
1072 : // CRS SetGeometry(grid_size, origin, offsets);
1073 0 : CPLError(CE_Failure, CPLE_AppDefined,
1074 : "CRS override not yet supported.");
1075 0 : CSLDestroy(metadata);
1076 0 : return false;
1077 : }
1078 :
1079 : // todo: ElevationDomain, DimensionDomain
1080 :
1081 : // rangeType
1082 :
1083 : // get the field metadata
1084 : // get the count of fields
1085 : // if Range is set in service that may limit the fields
1086 21 : int fields = ParseRange(coverage, range, &metadata);
1087 : // if fields is 0 an error message has been emitted
1088 : // but we let this go on since the user may be experimenting
1089 : // and she wants to see the resulting metadata and not just an error message
1090 : // situation is ~the same when bands == 0 when we exit here
1091 :
1092 : // todo: do this only if metadata is dirty
1093 21 : this->SetMetadata(metadata, md_domain);
1094 21 : CSLDestroy(metadata);
1095 21 : TrySaveXML();
1096 :
1097 : // determine the band count
1098 21 : int bands = 0;
1099 21 : if (dimensions_are_ok)
1100 : {
1101 21 : bands = fields;
1102 : }
1103 21 : bServiceDirty =
1104 35 : CPLUpdateXML(psService, "BandCount", CPLString().Printf("%d", bands)) ||
1105 14 : bServiceDirty;
1106 :
1107 : // set the Format value in service, unless it is set
1108 : // by the user, either through direct edit or options
1109 21 : if (!SetFormat(coverage))
1110 : {
1111 : // all attempts to find a format have failed...
1112 0 : CPLError(CE_Failure, CPLE_AppDefined,
1113 : "All attempts to find a format have failed, giving up.");
1114 0 : return false;
1115 : }
1116 :
1117 21 : return true;
1118 : }
|
