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