Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: PDS 4 Driver; Planetary Data System Format
4 : * Purpose: Implementation of PDS4Dataset
5 : * Author: Even Rouault, even.rouault at spatialys.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2017, Hobu Inc
9 : *
10 : * SPDX-License-Identifier: MIT
11 : ****************************************************************************/
12 :
13 : #include "cpl_vsi_error.h"
14 : #include "gdal_proxy.h"
15 : #include "gdal_frmts.h"
16 : #include "rawdataset.h"
17 : #include "vrtdataset.h"
18 : #include "ogrsf_frmts.h"
19 : #include "ogr_spatialref.h"
20 : #include "gdal_priv_templates.hpp"
21 : #include "ogreditablelayer.h"
22 : #include "pds4dataset.h"
23 : #include "pdsdrivercore.h"
24 :
25 : #ifdef EMBED_RESOURCE_FILES
26 : #include "embedded_resources.h"
27 : #endif
28 :
29 : #include <cstdlib>
30 : #include <vector>
31 : #include <algorithm>
32 : #include <optional>
33 :
34 : #define TIFF_GEOTIFF_STRING "TIFF 6.0"
35 : #define BIGTIFF_GEOTIFF_STRING "TIFF 6.0"
36 : #define PREEXISTING_BINARY_FILE \
37 : "Binary file pre-existing PDS4 label. This comment is used by GDAL to " \
38 : "avoid deleting the binary file when the label is deleted. Keep it to " \
39 : "preserve this behavior."
40 :
41 : #define CURRENT_CART_VERSION "1O00_1970"
42 :
43 : /************************************************************************/
44 : /* PDS4BrowseImageProxyRasterBand() */
45 : /************************************************************************/
46 :
47 2 : PDS4BrowseImageProxyRasterBand::PDS4BrowseImageProxyRasterBand(
48 2 : GDALRasterBand *poBaseBandIn)
49 2 : : m_poBaseBand(poBaseBandIn)
50 : {
51 2 : eDataType = m_poBaseBand->GetRasterDataType();
52 2 : m_poBaseBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
53 2 : }
54 :
55 : /************************************************************************/
56 : /* RefUnderlyingRasterBand() */
57 : /************************************************************************/
58 :
59 3 : GDALRasterBand *PDS4BrowseImageProxyRasterBand::RefUnderlyingRasterBand(
60 : bool /*bForceOpen*/) const
61 : {
62 3 : return m_poBaseBand;
63 : }
64 :
65 : /************************************************************************/
66 : /* PDS4WrapperRasterBand() */
67 : /************************************************************************/
68 :
69 24 : PDS4WrapperRasterBand::PDS4WrapperRasterBand(GDALRasterBand *poBaseBandIn)
70 24 : : m_poBaseBand(poBaseBandIn)
71 : {
72 24 : eDataType = m_poBaseBand->GetRasterDataType();
73 24 : m_poBaseBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
74 24 : }
75 :
76 : /************************************************************************/
77 : /* SetMaskBand() */
78 : /************************************************************************/
79 :
80 0 : void PDS4WrapperRasterBand::SetMaskBand(
81 : std::unique_ptr<GDALRasterBand> poMaskBand)
82 : {
83 0 : poMask.reset(std::move(poMaskBand));
84 0 : nMaskFlags = 0;
85 0 : }
86 :
87 : /************************************************************************/
88 : /* GetOffset() */
89 : /************************************************************************/
90 :
91 18 : double PDS4WrapperRasterBand::GetOffset(int *pbSuccess)
92 : {
93 18 : if (pbSuccess)
94 18 : *pbSuccess = m_bHasOffset;
95 18 : return m_dfOffset;
96 : }
97 :
98 : /************************************************************************/
99 : /* GetScale() */
100 : /************************************************************************/
101 :
102 18 : double PDS4WrapperRasterBand::GetScale(int *pbSuccess)
103 : {
104 18 : if (pbSuccess)
105 18 : *pbSuccess = m_bHasScale;
106 18 : return m_dfScale;
107 : }
108 :
109 : /************************************************************************/
110 : /* SetOffset() */
111 : /************************************************************************/
112 :
113 1 : CPLErr PDS4WrapperRasterBand::SetOffset(double dfNewOffset)
114 : {
115 1 : m_dfOffset = dfNewOffset;
116 1 : m_bHasOffset = true;
117 :
118 1 : PDS4Dataset *poGDS = cpl::down_cast<PDS4Dataset *>(poDS);
119 1 : if (poGDS->m_poExternalDS && eAccess == GA_Update)
120 1 : poGDS->m_poExternalDS->GetRasterBand(nBand)->SetOffset(dfNewOffset);
121 :
122 1 : return CE_None;
123 : }
124 :
125 : /************************************************************************/
126 : /* SetScale() */
127 : /************************************************************************/
128 :
129 1 : CPLErr PDS4WrapperRasterBand::SetScale(double dfNewScale)
130 : {
131 1 : m_dfScale = dfNewScale;
132 1 : m_bHasScale = true;
133 :
134 1 : PDS4Dataset *poGDS = cpl::down_cast<PDS4Dataset *>(poDS);
135 1 : if (poGDS->m_poExternalDS && eAccess == GA_Update)
136 1 : poGDS->m_poExternalDS->GetRasterBand(nBand)->SetScale(dfNewScale);
137 :
138 1 : return CE_None;
139 : }
140 :
141 : /************************************************************************/
142 : /* GetNoDataValue() */
143 : /************************************************************************/
144 :
145 35 : double PDS4WrapperRasterBand::GetNoDataValue(int *pbSuccess)
146 : {
147 35 : if (m_bHasNoDataInt64)
148 : {
149 0 : if (pbSuccess)
150 0 : *pbSuccess = true;
151 0 : return GDALGetNoDataValueCastToDouble(m_nNoDataInt64);
152 : }
153 :
154 35 : if (m_bHasNoDataUInt64)
155 : {
156 0 : if (pbSuccess)
157 0 : *pbSuccess = true;
158 0 : return GDALGetNoDataValueCastToDouble(m_nNoDataUInt64);
159 : }
160 :
161 35 : if (pbSuccess)
162 35 : *pbSuccess = m_bHasNoData;
163 35 : return m_dfNoData;
164 : }
165 :
166 : /************************************************************************/
167 : /* SetNoDataValue() */
168 : /************************************************************************/
169 :
170 9 : CPLErr PDS4WrapperRasterBand::SetNoDataValue(double dfNewNoData)
171 : {
172 9 : m_dfNoData = dfNewNoData;
173 9 : m_bHasNoData = true;
174 :
175 9 : PDS4Dataset *poGDS = cpl::down_cast<PDS4Dataset *>(poDS);
176 9 : if (poGDS->m_poExternalDS && eAccess == GA_Update)
177 9 : poGDS->m_poExternalDS->GetRasterBand(nBand)->SetNoDataValue(
178 9 : dfNewNoData);
179 :
180 9 : return CE_None;
181 : }
182 :
183 : /************************************************************************/
184 : /* GetNoDataValueAsInt64() */
185 : /************************************************************************/
186 :
187 3 : int64_t PDS4WrapperRasterBand::GetNoDataValueAsInt64(int *pbSuccess)
188 : {
189 3 : if (pbSuccess)
190 3 : *pbSuccess = m_bHasNoDataInt64;
191 3 : return m_nNoDataInt64;
192 : }
193 :
194 : /************************************************************************/
195 : /* SetNoDataValueAsInt64() */
196 : /************************************************************************/
197 :
198 1 : CPLErr PDS4WrapperRasterBand::SetNoDataValueAsInt64(int64_t nNewNoData)
199 : {
200 1 : m_nNoDataInt64 = nNewNoData;
201 1 : m_bHasNoDataInt64 = true;
202 :
203 1 : PDS4Dataset *poGDS = cpl::down_cast<PDS4Dataset *>(poDS);
204 1 : if (poGDS->m_poExternalDS && eAccess == GA_Update)
205 1 : poGDS->m_poExternalDS->GetRasterBand(nBand)->SetNoDataValueAsInt64(
206 1 : nNewNoData);
207 :
208 1 : return CE_None;
209 : }
210 :
211 : /************************************************************************/
212 : /* GetNoDataValueAsUInt64() */
213 : /************************************************************************/
214 :
215 3 : uint64_t PDS4WrapperRasterBand::GetNoDataValueAsUInt64(int *pbSuccess)
216 : {
217 3 : if (pbSuccess)
218 3 : *pbSuccess = m_bHasNoDataUInt64;
219 3 : return m_nNoDataUInt64;
220 : }
221 :
222 : /************************************************************************/
223 : /* SetNoDataValueAsUInt64() */
224 : /************************************************************************/
225 :
226 1 : CPLErr PDS4WrapperRasterBand::SetNoDataValueAsUInt64(uint64_t nNewNoData)
227 : {
228 1 : m_nNoDataUInt64 = nNewNoData;
229 1 : m_bHasNoDataUInt64 = true;
230 :
231 1 : PDS4Dataset *poGDS = cpl::down_cast<PDS4Dataset *>(poDS);
232 1 : if (poGDS->m_poExternalDS && eAccess == GA_Update)
233 1 : poGDS->m_poExternalDS->GetRasterBand(nBand)->SetNoDataValueAsUInt64(
234 1 : nNewNoData);
235 :
236 1 : return CE_None;
237 : }
238 :
239 : /************************************************************************/
240 : /* Fill() */
241 : /************************************************************************/
242 :
243 2 : CPLErr PDS4WrapperRasterBand::Fill(double dfRealValue, double dfImaginaryValue)
244 : {
245 2 : PDS4Dataset *poGDS = cpl::down_cast<PDS4Dataset *>(poDS);
246 2 : if (poGDS->m_bMustInitImageFile)
247 : {
248 2 : if (!poGDS->InitImageFile())
249 0 : return CE_Failure;
250 : }
251 2 : return GDALProxyRasterBand::Fill(dfRealValue, dfImaginaryValue);
252 : }
253 :
254 : /************************************************************************/
255 : /* IWriteBlock() */
256 : /************************************************************************/
257 :
258 0 : CPLErr PDS4WrapperRasterBand::IWriteBlock(int nXBlock, int nYBlock,
259 : void *pImage)
260 :
261 : {
262 0 : PDS4Dataset *poGDS = cpl::down_cast<PDS4Dataset *>(poDS);
263 0 : if (poGDS->m_bMustInitImageFile)
264 : {
265 0 : if (!poGDS->InitImageFile())
266 0 : return CE_Failure;
267 : }
268 0 : return GDALProxyRasterBand::IWriteBlock(nXBlock, nYBlock, pImage);
269 : }
270 :
271 : /************************************************************************/
272 : /* IRasterIO() */
273 : /************************************************************************/
274 :
275 115 : CPLErr PDS4WrapperRasterBand::IRasterIO(
276 : GDALRWFlag eRWFlag, int nXOff, int nYOff, int nXSize, int nYSize,
277 : void *pData, int nBufXSize, int nBufYSize, GDALDataType eBufType,
278 : GSpacing nPixelSpace, GSpacing nLineSpace, GDALRasterIOExtraArg *psExtraArg)
279 :
280 : {
281 115 : PDS4Dataset *poGDS = cpl::down_cast<PDS4Dataset *>(poDS);
282 115 : if (eRWFlag == GF_Write && poGDS->m_bMustInitImageFile)
283 : {
284 9 : if (!poGDS->InitImageFile())
285 0 : return CE_Failure;
286 : }
287 115 : return GDALProxyRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
288 : pData, nBufXSize, nBufYSize, eBufType,
289 115 : nPixelSpace, nLineSpace, psExtraArg);
290 : }
291 :
292 : /************************************************************************/
293 : /* PDS4RawRasterBand() */
294 : /************************************************************************/
295 :
296 460 : PDS4RawRasterBand::PDS4RawRasterBand(GDALDataset *l_poDS, int l_nBand,
297 : VSILFILE *l_fpRaw,
298 : vsi_l_offset l_nImgOffset,
299 : int l_nPixelOffset, int l_nLineOffset,
300 : GDALDataType l_eDataType,
301 460 : RawRasterBand::ByteOrder eByteOrderIn)
302 : : RawRasterBand(l_poDS, l_nBand, l_fpRaw, l_nImgOffset, l_nPixelOffset,
303 : l_nLineOffset, l_eDataType, eByteOrderIn,
304 460 : RawRasterBand::OwnFP::NO)
305 : {
306 460 : }
307 :
308 : /************************************************************************/
309 : /* SetMaskBand() */
310 : /************************************************************************/
311 :
312 43 : void PDS4RawRasterBand::SetMaskBand(std::unique_ptr<GDALRasterBand> poMaskBand)
313 : {
314 43 : poMask.reset(std::move(poMaskBand));
315 43 : nMaskFlags = 0;
316 43 : }
317 :
318 : /************************************************************************/
319 : /* GetOffset() */
320 : /************************************************************************/
321 :
322 127 : double PDS4RawRasterBand::GetOffset(int *pbSuccess)
323 : {
324 127 : if (pbSuccess)
325 119 : *pbSuccess = m_bHasOffset;
326 127 : return m_dfOffset;
327 : }
328 :
329 : /************************************************************************/
330 : /* GetScale() */
331 : /************************************************************************/
332 :
333 127 : double PDS4RawRasterBand::GetScale(int *pbSuccess)
334 : {
335 127 : if (pbSuccess)
336 119 : *pbSuccess = m_bHasScale;
337 127 : return m_dfScale;
338 : }
339 :
340 : /************************************************************************/
341 : /* SetOffset() */
342 : /************************************************************************/
343 :
344 289 : CPLErr PDS4RawRasterBand::SetOffset(double dfNewOffset)
345 : {
346 289 : m_dfOffset = dfNewOffset;
347 289 : m_bHasOffset = true;
348 289 : return CE_None;
349 : }
350 :
351 : /************************************************************************/
352 : /* SetScale() */
353 : /************************************************************************/
354 :
355 289 : CPLErr PDS4RawRasterBand::SetScale(double dfNewScale)
356 : {
357 289 : m_dfScale = dfNewScale;
358 289 : m_bHasScale = true;
359 289 : return CE_None;
360 : }
361 :
362 : /************************************************************************/
363 : /* GetNoDataValue() */
364 : /************************************************************************/
365 :
366 230 : double PDS4RawRasterBand::GetNoDataValue(int *pbSuccess)
367 : {
368 230 : if (m_bHasNoDataInt64)
369 : {
370 0 : if (pbSuccess)
371 0 : *pbSuccess = true;
372 0 : return GDALGetNoDataValueCastToDouble(m_nNoDataInt64);
373 : }
374 :
375 230 : if (m_bHasNoDataUInt64)
376 : {
377 0 : if (pbSuccess)
378 0 : *pbSuccess = true;
379 0 : return GDALGetNoDataValueCastToDouble(m_nNoDataUInt64);
380 : }
381 :
382 230 : if (pbSuccess)
383 230 : *pbSuccess = m_bHasNoData;
384 230 : return m_dfNoData;
385 : }
386 :
387 : /************************************************************************/
388 : /* SetNoDataValue() */
389 : /************************************************************************/
390 :
391 92 : CPLErr PDS4RawRasterBand::SetNoDataValue(double dfNewNoData)
392 : {
393 92 : m_dfNoData = dfNewNoData;
394 92 : m_bHasNoData = true;
395 92 : return CE_None;
396 : }
397 :
398 : /************************************************************************/
399 : /* GetNoDataValueAsInt64() */
400 : /************************************************************************/
401 :
402 9 : int64_t PDS4RawRasterBand::GetNoDataValueAsInt64(int *pbSuccess)
403 : {
404 9 : if (pbSuccess)
405 9 : *pbSuccess = m_bHasNoDataInt64;
406 9 : return m_nNoDataInt64;
407 : }
408 :
409 : /************************************************************************/
410 : /* SetNoDataValueAsInt64() */
411 : /************************************************************************/
412 :
413 3 : CPLErr PDS4RawRasterBand::SetNoDataValueAsInt64(int64_t nNewNoData)
414 : {
415 3 : m_nNoDataInt64 = nNewNoData;
416 3 : m_bHasNoDataInt64 = true;
417 :
418 3 : PDS4Dataset *poGDS = cpl::down_cast<PDS4Dataset *>(poDS);
419 3 : if (poGDS->m_poExternalDS && eAccess == GA_Update)
420 0 : poGDS->m_poExternalDS->GetRasterBand(nBand)->SetNoDataValueAsInt64(
421 0 : nNewNoData);
422 :
423 3 : return CE_None;
424 : }
425 :
426 : /************************************************************************/
427 : /* GetNoDataValueAsUInt64() */
428 : /************************************************************************/
429 :
430 9 : uint64_t PDS4RawRasterBand::GetNoDataValueAsUInt64(int *pbSuccess)
431 : {
432 9 : if (pbSuccess)
433 9 : *pbSuccess = m_bHasNoDataUInt64;
434 9 : return m_nNoDataUInt64;
435 : }
436 :
437 : /************************************************************************/
438 : /* SetNoDataValueAsUInt64() */
439 : /************************************************************************/
440 :
441 3 : CPLErr PDS4RawRasterBand::SetNoDataValueAsUInt64(uint64_t nNewNoData)
442 : {
443 3 : m_nNoDataUInt64 = nNewNoData;
444 3 : m_bHasNoDataUInt64 = true;
445 :
446 3 : PDS4Dataset *poGDS = cpl::down_cast<PDS4Dataset *>(poDS);
447 3 : if (poGDS->m_poExternalDS && eAccess == GA_Update)
448 0 : poGDS->m_poExternalDS->GetRasterBand(nBand)->SetNoDataValueAsUInt64(
449 0 : nNewNoData);
450 :
451 3 : return CE_None;
452 : }
453 :
454 : /************************************************************************/
455 : /* IReadBlock() */
456 : /************************************************************************/
457 :
458 132 : CPLErr PDS4RawRasterBand::IWriteBlock(int nXBlock, int nYBlock, void *pImage)
459 :
460 : {
461 132 : PDS4Dataset *poGDS = cpl::down_cast<PDS4Dataset *>(poDS);
462 132 : if (poGDS->m_bMustInitImageFile)
463 : {
464 0 : if (!poGDS->InitImageFile())
465 0 : return CE_Failure;
466 : }
467 :
468 132 : return RawRasterBand::IWriteBlock(nXBlock, nYBlock, pImage);
469 : }
470 :
471 : /************************************************************************/
472 : /* IRasterIO() */
473 : /************************************************************************/
474 :
475 1106 : CPLErr PDS4RawRasterBand::IRasterIO(GDALRWFlag eRWFlag, int nXOff, int nYOff,
476 : int nXSize, int nYSize, void *pData,
477 : int nBufXSize, int nBufYSize,
478 : GDALDataType eBufType, GSpacing nPixelSpace,
479 : GSpacing nLineSpace,
480 : GDALRasterIOExtraArg *psExtraArg)
481 :
482 : {
483 1106 : PDS4Dataset *poGDS = cpl::down_cast<PDS4Dataset *>(poDS);
484 1106 : if (eRWFlag == GF_Write && poGDS->m_bMustInitImageFile)
485 : {
486 10 : if (!poGDS->InitImageFile())
487 0 : return CE_Failure;
488 : }
489 :
490 1106 : return RawRasterBand::IRasterIO(eRWFlag, nXOff, nYOff, nXSize, nYSize,
491 : pData, nBufXSize, nBufYSize, eBufType,
492 1106 : nPixelSpace, nLineSpace, psExtraArg);
493 : }
494 :
495 : /************************************************************************/
496 : /* PDS4MaskBand() */
497 : /************************************************************************/
498 :
499 43 : PDS4MaskBand::PDS4MaskBand(GDALRasterBand *poBaseBand,
500 43 : const std::vector<double> &adfConstants)
501 43 : : m_poBaseBand(poBaseBand), m_pBuffer(nullptr), m_adfConstants(adfConstants)
502 : {
503 43 : eDataType = GDT_UInt8;
504 43 : poBaseBand->GetBlockSize(&nBlockXSize, &nBlockYSize);
505 43 : nRasterXSize = poBaseBand->GetXSize();
506 43 : nRasterYSize = poBaseBand->GetYSize();
507 43 : }
508 :
509 : /************************************************************************/
510 : /* ~PDS4MaskBand() */
511 : /************************************************************************/
512 :
513 86 : PDS4MaskBand::~PDS4MaskBand()
514 : {
515 43 : VSIFree(m_pBuffer);
516 86 : }
517 :
518 : /************************************************************************/
519 : /* FillMask() */
520 : /************************************************************************/
521 :
522 : template <class T>
523 88 : static void FillMask(void *pvBuffer, GByte *pabyDst, int nReqXSize,
524 : int nReqYSize, int nBlockXSize,
525 : const std::vector<double> &adfConstants)
526 : {
527 88 : const T *pSrc = static_cast<T *>(pvBuffer);
528 176 : std::vector<T> aConstants;
529 248 : for (size_t i = 0; i < adfConstants.size(); i++)
530 : {
531 : T cst;
532 160 : GDALCopyWord(adfConstants[i], cst);
533 160 : aConstants.push_back(cst);
534 : }
535 :
536 176 : for (int y = 0; y < nReqYSize; y++)
537 : {
538 1696 : for (int x = 0; x < nReqXSize; x++)
539 : {
540 1608 : const T nSrc = pSrc[y * nBlockXSize + x];
541 1608 : if (std::find(aConstants.begin(), aConstants.end(), nSrc) !=
542 : aConstants.end())
543 : {
544 6 : pabyDst[y * nBlockXSize + x] = 0;
545 : }
546 : else
547 : {
548 1602 : pabyDst[y * nBlockXSize + x] = 255;
549 : }
550 : }
551 : }
552 88 : }
553 :
554 : /************************************************************************/
555 : /* IReadBlock() */
556 : /************************************************************************/
557 :
558 88 : CPLErr PDS4MaskBand::IReadBlock(int nBlockXOff, int nBlockYOff, void *pImage)
559 :
560 : {
561 88 : const GDALDataType eSrcDT = m_poBaseBand->GetRasterDataType();
562 88 : const int nSrcDTSize = GDALGetDataTypeSizeBytes(eSrcDT);
563 88 : if (m_pBuffer == nullptr)
564 : {
565 12 : m_pBuffer = VSI_MALLOC3_VERBOSE(nBlockXSize, nBlockYSize, nSrcDTSize);
566 12 : if (m_pBuffer == nullptr)
567 0 : return CE_Failure;
568 : }
569 :
570 88 : const int nXOff = nBlockXOff * nBlockXSize;
571 88 : const int nReqXSize = std::min(nBlockXSize, nRasterXSize - nXOff);
572 88 : const int nYOff = nBlockYOff * nBlockYSize;
573 88 : const int nReqYSize = std::min(nBlockYSize, nRasterYSize - nYOff);
574 :
575 176 : if (m_poBaseBand->RasterIO(GF_Read, nXOff, nYOff, nReqXSize, nReqYSize,
576 : m_pBuffer, nReqXSize, nReqYSize, eSrcDT,
577 : nSrcDTSize,
578 88 : static_cast<GSpacing>(nSrcDTSize) * nBlockXSize,
579 88 : nullptr) != CE_None)
580 : {
581 0 : return CE_Failure;
582 : }
583 :
584 88 : GByte *pabyDst = static_cast<GByte *>(pImage);
585 88 : if (eSrcDT == GDT_UInt8)
586 : {
587 81 : FillMask<GByte>(m_pBuffer, pabyDst, nReqXSize, nReqYSize, nBlockXSize,
588 81 : m_adfConstants);
589 : }
590 7 : else if (eSrcDT == GDT_Int8)
591 : {
592 1 : FillMask<GInt8>(m_pBuffer, pabyDst, nReqXSize, nReqYSize, nBlockXSize,
593 1 : m_adfConstants);
594 : }
595 6 : else if (eSrcDT == GDT_UInt16)
596 : {
597 1 : FillMask<GUInt16>(m_pBuffer, pabyDst, nReqXSize, nReqYSize, nBlockXSize,
598 1 : m_adfConstants);
599 : }
600 5 : else if (eSrcDT == GDT_Int16)
601 : {
602 1 : FillMask<GInt16>(m_pBuffer, pabyDst, nReqXSize, nReqYSize, nBlockXSize,
603 1 : m_adfConstants);
604 : }
605 4 : else if (eSrcDT == GDT_UInt32)
606 : {
607 1 : FillMask<GUInt32>(m_pBuffer, pabyDst, nReqXSize, nReqYSize, nBlockXSize,
608 1 : m_adfConstants);
609 : }
610 3 : else if (eSrcDT == GDT_Int32)
611 : {
612 1 : FillMask<GInt32>(m_pBuffer, pabyDst, nReqXSize, nReqYSize, nBlockXSize,
613 1 : m_adfConstants);
614 : }
615 2 : else if (eSrcDT == GDT_UInt64)
616 : {
617 0 : FillMask<uint64_t>(m_pBuffer, pabyDst, nReqXSize, nReqYSize,
618 0 : nBlockXSize, m_adfConstants);
619 : }
620 2 : else if (eSrcDT == GDT_Int64)
621 : {
622 0 : FillMask<int64_t>(m_pBuffer, pabyDst, nReqXSize, nReqYSize, nBlockXSize,
623 0 : m_adfConstants);
624 : }
625 2 : else if (eSrcDT == GDT_Float32)
626 : {
627 1 : FillMask<float>(m_pBuffer, pabyDst, nReqXSize, nReqYSize, nBlockXSize,
628 1 : m_adfConstants);
629 : }
630 1 : else if (eSrcDT == GDT_Float64)
631 : {
632 1 : FillMask<double>(m_pBuffer, pabyDst, nReqXSize, nReqYSize, nBlockXSize,
633 1 : m_adfConstants);
634 : }
635 :
636 88 : return CE_None;
637 : }
638 :
639 : /************************************************************************/
640 : /* PDS4Dataset() */
641 : /************************************************************************/
642 :
643 490 : PDS4Dataset::PDS4Dataset()
644 : {
645 490 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
646 490 : }
647 :
648 : /************************************************************************/
649 : /* ~PDS4Dataset() */
650 : /************************************************************************/
651 :
652 980 : PDS4Dataset::~PDS4Dataset()
653 : {
654 490 : PDS4Dataset::Close();
655 980 : }
656 :
657 : /************************************************************************/
658 : /* Close() */
659 : /************************************************************************/
660 :
661 844 : CPLErr PDS4Dataset::Close(GDALProgressFunc, void *)
662 : {
663 844 : CPLErr eErr = CE_None;
664 844 : if (nOpenFlags != OPEN_FLAGS_CLOSED)
665 : {
666 490 : if (m_bMustInitImageFile)
667 : {
668 72 : if (!InitImageFile())
669 0 : eErr = CE_Failure;
670 : }
671 :
672 490 : if (PDS4Dataset::FlushCache(true) != CE_None)
673 0 : eErr = CE_Failure;
674 :
675 490 : if (m_bCreateHeader || m_bDirtyHeader)
676 193 : WriteHeader();
677 490 : if (m_fpImage)
678 337 : VSIFCloseL(m_fpImage);
679 490 : CSLDestroy(m_papszCreationOptions);
680 490 : PDS4Dataset::CloseDependentDatasets();
681 :
682 490 : if (GDALPamDataset::Close() != CE_None)
683 0 : eErr = CE_Failure;
684 : }
685 844 : return eErr;
686 : }
687 :
688 : /************************************************************************/
689 : /* GetRawBinaryLayout() */
690 : /************************************************************************/
691 :
692 1 : bool PDS4Dataset::GetRawBinaryLayout(GDALDataset::RawBinaryLayout &sLayout)
693 : {
694 1 : if (!RawDataset::GetRawBinaryLayout(sLayout))
695 0 : return false;
696 1 : sLayout.osRawFilename = m_osImageFilename;
697 1 : return true;
698 : }
699 :
700 : /************************************************************************/
701 : /* CloseDependentDatasets() */
702 : /************************************************************************/
703 :
704 490 : int PDS4Dataset::CloseDependentDatasets()
705 : {
706 490 : int bHasDroppedRef = GDALPamDataset::CloseDependentDatasets();
707 :
708 490 : if (m_poExternalDS)
709 : {
710 20 : bHasDroppedRef = FALSE;
711 20 : delete m_poExternalDS;
712 20 : m_poExternalDS = nullptr;
713 :
714 46 : for (int iBand = 0; iBand < nBands; iBand++)
715 : {
716 26 : delete papoBands[iBand];
717 26 : papoBands[iBand] = nullptr;
718 : }
719 20 : nBands = 0;
720 : }
721 :
722 490 : return bHasDroppedRef;
723 : }
724 :
725 : /************************************************************************/
726 : /* GetSpatialRef() */
727 : /************************************************************************/
728 :
729 46 : const OGRSpatialReference *PDS4Dataset::GetSpatialRef() const
730 : {
731 46 : if (!m_oSRS.IsEmpty())
732 41 : return &m_oSRS;
733 5 : return GDALPamDataset::GetSpatialRef();
734 : }
735 :
736 : /************************************************************************/
737 : /* SetSpatialRef() */
738 : /************************************************************************/
739 :
740 101 : CPLErr PDS4Dataset::SetSpatialRef(const OGRSpatialReference *poSRS)
741 :
742 : {
743 101 : if (eAccess == GA_ReadOnly)
744 0 : return CE_Failure;
745 101 : m_oSRS.Clear();
746 101 : if (poSRS)
747 101 : m_oSRS = *poSRS;
748 101 : if (m_poExternalDS)
749 10 : m_poExternalDS->SetSpatialRef(poSRS);
750 101 : return CE_None;
751 : }
752 :
753 : /************************************************************************/
754 : /* GetGeoTransform() */
755 : /************************************************************************/
756 :
757 50 : CPLErr PDS4Dataset::GetGeoTransform(GDALGeoTransform >) const
758 :
759 : {
760 50 : if (m_bGotTransform)
761 : {
762 47 : gt = m_gt;
763 47 : return CE_None;
764 : }
765 :
766 3 : return GDALPamDataset::GetGeoTransform(gt);
767 : }
768 :
769 : /************************************************************************/
770 : /* SetGeoTransform() */
771 : /************************************************************************/
772 :
773 101 : CPLErr PDS4Dataset::SetGeoTransform(const GDALGeoTransform >)
774 :
775 : {
776 101 : if (!((gt.xscale > 0.0 && gt.xrot == 0.0 && gt.yrot == 0.0 &&
777 100 : gt.yscale < 0.0) ||
778 1 : (gt.xscale == 0.0 && gt.xrot > 0.0 && gt.yrot > 0.0 &&
779 1 : gt.yscale == 0.0)))
780 : {
781 0 : CPLError(CE_Failure, CPLE_NotSupported,
782 : "Only north-up geotransform or map_projection_rotation=90 "
783 : "supported");
784 0 : return CE_Failure;
785 : }
786 101 : m_gt = gt;
787 101 : m_bGotTransform = true;
788 101 : if (m_poExternalDS)
789 10 : m_poExternalDS->SetGeoTransform(m_gt);
790 101 : return CE_None;
791 : }
792 :
793 : /************************************************************************/
794 : /* SetMetadata() */
795 : /************************************************************************/
796 :
797 8 : CPLErr PDS4Dataset::SetMetadata(CSLConstList papszMD, const char *pszDomain)
798 : {
799 8 : if (m_bUseSrcLabel && eAccess == GA_Update && pszDomain != nullptr &&
800 8 : EQUAL(pszDomain, "xml:PDS4"))
801 : {
802 6 : if (papszMD != nullptr && papszMD[0] != nullptr)
803 : {
804 6 : m_osXMLPDS4 = papszMD[0];
805 : }
806 6 : return CE_None;
807 : }
808 2 : return GDALPamDataset::SetMetadata(papszMD, pszDomain);
809 : }
810 :
811 : /************************************************************************/
812 : /* GetFileList() */
813 : /************************************************************************/
814 :
815 139 : char **PDS4Dataset::GetFileList()
816 : {
817 139 : char **papszFileList = GDALPamDataset::GetFileList();
818 276 : if (!m_osXMLFilename.empty() &&
819 137 : CSLFindString(papszFileList, m_osXMLFilename) < 0)
820 : {
821 0 : papszFileList = CSLAddString(papszFileList, m_osXMLFilename);
822 : }
823 139 : if (!m_osImageFilename.empty())
824 : {
825 90 : papszFileList = CSLAddString(papszFileList, m_osImageFilename);
826 : }
827 203 : for (const auto &poLayer : m_apoLayers)
828 : {
829 64 : auto papszTemp = poLayer->GetFileList();
830 64 : papszFileList = CSLInsertStrings(papszFileList, -1, papszTemp);
831 64 : CSLDestroy(papszTemp);
832 : }
833 139 : return papszFileList;
834 : }
835 :
836 : /************************************************************************/
837 : /* GetLinearValue() */
838 : /************************************************************************/
839 :
840 : static const struct
841 : {
842 : const char *pszUnit;
843 : double dfToMeter;
844 : } apsLinearUnits[] = {
845 : {"AU", 149597870700.0}, {"Angstrom", 1e-10}, {"cm", 1e-2}, {"km", 1e3},
846 : {"micrometer", 1e-6}, {"mm", 1e-3}, {"nm", 1e-9}};
847 :
848 802 : static double GetLinearValue(const CPLXMLNode *psParent,
849 : const char *pszElementName)
850 : {
851 802 : const CPLXMLNode *psNode = CPLGetXMLNode(psParent, pszElementName);
852 802 : if (psNode == nullptr)
853 0 : return 0.0;
854 802 : double dfVal = CPLAtof(CPLGetXMLValue(psNode, nullptr, ""));
855 802 : const char *pszUnit = CPLGetXMLValue(psNode, "unit", nullptr);
856 802 : if (pszUnit && !EQUAL(pszUnit, "m"))
857 : {
858 21 : bool bFound = false;
859 84 : for (size_t i = 0; i < CPL_ARRAYSIZE(apsLinearUnits); i++)
860 : {
861 84 : if (EQUAL(pszUnit, apsLinearUnits[i].pszUnit))
862 : {
863 21 : dfVal *= apsLinearUnits[i].dfToMeter;
864 21 : bFound = true;
865 21 : break;
866 : }
867 : }
868 21 : if (!bFound)
869 : {
870 0 : CPLError(CE_Warning, CPLE_AppDefined, "Unknown unit '%s' for '%s'",
871 : pszUnit, pszElementName);
872 : }
873 : }
874 802 : return dfVal;
875 : }
876 :
877 : /************************************************************************/
878 : /* GetResolutionValue() */
879 : /************************************************************************/
880 :
881 : static const struct
882 : {
883 : const char *pszUnit;
884 : double dfToMeter;
885 : } apsResolutionUnits[] = {
886 : {"km/pixel", 1e3},
887 : {"mm/pixel", 1e-3},
888 : };
889 :
890 316 : static double GetResolutionValue(CPLXMLNode *psParent,
891 : const char *pszElementName)
892 : {
893 316 : CPLXMLNode *psNode = CPLGetXMLNode(psParent, pszElementName);
894 316 : if (psNode == nullptr)
895 0 : return 0.0;
896 316 : double dfVal = CPLAtof(CPLGetXMLValue(psNode, nullptr, ""));
897 316 : const char *pszUnit = CPLGetXMLValue(psNode, "unit", nullptr);
898 316 : if (pszUnit && !EQUAL(pszUnit, "m/pixel"))
899 : {
900 21 : bool bFound = false;
901 21 : for (size_t i = 0; i < CPL_ARRAYSIZE(apsResolutionUnits); i++)
902 : {
903 21 : if (EQUAL(pszUnit, apsResolutionUnits[i].pszUnit))
904 : {
905 21 : dfVal *= apsResolutionUnits[i].dfToMeter;
906 21 : bFound = true;
907 21 : break;
908 : }
909 : }
910 21 : if (!bFound)
911 : {
912 0 : CPLError(CE_Warning, CPLE_AppDefined, "Unknown unit '%s' for '%s'",
913 : pszUnit, pszElementName);
914 : }
915 : }
916 316 : return dfVal;
917 : }
918 :
919 : /************************************************************************/
920 : /* GetAngularValue() */
921 : /************************************************************************/
922 :
923 : static const struct
924 : {
925 : const char *pszUnit;
926 : double dfToDeg;
927 : } apsAngularUnits[] = {{"arcmin", 1. / 60.},
928 : {"arcsec", 1. / 3600},
929 : {"hr", 15.0},
930 : {"mrad", 180.0 / M_PI / 1000.},
931 : {"rad", 180.0 / M_PI}};
932 :
933 819 : static double GetAngularValue(CPLXMLNode *psParent, const char *pszElementName,
934 : bool *pbGotVal = nullptr)
935 : {
936 819 : CPLXMLNode *psNode = CPLGetXMLNode(psParent, pszElementName);
937 819 : if (psNode == nullptr)
938 : {
939 424 : if (pbGotVal)
940 265 : *pbGotVal = false;
941 424 : return 0.0;
942 : }
943 395 : double dfVal = CPLAtof(CPLGetXMLValue(psNode, nullptr, ""));
944 395 : const char *pszUnit = CPLGetXMLValue(psNode, "unit", nullptr);
945 395 : if (pszUnit && !EQUAL(pszUnit, "deg"))
946 : {
947 0 : bool bFound = false;
948 0 : for (size_t i = 0; i < CPL_ARRAYSIZE(apsAngularUnits); i++)
949 : {
950 0 : if (EQUAL(pszUnit, apsAngularUnits[i].pszUnit))
951 : {
952 0 : dfVal *= apsAngularUnits[i].dfToDeg;
953 0 : bFound = true;
954 0 : break;
955 : }
956 : }
957 0 : if (!bFound)
958 : {
959 0 : CPLError(CE_Warning, CPLE_AppDefined, "Unknown unit '%s' for '%s'",
960 : pszUnit, pszElementName);
961 : }
962 : }
963 395 : if (pbGotVal)
964 220 : *pbGotVal = true;
965 395 : return dfVal;
966 : }
967 :
968 : /************************************************************************/
969 : /* ReadGeoreferencing() */
970 : /************************************************************************/
971 :
972 : // See https://pds.nasa.gov/pds4/cart/v1/PDS4_CART_1G00_1950.xsd, (GDAL 3.4)
973 : // https://pds.nasa.gov/pds4/cart/v1/PDS4_CART_1D00_1933.xsd,
974 : // https://raw.githubusercontent.com/nasa-pds-data-dictionaries/ldd-cart/master/build/1.B.0.0/PDS4_CART_1B00.xsd,
975 : // https://pds.nasa.gov/pds4/cart/v1/PDS4_CART_1700.xsd
976 : // and the corresponding .sch files
977 295 : void PDS4Dataset::ReadGeoreferencing(CPLXMLNode *psProduct)
978 : {
979 295 : CPLXMLNode *psCart = CPLGetXMLNode(
980 : psProduct, "Observation_Area.Discipline_Area.Cartography");
981 295 : if (psCart == nullptr)
982 : {
983 133 : CPLDebug("PDS4",
984 : "Did not find Observation_Area.Discipline_Area.Cartography");
985 133 : return;
986 : }
987 :
988 : // Bounding box: informative only
989 : CPLXMLNode *psBounding =
990 162 : CPLGetXMLNode(psCart, "Spatial_Domain.Bounding_Coordinates");
991 162 : if (psBounding)
992 : {
993 : const char *pszWest =
994 162 : CPLGetXMLValue(psBounding, "west_bounding_coordinate", nullptr);
995 : const char *pszEast =
996 162 : CPLGetXMLValue(psBounding, "east_bounding_coordinate", nullptr);
997 : const char *pszNorth =
998 162 : CPLGetXMLValue(psBounding, "north_bounding_coordinate", nullptr);
999 : const char *pszSouth =
1000 162 : CPLGetXMLValue(psBounding, "south_bounding_coordinate", nullptr);
1001 162 : if (pszWest)
1002 162 : CPLDebug("PDS4", "West: %s", pszWest);
1003 162 : if (pszEast)
1004 162 : CPLDebug("PDS4", "East: %s", pszEast);
1005 162 : if (pszNorth)
1006 162 : CPLDebug("PDS4", "North: %s", pszNorth);
1007 162 : if (pszSouth)
1008 162 : CPLDebug("PDS4", "South: %s", pszSouth);
1009 : }
1010 :
1011 : CPLXMLNode *psSR =
1012 162 : CPLGetXMLNode(psCart, "Spatial_Reference_Information.Horizontal_"
1013 : "Coordinate_System_Definition");
1014 162 : if (psSR == nullptr)
1015 : {
1016 0 : CPLDebug("PDS4", "Did not find Spatial_Reference_Information."
1017 : "Horizontal_Coordinate_System_Definition");
1018 0 : return;
1019 : }
1020 :
1021 162 : double dfLongitudeMultiplier = 1;
1022 162 : const CPLXMLNode *psGeodeticModel = CPLGetXMLNode(psSR, "Geodetic_Model");
1023 162 : if (psGeodeticModel != nullptr)
1024 : {
1025 162 : if (EQUAL(CPLGetXMLValue(psGeodeticModel, "longitude_direction", ""),
1026 : "Positive West"))
1027 : {
1028 3 : dfLongitudeMultiplier = -1;
1029 : }
1030 : }
1031 :
1032 324 : OGRSpatialReference oSRS;
1033 : CPLXMLNode *psGridCoordinateSystem =
1034 162 : CPLGetXMLNode(psSR, "Planar.Grid_Coordinate_System");
1035 162 : CPLXMLNode *psMapProjection = CPLGetXMLNode(psSR, "Planar.Map_Projection");
1036 324 : CPLString osProjName;
1037 162 : double dfCenterLon = 0.0;
1038 162 : double dfCenterLat = 0.0;
1039 162 : double dfStdParallel1 = 0.0;
1040 162 : double dfStdParallel2 = 0.0;
1041 162 : double dfScale = 1.0;
1042 162 : double dfMapProjectionRotation = 0.0;
1043 162 : if (psGridCoordinateSystem != nullptr)
1044 : {
1045 : osProjName = CPLGetXMLValue(psGridCoordinateSystem,
1046 0 : "grid_coordinate_system_name", "");
1047 0 : if (!osProjName.empty())
1048 : {
1049 0 : if (osProjName == "Universal Transverse Mercator")
1050 : {
1051 0 : CPLXMLNode *psUTMZoneNumber = CPLGetXMLNode(
1052 : psGridCoordinateSystem,
1053 : "Universal_Transverse_Mercator.utm_zone_number");
1054 0 : if (psUTMZoneNumber)
1055 : {
1056 : int nZone =
1057 0 : atoi(CPLGetXMLValue(psUTMZoneNumber, nullptr, ""));
1058 0 : oSRS.SetUTM(std::abs(nZone), nZone >= 0);
1059 : }
1060 : }
1061 0 : else if (osProjName == "Universal Polar Stereographic")
1062 : {
1063 0 : CPLXMLNode *psProjParamNode = CPLGetXMLNode(
1064 : psGridCoordinateSystem,
1065 : "Universal_Polar_Stereographic.Polar_Stereographic");
1066 0 : if (psProjParamNode)
1067 : {
1068 0 : dfCenterLon =
1069 0 : GetAngularValue(psProjParamNode,
1070 : "longitude_of_central_meridian") *
1071 : dfLongitudeMultiplier;
1072 0 : dfCenterLat = GetAngularValue(
1073 : psProjParamNode, "latitude_of_projection_origin");
1074 0 : dfScale = CPLAtof(CPLGetXMLValue(
1075 : psProjParamNode, "scale_factor_at_projection_origin",
1076 : "1"));
1077 0 : oSRS.SetPS(dfCenterLat, dfCenterLon, dfScale, 0, 0);
1078 : }
1079 : }
1080 : else
1081 : {
1082 0 : CPLError(CE_Warning, CPLE_NotSupported,
1083 : "grid_coordinate_system_name = %s not supported",
1084 : osProjName.c_str());
1085 : }
1086 : }
1087 : }
1088 162 : else if (psMapProjection != nullptr)
1089 : {
1090 161 : osProjName = CPLGetXMLValue(psMapProjection, "map_projection_name", "");
1091 161 : if (!osProjName.empty())
1092 : {
1093 161 : CPLXMLNode *psProjParamNode = CPLGetXMLNode(
1094 : psMapProjection,
1095 322 : CPLString(osProjName).replaceAll(' ', '_').c_str());
1096 161 : if (psProjParamNode == nullptr &&
1097 : // typo in https://pds.nasa.gov/pds4/cart/v1/PDS4_CART_1700.sch
1098 0 : EQUAL(osProjName, "Orothographic"))
1099 : {
1100 : psProjParamNode =
1101 0 : CPLGetXMLNode(psMapProjection, "Orthographic");
1102 : }
1103 161 : bool bGotStdParallel1 = false;
1104 161 : bool bGotStdParallel2 = false;
1105 161 : bool bGotScale = false;
1106 161 : if (psProjParamNode)
1107 : {
1108 161 : bool bGotCenterLon = false;
1109 161 : dfCenterLon = GetAngularValue(psProjParamNode,
1110 : "longitude_of_central_meridian",
1111 : &bGotCenterLon) *
1112 : dfLongitudeMultiplier;
1113 161 : if (!bGotCenterLon)
1114 : {
1115 2 : dfCenterLon =
1116 2 : GetAngularValue(psProjParamNode,
1117 : "straight_vertical_longitude_from_pole",
1118 : &bGotCenterLon) *
1119 : dfLongitudeMultiplier;
1120 : }
1121 161 : dfCenterLat = GetAngularValue(psProjParamNode,
1122 : "latitude_of_projection_origin");
1123 161 : dfStdParallel1 = GetAngularValue(
1124 : psProjParamNode, "standard_parallel_1", &bGotStdParallel1);
1125 161 : dfStdParallel2 = GetAngularValue(
1126 : psProjParamNode, "standard_parallel_2", &bGotStdParallel2);
1127 : const char *pszScaleParam =
1128 161 : (osProjName == "Transverse Mercator")
1129 161 : ? "scale_factor_at_central_meridian"
1130 161 : : "scale_factor_at_projection_origin";
1131 : const char *pszScaleVal =
1132 161 : CPLGetXMLValue(psProjParamNode, pszScaleParam, nullptr);
1133 161 : bGotScale = pszScaleVal != nullptr;
1134 161 : dfScale = (pszScaleVal) ? CPLAtof(pszScaleVal) : 1.0;
1135 :
1136 : dfMapProjectionRotation =
1137 161 : GetAngularValue(psProjParamNode, "map_projection_rotation");
1138 : }
1139 :
1140 : CPLXMLNode *psObliqueAzimuth =
1141 161 : CPLGetXMLNode(psProjParamNode, "Oblique_Line_Azimuth");
1142 : CPLXMLNode *psObliquePoint =
1143 161 : CPLGetXMLNode(psProjParamNode, "Oblique_Line_Point");
1144 :
1145 161 : if (EQUAL(osProjName, "Equirectangular"))
1146 : {
1147 53 : oSRS.SetEquirectangular2(dfCenterLat, dfCenterLon,
1148 : dfStdParallel1, 0, 0);
1149 : }
1150 108 : else if (EQUAL(osProjName, "Lambert Conformal Conic"))
1151 : {
1152 4 : if (bGotScale)
1153 : {
1154 2 : if ((bGotStdParallel1 && dfStdParallel1 != dfCenterLat) ||
1155 0 : (bGotStdParallel2 && dfStdParallel2 != dfCenterLat))
1156 : {
1157 0 : CPLError(
1158 : CE_Warning, CPLE_AppDefined,
1159 : "Ignoring standard_parallel_1 and/or "
1160 : "standard_parallel_2 with LCC_1SP formulation");
1161 : }
1162 2 : oSRS.SetLCC1SP(dfCenterLat, dfCenterLon, dfScale, 0, 0);
1163 : }
1164 : else
1165 : {
1166 2 : oSRS.SetLCC(dfStdParallel1, dfStdParallel2, dfCenterLat,
1167 : dfCenterLon, 0, 0);
1168 : }
1169 : }
1170 104 : else if (EQUAL(osProjName, "Mercator"))
1171 : {
1172 6 : if (bGotScale)
1173 : {
1174 4 : oSRS.SetMercator(dfCenterLat, // should be 0 normally
1175 : dfCenterLon, dfScale, 0.0, 0.0);
1176 : }
1177 : else
1178 : {
1179 2 : oSRS.SetMercator2SP(dfStdParallel1,
1180 : dfCenterLat, // should be 0 normally
1181 : dfCenterLon, 0.0, 0.0);
1182 : }
1183 : }
1184 98 : else if (EQUAL(osProjName, "Orthographic"))
1185 : {
1186 2 : oSRS.SetOrthographic(dfCenterLat, dfCenterLon, 0.0, 0.0);
1187 : }
1188 98 : else if (EQUAL(osProjName, "Oblique Mercator") &&
1189 2 : (psObliqueAzimuth != nullptr || psObliquePoint != nullptr))
1190 : {
1191 4 : if (psObliqueAzimuth)
1192 : {
1193 : // Not sure of this
1194 2 : dfCenterLon = CPLAtof(
1195 : CPLGetXMLValue(psObliqueAzimuth,
1196 : "azimuth_measure_point_longitude", "0"));
1197 :
1198 2 : double dfAzimuth = CPLAtof(CPLGetXMLValue(
1199 : psObliqueAzimuth, "azimuthal_angle", "0"));
1200 2 : oSRS.SetProjection(
1201 : SRS_PT_HOTINE_OBLIQUE_MERCATOR_AZIMUTH_CENTER);
1202 2 : oSRS.SetNormProjParm(SRS_PP_LATITUDE_OF_CENTER,
1203 : dfCenterLat);
1204 2 : oSRS.SetNormProjParm(SRS_PP_LONGITUDE_OF_CENTER,
1205 : dfCenterLon);
1206 2 : oSRS.SetNormProjParm(SRS_PP_AZIMUTH, dfAzimuth);
1207 : // SetNormProjParm( SRS_PP_RECTIFIED_GRID_ANGLE,
1208 : // dfRectToSkew );
1209 2 : oSRS.SetNormProjParm(SRS_PP_SCALE_FACTOR, dfScale);
1210 2 : oSRS.SetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
1211 2 : oSRS.SetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
1212 : }
1213 : else
1214 : {
1215 2 : double dfLat1 = 0.0;
1216 2 : double dfLong1 = 0.0;
1217 2 : double dfLat2 = 0.0;
1218 2 : double dfLong2 = 0.0;
1219 2 : CPLXMLNode *psPoint = CPLGetXMLNode(
1220 : psObliquePoint, "Oblique_Line_Point_Group");
1221 2 : if (psPoint)
1222 : {
1223 2 : dfLat1 = CPLAtof(CPLGetXMLValue(
1224 : psPoint, "oblique_line_latitude", "0.0"));
1225 2 : dfLong1 = CPLAtof(CPLGetXMLValue(
1226 : psPoint, "oblique_line_longitude", "0.0"));
1227 2 : psPoint = psPoint->psNext;
1228 2 : if (psPoint && psPoint->eType == CXT_Element &&
1229 2 : EQUAL(psPoint->pszValue,
1230 : "Oblique_Line_Point_Group"))
1231 : {
1232 2 : dfLat2 = CPLAtof(CPLGetXMLValue(
1233 : psPoint, "oblique_line_latitude", "0.0"));
1234 2 : dfLong2 = CPLAtof(CPLGetXMLValue(
1235 : psPoint, "oblique_line_longitude", "0.0"));
1236 : }
1237 : }
1238 2 : oSRS.SetHOM2PNO(dfCenterLat, dfLat1, dfLong1, dfLat2,
1239 : dfLong2, dfScale, 0.0, 0.0);
1240 : }
1241 : }
1242 92 : else if (EQUAL(osProjName, "Polar Stereographic"))
1243 : {
1244 2 : oSRS.SetPS(dfCenterLat, dfCenterLon, dfScale, 0, 0);
1245 : }
1246 90 : else if (EQUAL(osProjName, "Polyconic"))
1247 : {
1248 2 : oSRS.SetPolyconic(dfCenterLat, dfCenterLon, 0, 0);
1249 : }
1250 88 : else if (EQUAL(osProjName, "Sinusoidal"))
1251 : {
1252 4 : oSRS.SetSinusoidal(dfCenterLon, 0, 0);
1253 : }
1254 84 : else if (EQUAL(osProjName, "Transverse Mercator"))
1255 : {
1256 78 : oSRS.SetTM(dfCenterLat, dfCenterLon, dfScale, 0, 0);
1257 : }
1258 :
1259 : // Below values are valid map_projection_name according to
1260 : // the schematron but they don't have a dedicated element to
1261 : // hold the projection parameter. Assumed the schema is extended
1262 : // similarly to the existing for a few obvious ones
1263 6 : else if (EQUAL(osProjName, "Albers Conical Equal Area"))
1264 : {
1265 0 : oSRS.SetACEA(dfStdParallel1, dfStdParallel2, dfCenterLat,
1266 : dfCenterLon, 0.0, 0.0);
1267 : }
1268 6 : else if (EQUAL(osProjName, "Azimuthal Equidistant"))
1269 : {
1270 0 : oSRS.SetAE(dfCenterLat, dfCenterLon, 0, 0);
1271 : }
1272 6 : else if (EQUAL(osProjName, "Equidistant Conic"))
1273 : {
1274 0 : oSRS.SetEC(dfStdParallel1, dfStdParallel2, dfCenterLat,
1275 : dfCenterLon, 0.0, 0.0);
1276 : }
1277 : // Unhandled: General Vertical Near-sided Projection
1278 6 : else if (EQUAL(osProjName, "Gnomonic"))
1279 : {
1280 0 : oSRS.SetGnomonic(dfCenterLat, dfCenterLon, 0, 0);
1281 : }
1282 6 : else if (EQUAL(osProjName, "Lambert Azimuthal Equal Area"))
1283 : {
1284 2 : oSRS.SetLAEA(dfCenterLat, dfCenterLon, 0, 0);
1285 : }
1286 4 : else if (EQUAL(osProjName, "Miller Cylindrical"))
1287 : {
1288 0 : oSRS.SetMC(dfCenterLat, dfCenterLon, 0, 0);
1289 : }
1290 4 : else if (EQUAL(osProjName, "Orothographic") // typo
1291 4 : || EQUAL(osProjName, "Orthographic"))
1292 : {
1293 0 : osProjName = "Orthographic";
1294 0 : oSRS.SetOrthographic(dfCenterLat, dfCenterLon, 0, 0);
1295 : }
1296 4 : else if (EQUAL(osProjName, "Robinson"))
1297 : {
1298 0 : oSRS.SetRobinson(dfCenterLon, 0, 0);
1299 : }
1300 : // Unhandled: Space Oblique Mercator
1301 4 : else if (EQUAL(osProjName, "Stereographic"))
1302 : {
1303 0 : oSRS.SetStereographic(dfCenterLat, dfCenterLon, dfScale, 0, 0);
1304 : }
1305 4 : else if (EQUAL(osProjName, "van der Grinten"))
1306 : {
1307 0 : oSRS.SetVDG(dfCenterLon, 0, 0);
1308 : }
1309 4 : else if (EQUAL(osProjName, "Oblique Cylindrical"))
1310 : {
1311 4 : const double poleLatitude = GetAngularValue(
1312 : psProjParamNode, "oblique_proj_pole_latitude");
1313 : const double poleLongitude =
1314 4 : GetAngularValue(psProjParamNode,
1315 : "oblique_proj_pole_longitude") *
1316 4 : dfLongitudeMultiplier;
1317 4 : const double poleRotation = GetAngularValue(
1318 : psProjParamNode, "oblique_proj_pole_rotation");
1319 :
1320 8 : CPLString oProj4String;
1321 : // Cf isis3dataset.cpp comments for ObliqueCylindrical
1322 : oProj4String.Printf("+proj=ob_tran +o_proj=eqc +o_lon_p=%.17g "
1323 : "+o_lat_p=%.17g +lon_0=%.17g",
1324 : -poleRotation, 180 - poleLatitude,
1325 4 : poleLongitude);
1326 4 : oSRS.SetFromUserInput(oProj4String);
1327 : }
1328 : else
1329 : {
1330 0 : CPLError(CE_Warning, CPLE_NotSupported,
1331 : "map_projection_name = %s not supported",
1332 : osProjName.c_str());
1333 : }
1334 : }
1335 : }
1336 : else
1337 : {
1338 1 : CPLXMLNode *psGeographic = CPLGetXMLNode(psSR, "Geographic");
1339 1 : if (GetLayerCount() && psGeographic)
1340 : {
1341 : // do nothing
1342 : }
1343 : else
1344 : {
1345 0 : CPLError(CE_Warning, CPLE_AppDefined,
1346 : "Planar.Map_Projection not found");
1347 : }
1348 : }
1349 :
1350 162 : if (oSRS.IsProjected())
1351 : {
1352 161 : oSRS.SetLinearUnits("Metre", 1.0);
1353 : }
1354 :
1355 162 : if (psGeodeticModel != nullptr)
1356 : {
1357 : const char *pszLatitudeType =
1358 162 : CPLGetXMLValue(psGeodeticModel, "latitude_type", "");
1359 162 : bool bIsOgraphic = EQUAL(pszLatitudeType, "Planetographic");
1360 :
1361 : const bool bUseLDD1930RadiusNames =
1362 162 : CPLGetXMLNode(psGeodeticModel, "a_axis_radius") != nullptr;
1363 :
1364 : // Before PDS CART schema pre-1.B.10.0 (pre LDD version 1.9.3.0),
1365 : // the confusing semi_major_radius, semi_minor_radius and polar_radius
1366 : // were used but did not follow the recommended
1367 : // FGDC names. Using both "semi" and "radius" in the same keyword,
1368 : // which both mean half, does not make sense.
1369 162 : const char *pszAAxis =
1370 162 : bUseLDD1930RadiusNames ? "a_axis_radius" : "semi_major_radius";
1371 162 : const char *pszBAxis =
1372 162 : bUseLDD1930RadiusNames ? "b_axis_radius" : "semi_minor_radius";
1373 162 : const char *pszCAxis =
1374 162 : bUseLDD1930RadiusNames ? "c_axis_radius" : "polar_radius";
1375 :
1376 162 : const double dfSemiMajor = GetLinearValue(psGeodeticModel, pszAAxis);
1377 :
1378 : // a_axis_radius and b_axis_radius should be the same in most cases
1379 : // unless a triaxial body is being defined. This should be extremely
1380 : // rare (and not used) since the IAU generally defines a best-fit sphere
1381 : // for triaxial bodies: https://astrogeology.usgs.gov/groups/IAU-WGCCRE
1382 162 : const double dfBValue = GetLinearValue(psGeodeticModel, pszBAxis);
1383 162 : if (dfSemiMajor != dfBValue)
1384 : {
1385 0 : CPLError(CE_Warning, CPLE_AppDefined,
1386 : "%s = %f m, different from "
1387 : "%s = %f, will be ignored",
1388 : pszBAxis, dfBValue, pszAAxis, dfSemiMajor);
1389 : }
1390 :
1391 162 : const double dfPolarRadius = GetLinearValue(psGeodeticModel, pszCAxis);
1392 : // Use the polar_radius as the actual semi minor
1393 162 : const double dfSemiMinor = dfPolarRadius;
1394 :
1395 : // Compulsory
1396 162 : const char *pszTargetName = CPLGetXMLValue(
1397 : psProduct, "Observation_Area.Target_Identification.name",
1398 : "unknown");
1399 :
1400 162 : if (oSRS.IsProjected())
1401 : {
1402 483 : CPLString osProjTargetName = osProjName + " " + pszTargetName;
1403 161 : oSRS.SetProjCS(osProjTargetName);
1404 : }
1405 :
1406 486 : CPLString osGeogName = CPLString("GCS_") + pszTargetName;
1407 :
1408 : CPLString osSphereName =
1409 324 : CPLGetXMLValue(psGeodeticModel, "spheroid_name", pszTargetName);
1410 324 : CPLString osDatumName = "D_" + osSphereName;
1411 :
1412 : // calculate inverse flattening from major and minor axis: 1/f = a/(a-b)
1413 162 : double dfInvFlattening = 0;
1414 162 : if ((dfSemiMajor - dfSemiMinor) >= 0.00000001)
1415 : {
1416 82 : dfInvFlattening = dfSemiMajor / (dfSemiMajor - dfSemiMinor);
1417 : }
1418 :
1419 : //(if stereographic with center lat ==90) or (polar stereographic )
1420 324 : if ((EQUAL(osProjName, "STEREOGRAPHIC") && fabs(dfCenterLat) == 90) ||
1421 162 : (EQUAL(osProjName, "POLAR STEREOGRAPHIC")))
1422 : {
1423 2 : if (bIsOgraphic)
1424 : {
1425 0 : oSRS.SetGeogCS(osGeogName, osDatumName, osSphereName,
1426 : dfSemiMajor, dfInvFlattening,
1427 : "Reference_Meridian", 0.0);
1428 : }
1429 : else
1430 : {
1431 2 : osSphereName += "_polarRadius";
1432 2 : oSRS.SetGeogCS(osGeogName, osDatumName, osSphereName,
1433 : dfPolarRadius, 0.0, "Reference_Meridian", 0.0);
1434 : }
1435 : }
1436 160 : else if ((EQUAL(osProjName, "EQUIRECTANGULAR")) ||
1437 107 : (EQUAL(osProjName, "ORTHOGRAPHIC")) ||
1438 372 : (EQUAL(osProjName, "STEREOGRAPHIC")) ||
1439 105 : (EQUAL(osProjName, "SINUSOIDAL")))
1440 : {
1441 59 : oSRS.SetGeogCS(osGeogName, osDatumName, osSphereName, dfSemiMajor,
1442 : 0.0, "Reference_Meridian", 0.0);
1443 : }
1444 : else
1445 : {
1446 101 : if (bIsOgraphic)
1447 : {
1448 2 : oSRS.SetGeogCS(osGeogName, osDatumName, osSphereName,
1449 : dfSemiMajor, dfInvFlattening,
1450 : "Reference_Meridian", 0.0);
1451 : }
1452 : else
1453 : {
1454 99 : oSRS.SetGeogCS(osGeogName, osDatumName, osSphereName,
1455 : dfSemiMajor, 0.0, "Reference_Meridian", 0.0);
1456 : }
1457 : }
1458 : }
1459 :
1460 : CPLXMLNode *psPCI =
1461 162 : CPLGetXMLNode(psSR, "Planar.Planar_Coordinate_Information");
1462 162 : CPLXMLNode *psGT = CPLGetXMLNode(psSR, "Planar.Geo_Transformation");
1463 162 : if (psPCI && psGT)
1464 : {
1465 : const char *pszPCIEncoding =
1466 158 : CPLGetXMLValue(psPCI, "planar_coordinate_encoding_method", "");
1467 158 : CPLXMLNode *psCR = CPLGetXMLNode(psPCI, "Coordinate_Representation");
1468 158 : if (!EQUAL(pszPCIEncoding, "Coordinate Pair"))
1469 : {
1470 0 : CPLError(CE_Warning, CPLE_NotSupported,
1471 : "planar_coordinate_encoding_method = %s not supported",
1472 : pszPCIEncoding);
1473 : }
1474 158 : else if (psCR != nullptr)
1475 : {
1476 158 : double dfXRes = GetResolutionValue(psCR, "pixel_resolution_x");
1477 158 : double dfYRes = GetResolutionValue(psCR, "pixel_resolution_y");
1478 158 : double dfULX = GetLinearValue(psGT, "upperleft_corner_x");
1479 158 : double dfULY = GetLinearValue(psGT, "upperleft_corner_y");
1480 :
1481 : // The PDS4 specification is not really clear about the
1482 : // origin convention, but it appears from
1483 : // https://github.com/OSGeo/gdal/issues/735 that it matches GDAL
1484 : // top-left corner of top-left pixel
1485 158 : m_gt.xorig = dfULX;
1486 158 : m_gt.xscale = dfXRes;
1487 158 : m_gt.xrot = 0.0;
1488 158 : m_gt.yorig = dfULY;
1489 158 : m_gt.yrot = 0.0;
1490 158 : m_gt.yscale = -dfYRes;
1491 158 : m_bGotTransform = true;
1492 :
1493 158 : if (dfMapProjectionRotation != 0)
1494 : {
1495 4 : const double sin_rot =
1496 : dfMapProjectionRotation == 90
1497 4 : ? 1.0
1498 0 : : sin(dfMapProjectionRotation / 180 * M_PI);
1499 4 : const double cos_rot =
1500 : dfMapProjectionRotation == 90
1501 4 : ? 0.0
1502 0 : : cos(dfMapProjectionRotation / 180 * M_PI);
1503 4 : const double gt_1 = cos_rot * m_gt.xscale - sin_rot * m_gt.yrot;
1504 4 : const double gt_2 = cos_rot * m_gt.xrot - sin_rot * m_gt.yscale;
1505 4 : const double gt_0 = cos_rot * m_gt.xorig - sin_rot * m_gt.yorig;
1506 4 : const double gt_4 = sin_rot * m_gt.xscale + cos_rot * m_gt.yrot;
1507 4 : const double gt_5 = sin_rot * m_gt.xrot + cos_rot * m_gt.yscale;
1508 4 : const double gt_3 = sin_rot * m_gt.xorig + cos_rot * m_gt.yorig;
1509 4 : m_gt.xscale = gt_1;
1510 4 : m_gt.xrot = gt_2;
1511 4 : m_gt.xorig = gt_0;
1512 4 : m_gt.yrot = gt_4;
1513 4 : m_gt.yscale = gt_5;
1514 4 : m_gt.yorig = gt_3;
1515 : }
1516 : }
1517 : }
1518 :
1519 162 : if (!oSRS.IsEmpty())
1520 : {
1521 162 : if (GetRasterCount())
1522 : {
1523 158 : m_oSRS = std::move(oSRS);
1524 : }
1525 4 : else if (GetLayerCount())
1526 : {
1527 8 : for (auto &poLayer : m_apoLayers)
1528 : {
1529 4 : if (poLayer->GetGeomType() != wkbNone)
1530 : {
1531 4 : auto poSRSClone = oSRS.Clone();
1532 4 : poLayer->SetSpatialRef(poSRSClone);
1533 4 : poSRSClone->Release();
1534 : }
1535 : }
1536 : }
1537 : }
1538 : }
1539 :
1540 : /************************************************************************/
1541 : /* GetLayer() */
1542 : /************************************************************************/
1543 :
1544 319 : const OGRLayer *PDS4Dataset::GetLayer(int nIndex) const
1545 : {
1546 319 : if (nIndex < 0 || nIndex >= GetLayerCount())
1547 8 : return nullptr;
1548 311 : return dynamic_cast<const OGRLayer *>(m_apoLayers[nIndex].get());
1549 : }
1550 :
1551 : /************************************************************************/
1552 : /* FixupTableFilename() */
1553 : /************************************************************************/
1554 :
1555 97 : static std::string FixupTableFilename(const std::string &osFilename)
1556 : {
1557 : VSIStatBufL sStat;
1558 97 : if (VSIStatL(osFilename.c_str(), &sStat) == 0)
1559 : {
1560 97 : return osFilename;
1561 : }
1562 0 : const std::string osExt = CPLGetExtensionSafe(osFilename.c_str());
1563 0 : if (!osExt.empty())
1564 : {
1565 0 : std::string osTry(osFilename);
1566 0 : if (osExt[0] >= 'a' && osExt[0] <= 'z')
1567 : {
1568 0 : osTry = CPLResetExtensionSafe(osFilename.c_str(),
1569 0 : CPLString(osExt).toupper());
1570 : }
1571 : else
1572 : {
1573 0 : osTry = CPLResetExtensionSafe(osFilename.c_str(),
1574 0 : CPLString(osExt).tolower());
1575 : }
1576 0 : if (VSIStatL(osTry.c_str(), &sStat) == 0)
1577 : {
1578 0 : return osTry;
1579 : }
1580 : }
1581 0 : return osFilename;
1582 : }
1583 :
1584 : /************************************************************************/
1585 : /* OpenTableCharacter() */
1586 : /************************************************************************/
1587 :
1588 22 : bool PDS4Dataset::OpenTableCharacter(const char *pszFilename,
1589 : const CPLXMLNode *psTable)
1590 : {
1591 22 : if (CPLHasPathTraversal(pszFilename))
1592 : {
1593 0 : CPLError(CE_Failure, CPLE_NotSupported,
1594 : "OpenTableCharacter(): path traversal detected: %s",
1595 : pszFilename);
1596 0 : return false;
1597 : }
1598 44 : std::string osLayerName(CPLGetBasenameSafe(pszFilename));
1599 22 : if (cpl::starts_with(osLayerName,
1600 44 : CPLGetBasenameSafe(m_osXMLFilename.c_str()) + "_"))
1601 26 : osLayerName = osLayerName.substr(
1602 39 : CPLGetBasenameSafe(m_osXMLFilename.c_str()).size() + 1);
1603 :
1604 44 : const std::string osFullFilename = FixupTableFilename(CPLFormFilenameSafe(
1605 66 : CPLGetPathSafe(m_osXMLFilename.c_str()).c_str(), pszFilename, nullptr));
1606 : auto poLayer = std::make_unique<PDS4TableCharacter>(
1607 22 : this, osLayerName.c_str(), osFullFilename.c_str(),
1608 66 : eAccess == GA_Update);
1609 22 : if (!poLayer->ReadTableDef(psTable))
1610 : {
1611 0 : return false;
1612 : }
1613 22 : if (eAccess == GA_Update)
1614 : {
1615 6 : m_apoLayers.push_back(
1616 12 : std::make_unique<PDS4EditableLayer>(std::move(poLayer)));
1617 : }
1618 : else
1619 : {
1620 16 : m_apoLayers.push_back(std::move(poLayer));
1621 : }
1622 22 : return true;
1623 : }
1624 :
1625 : /************************************************************************/
1626 : /* OpenTableBinary() */
1627 : /************************************************************************/
1628 :
1629 11 : bool PDS4Dataset::OpenTableBinary(const char *pszFilename,
1630 : const CPLXMLNode *psTable)
1631 : {
1632 11 : if (CPLHasPathTraversal(pszFilename))
1633 : {
1634 0 : CPLError(CE_Failure, CPLE_NotSupported,
1635 : "OpenTableBinary(): path traversal detected: %s", pszFilename);
1636 0 : return false;
1637 : }
1638 :
1639 22 : std::string osLayerName(CPLGetBasenameSafe(pszFilename));
1640 11 : if (cpl::starts_with(osLayerName,
1641 22 : CPLGetBasenameSafe(m_osXMLFilename.c_str()) + "_"))
1642 20 : osLayerName = osLayerName.substr(
1643 30 : CPLGetBasenameSafe(m_osXMLFilename.c_str()).size() + 1);
1644 :
1645 22 : const std::string osFullFilename = FixupTableFilename(CPLFormFilenameSafe(
1646 33 : CPLGetPathSafe(m_osXMLFilename.c_str()).c_str(), pszFilename, nullptr));
1647 0 : auto poLayer = std::make_unique<PDS4TableBinary>(this, osLayerName.c_str(),
1648 11 : osFullFilename.c_str(),
1649 33 : eAccess == GA_Update);
1650 11 : if (!poLayer->ReadTableDef(psTable))
1651 : {
1652 0 : return false;
1653 : }
1654 11 : if (eAccess == GA_Update)
1655 : {
1656 1 : m_apoLayers.push_back(
1657 2 : std::make_unique<PDS4EditableLayer>(std::move(poLayer)));
1658 : }
1659 : else
1660 : {
1661 10 : m_apoLayers.push_back(std::move(poLayer));
1662 : }
1663 11 : return true;
1664 : }
1665 :
1666 : /************************************************************************/
1667 : /* OpenTableDelimited() */
1668 : /************************************************************************/
1669 :
1670 64 : bool PDS4Dataset::OpenTableDelimited(const char *pszFilename,
1671 : const CPLXMLNode *psTable)
1672 : {
1673 64 : if (CPLHasPathTraversal(pszFilename))
1674 : {
1675 0 : CPLError(CE_Failure, CPLE_NotSupported,
1676 : "OpenTableDelimited(): path traversal detected: %s",
1677 : pszFilename);
1678 0 : return false;
1679 : }
1680 :
1681 128 : std::string osLayerName(CPLGetBasenameSafe(pszFilename));
1682 64 : if (cpl::starts_with(osLayerName,
1683 128 : CPLGetBasenameSafe(m_osXMLFilename.c_str()) + "_"))
1684 120 : osLayerName = osLayerName.substr(
1685 180 : CPLGetBasenameSafe(m_osXMLFilename.c_str()).size() + 1);
1686 :
1687 128 : const std::string osFullFilename = FixupTableFilename(CPLFormFilenameSafe(
1688 192 : CPLGetPathSafe(m_osXMLFilename.c_str()).c_str(), pszFilename, nullptr));
1689 : auto poLayer = std::make_unique<PDS4DelimitedTable>(
1690 64 : this, osLayerName.c_str(), osFullFilename.c_str(),
1691 192 : eAccess == GA_Update);
1692 64 : if (!poLayer->ReadTableDef(psTable))
1693 : {
1694 0 : return false;
1695 : }
1696 64 : if (eAccess == GA_Update)
1697 : {
1698 2 : m_apoLayers.push_back(
1699 4 : std::make_unique<PDS4EditableLayer>(std::move(poLayer)));
1700 : }
1701 : else
1702 : {
1703 62 : m_apoLayers.push_back(std::move(poLayer));
1704 : }
1705 64 : return true;
1706 : }
1707 :
1708 : /************************************************************************/
1709 : /* ConstantToDouble() */
1710 : /************************************************************************/
1711 :
1712 180 : static std::optional<double> ConstantToDouble(const char *pszItem,
1713 : const char *pszVal)
1714 : {
1715 180 : if (STARTS_WITH(pszVal, "0x"))
1716 : {
1717 21 : if (strlen(pszVal) == strlen("0x") + 2 * sizeof(float))
1718 : {
1719 15 : char *endptr = nullptr;
1720 : const uint32_t nVal =
1721 15 : static_cast<uint32_t>(std::strtoull(pszVal, &endptr, 0));
1722 15 : if (endptr == pszVal + strlen(pszVal))
1723 : {
1724 : float fVal;
1725 15 : memcpy(&fVal, &nVal, sizeof(nVal));
1726 15 : return fVal;
1727 : }
1728 : }
1729 6 : else if (strlen(pszVal) == strlen("0x") + 2 * sizeof(double))
1730 : {
1731 6 : char *endptr = nullptr;
1732 : const uint64_t nVal =
1733 6 : static_cast<uint64_t>(std::strtoull(pszVal, &endptr, 0));
1734 6 : if (endptr == pszVal + strlen(pszVal))
1735 : {
1736 : double dfVal;
1737 6 : memcpy(&dfVal, &nVal, sizeof(nVal));
1738 6 : return dfVal;
1739 : }
1740 : }
1741 0 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid value for '%s': '%s'",
1742 : pszItem, pszVal);
1743 0 : return std::nullopt;
1744 : }
1745 : else
1746 : {
1747 159 : char *endptr = nullptr;
1748 159 : double dfVal = std::strtod(pszVal, &endptr);
1749 159 : if (endptr == pszVal + strlen(pszVal))
1750 : {
1751 159 : return dfVal;
1752 : }
1753 : else
1754 : {
1755 0 : CPLError(CE_Failure, CPLE_AppDefined,
1756 : "Invalid value for '%s': '%s'", pszItem, pszVal);
1757 0 : return std::nullopt;
1758 : }
1759 : }
1760 : }
1761 :
1762 : /************************************************************************/
1763 : /* OpenBrowse() */
1764 : /************************************************************************/
1765 :
1766 : /* static */ std::unique_ptr<PDS4Dataset>
1767 2 : PDS4Dataset::OpenBrowse(GDALOpenInfo *poOpenInfo, const CPLXMLNode *psProduct)
1768 : {
1769 : const CPLXMLNode *psFileAreaBrowse =
1770 2 : CPLGetXMLNode(psProduct, "File_Area_Browse");
1771 2 : if (!psFileAreaBrowse)
1772 0 : return nullptr;
1773 : const char *pszFilename =
1774 2 : CPLGetXMLValue(psFileAreaBrowse, "File.file_name", nullptr);
1775 2 : if (!pszFilename)
1776 0 : return nullptr;
1777 2 : if (CPLHasPathTraversal(pszFilename))
1778 : {
1779 0 : CPLError(CE_Failure, CPLE_NotSupported, "Path traversal detected in %s",
1780 : pszFilename);
1781 0 : return nullptr;
1782 : }
1783 2 : const char *pszFormat = CPLGetXMLValue(
1784 : psFileAreaBrowse, "Encoded_Image.encoding_standard_id", "");
1785 3 : const char *const apszAllowedDrivers[] = {EQUAL(pszFormat, "TIFF") ? "GTiff"
1786 1 : : EQUAL(pszFormat, "PNG")
1787 1 : ? "PNG"
1788 : : nullptr,
1789 2 : nullptr};
1790 2 : if (apszAllowedDrivers[0] == nullptr)
1791 : {
1792 0 : CPLError(CE_Failure, CPLE_NotSupported,
1793 : "Unhandled encoding_standard_id=%s", pszFormat);
1794 0 : return nullptr;
1795 : }
1796 : const std::string osImageFullFilename = CPLFormFilenameSafe(
1797 4 : CPLGetPathSafe(poOpenInfo->pszFilename).c_str(), pszFilename, nullptr);
1798 : auto poExternalDS = std::unique_ptr<GDALDataset>(GDALDataset::Open(
1799 : osImageFullFilename.c_str(), GDAL_OF_RASTER | GDAL_OF_VERBOSE_ERROR,
1800 4 : apszAllowedDrivers));
1801 2 : if (!poExternalDS)
1802 0 : return nullptr;
1803 :
1804 4 : auto poDS = std::make_unique<PDS4Dataset>();
1805 2 : poDS->SetDescription(poOpenInfo->pszFilename);
1806 2 : poDS->m_poExternalDS = poExternalDS.release();
1807 2 : poDS->nRasterXSize = poDS->m_poExternalDS->GetRasterXSize();
1808 2 : poDS->nRasterYSize = poDS->m_poExternalDS->GetRasterYSize();
1809 2 : poDS->eAccess = GA_ReadOnly;
1810 2 : poDS->m_osImageFilename = osImageFullFilename;
1811 :
1812 4 : for (int i = 0; i < poDS->m_poExternalDS->GetRasterCount(); i++)
1813 : {
1814 : auto poBand = std::make_unique<PDS4BrowseImageProxyRasterBand>(
1815 2 : poDS->m_poExternalDS->GetRasterBand(i + 1));
1816 2 : poDS->SetBand(i + 1, std::move(poBand));
1817 : }
1818 :
1819 2 : if (poDS->m_poExternalDS->GetGeoTransform(poDS->m_gt) == CE_None)
1820 : {
1821 1 : poDS->m_bGotTransform = true;
1822 : }
1823 2 : const auto poSRS = poDS->m_poExternalDS->GetSpatialRef();
1824 2 : if (poSRS)
1825 1 : poDS->m_oSRS = *poSRS;
1826 :
1827 6 : for (const char *pszDomain : {GDAL_MDD_DEFAULT, GDAL_MDD_IMAGE_STRUCTURE})
1828 : {
1829 8 : poDS->GDALDataset::SetMetadata(
1830 4 : poDS->m_poExternalDS->GetMetadata(pszDomain), pszDomain);
1831 : }
1832 :
1833 : // Expose XML content in xml:PDS4 metadata domain
1834 2 : GByte *pabyRet = nullptr;
1835 2 : CPL_IGNORE_RET_VAL(VSIIngestFile(nullptr, poOpenInfo->pszFilename, &pabyRet,
1836 : nullptr, 10 * 1024 * 1024));
1837 2 : if (pabyRet)
1838 : {
1839 2 : char *apszXML[2] = {reinterpret_cast<char *>(pabyRet), nullptr};
1840 2 : poDS->GDALDataset::SetMetadata(apszXML, "xml:PDS4");
1841 : }
1842 2 : VSIFree(pabyRet);
1843 :
1844 : /*--------------------------------------------------------------------------*/
1845 : /* Initialize any PAM information */
1846 : /*--------------------------------------------------------------------------*/
1847 2 : poDS->SetDescription(poOpenInfo->pszFilename);
1848 2 : poDS->TryLoadXML();
1849 :
1850 2 : return poDS;
1851 : }
1852 :
1853 : /************************************************************************/
1854 : /* Open() */
1855 : /************************************************************************/
1856 :
1857 : // See https://pds.nasa.gov/pds4/pds/v1/PDS4_PDS_1800.xsd
1858 : // and https://pds.nasa.gov/pds4/pds/v1/PDS4_PDS_1800.sch
1859 308 : std::unique_ptr<PDS4Dataset> PDS4Dataset::OpenInternal(GDALOpenInfo *poOpenInfo)
1860 : {
1861 308 : if (!PDS4DriverIdentify(poOpenInfo))
1862 0 : return nullptr;
1863 :
1864 616 : CPLString osXMLFilename(poOpenInfo->pszFilename);
1865 308 : int nFAOIdxLookup = -1;
1866 308 : int nArrayIdxLookup = -1;
1867 308 : if (STARTS_WITH_CI(poOpenInfo->pszFilename, "PDS4:"))
1868 : {
1869 : char **papszTokens =
1870 15 : CSLTokenizeString2(poOpenInfo->pszFilename, ":", 0);
1871 15 : int nCount = CSLCount(papszTokens);
1872 15 : if (nCount == 5 && strlen(papszTokens[1]) == 1 &&
1873 1 : (papszTokens[2][0] == '\\' || papszTokens[2][0] == '/'))
1874 : {
1875 1 : osXMLFilename = CPLString(papszTokens[1]) + ":" + papszTokens[2];
1876 1 : nFAOIdxLookup = atoi(papszTokens[3]);
1877 1 : nArrayIdxLookup = atoi(papszTokens[4]);
1878 : }
1879 14 : else if (nCount == 5 && (EQUAL(papszTokens[1], "/vsicurl/http") ||
1880 0 : EQUAL(papszTokens[1], "/vsicurl/https")))
1881 : {
1882 0 : osXMLFilename = CPLString(papszTokens[1]) + ":" + papszTokens[2];
1883 0 : nFAOIdxLookup = atoi(papszTokens[3]);
1884 0 : nArrayIdxLookup = atoi(papszTokens[4]);
1885 : }
1886 14 : else if (nCount == 4)
1887 : {
1888 13 : osXMLFilename = papszTokens[1];
1889 13 : nFAOIdxLookup = atoi(papszTokens[2]);
1890 13 : nArrayIdxLookup = atoi(papszTokens[3]);
1891 : }
1892 : else
1893 : {
1894 1 : CPLError(CE_Failure, CPLE_AppDefined,
1895 : "Invalid syntax for PDS4 subdataset name");
1896 1 : CSLDestroy(papszTokens);
1897 1 : return nullptr;
1898 : }
1899 14 : CSLDestroy(papszTokens);
1900 : }
1901 :
1902 614 : CPLXMLTreeCloser oCloser(CPLParseXMLFile(osXMLFilename));
1903 307 : CPLXMLNode *psRoot = oCloser.get();
1904 307 : CPLStripXMLNamespace(psRoot, nullptr, TRUE);
1905 :
1906 614 : GDALAccess eAccess = STARTS_WITH_CI(poOpenInfo->pszFilename, "PDS4:")
1907 307 : ? GA_ReadOnly
1908 : : poOpenInfo->eAccess;
1909 :
1910 307 : CPLXMLNode *psProduct = CPLGetXMLNode(psRoot, "=Product_Observational");
1911 307 : if (psProduct == nullptr)
1912 : {
1913 6 : eAccess = GA_ReadOnly;
1914 6 : psProduct = CPLGetXMLNode(psRoot, "=Product_Ancillary");
1915 6 : if (psProduct == nullptr)
1916 : {
1917 6 : psProduct = CPLGetXMLNode(psRoot, "=Product_Collection");
1918 : }
1919 6 : if (psProduct == nullptr)
1920 : {
1921 5 : psProduct = CPLGetXMLNode(psRoot, "=Product_Browse");
1922 5 : if (psProduct)
1923 2 : return OpenBrowse(poOpenInfo, psProduct);
1924 : }
1925 : }
1926 305 : if (psProduct == nullptr)
1927 : {
1928 3 : return nullptr;
1929 : }
1930 :
1931 : // Test case:
1932 : // https://starbase.jpl.nasa.gov/pds4/1700/dph_example_products/test_Images_DisplaySettings/TestPattern_Image/TestPattern.xml
1933 302 : const char *pszVertDir = CPLGetXMLValue(
1934 : psProduct,
1935 : "Observation_Area.Discipline_Area.Display_Settings.Display_Direction."
1936 : "vertical_display_direction",
1937 : "");
1938 302 : const bool bBottomToTop = EQUAL(pszVertDir, "Bottom to Top");
1939 :
1940 302 : const char *pszHorizDir = CPLGetXMLValue(
1941 : psProduct,
1942 : "Observation_Area.Discipline_Area.Display_Settings.Display_Direction."
1943 : "horizontal_display_direction",
1944 : "");
1945 302 : const bool bRightToLeft = EQUAL(pszHorizDir, "Right to Left");
1946 :
1947 604 : auto poDS = std::make_unique<PDS4Dataset>();
1948 302 : poDS->m_osXMLFilename = osXMLFilename;
1949 302 : poDS->eAccess = eAccess;
1950 302 : poDS->papszOpenOptions = CSLDuplicate(poOpenInfo->papszOpenOptions);
1951 :
1952 604 : CPLStringList aosSubdatasets;
1953 302 : int nFAOIdx = 0;
1954 2946 : for (CPLXMLNode *psIter = psProduct->psChild; psIter != nullptr;
1955 2644 : psIter = psIter->psNext)
1956 : {
1957 2646 : if (psIter->eType != CXT_Element ||
1958 933 : (strcmp(psIter->pszValue, "File_Area_Observational") != 0 &&
1959 602 : strcmp(psIter->pszValue, "File_Area_Ancillary") != 0 &&
1960 602 : strcmp(psIter->pszValue, "File_Area_Inventory") != 0))
1961 : {
1962 2314 : continue;
1963 : }
1964 :
1965 332 : nFAOIdx++;
1966 332 : CPLXMLNode *psFile = CPLGetXMLNode(psIter, "File");
1967 332 : if (psFile == nullptr)
1968 : {
1969 1 : continue;
1970 : }
1971 331 : const char *pszFilename = CPLGetXMLValue(psFile, "file_name", nullptr);
1972 331 : if (pszFilename == nullptr)
1973 : {
1974 1 : continue;
1975 : }
1976 :
1977 693 : for (CPLXMLNode *psSubIter = psFile->psChild; psSubIter != nullptr;
1978 363 : psSubIter = psSubIter->psNext)
1979 : {
1980 363 : if (psSubIter->eType == CXT_Comment &&
1981 14 : EQUAL(psSubIter->pszValue, PREEXISTING_BINARY_FILE))
1982 : {
1983 14 : poDS->m_bCreatedFromExistingBinaryFile = true;
1984 : }
1985 : }
1986 :
1987 330 : int nArrayIdx = 0;
1988 330 : for (CPLXMLNode *psSubIter = psIter->psChild;
1989 1050 : (nFAOIdxLookup < 0 || nFAOIdxLookup == nFAOIdx) &&
1990 : psSubIter != nullptr;
1991 720 : psSubIter = psSubIter->psNext)
1992 : {
1993 722 : if (psSubIter->eType != CXT_Element)
1994 : {
1995 504 : continue;
1996 : }
1997 722 : int nDIM = 0;
1998 722 : if (STARTS_WITH(psSubIter->pszValue, "Array_1D"))
1999 : {
2000 0 : nDIM = 1;
2001 : }
2002 722 : else if (STARTS_WITH(psSubIter->pszValue, "Array_2D"))
2003 : {
2004 6 : nDIM = 2;
2005 : }
2006 716 : else if (STARTS_WITH(psSubIter->pszValue, "Array_3D"))
2007 : {
2008 241 : nDIM = 3;
2009 : }
2010 475 : else if (strcmp(psSubIter->pszValue, "Array") == 0)
2011 : {
2012 3 : nDIM = atoi(CPLGetXMLValue(psSubIter, "axes", "0"));
2013 : }
2014 472 : else if (strcmp(psSubIter->pszValue, "Table_Character") == 0)
2015 : {
2016 22 : poDS->OpenTableCharacter(pszFilename, psSubIter);
2017 22 : continue;
2018 : }
2019 450 : else if (strcmp(psSubIter->pszValue, "Table_Binary") == 0)
2020 : {
2021 11 : poDS->OpenTableBinary(pszFilename, psSubIter);
2022 11 : continue;
2023 : }
2024 439 : else if (strcmp(psSubIter->pszValue, "Table_Delimited") == 0 ||
2025 376 : strcmp(psSubIter->pszValue, "Inventory") == 0)
2026 : {
2027 64 : poDS->OpenTableDelimited(pszFilename, psSubIter);
2028 64 : continue;
2029 : }
2030 625 : if (nDIM == 0)
2031 : {
2032 375 : continue;
2033 : }
2034 250 : if (!(nDIM >= 1 && nDIM <= 3))
2035 : {
2036 0 : CPLError(CE_Warning, CPLE_NotSupported,
2037 : "Array with %d dimensions not supported", nDIM);
2038 0 : continue;
2039 : }
2040 :
2041 250 : nArrayIdx++;
2042 : // Does it match a selected subdataset ?
2043 250 : if (nArrayIdxLookup > 0 && nArrayIdx != nArrayIdxLookup)
2044 : {
2045 13 : continue;
2046 : }
2047 :
2048 : const char *pszArrayName =
2049 237 : CPLGetXMLValue(psSubIter, "name", nullptr);
2050 : const char *pszArrayId =
2051 237 : CPLGetXMLValue(psSubIter, "local_identifier", nullptr);
2052 : vsi_l_offset nOffset = static_cast<vsi_l_offset>(
2053 237 : CPLAtoGIntBig(CPLGetXMLValue(psSubIter, "offset", "0")));
2054 :
2055 : const char *pszAxisIndexOrder =
2056 237 : CPLGetXMLValue(psSubIter, "axis_index_order", "");
2057 237 : if (!EQUAL(pszAxisIndexOrder, "Last Index Fastest"))
2058 : {
2059 1 : CPLError(CE_Warning, CPLE_NotSupported,
2060 : "axis_index_order = '%s' unhandled",
2061 : pszAxisIndexOrder);
2062 1 : continue;
2063 : }
2064 :
2065 : // Figure out data type
2066 : const char *pszDataType =
2067 236 : CPLGetXMLValue(psSubIter, "Element_Array.data_type", "");
2068 236 : GDALDataType eDT = GDT_UInt8;
2069 236 : bool bLSBOrder = strstr(pszDataType, "LSB") != nullptr;
2070 :
2071 : // ComplexLSB16', 'ComplexLSB8', 'ComplexMSB16', 'ComplexMSB8',
2072 : // 'IEEE754LSBDouble', 'IEEE754LSBSingle', 'IEEE754MSBDouble',
2073 : // 'IEEE754MSBSingle', 'SignedBitString', 'SignedByte',
2074 : // 'SignedLSB2', 'SignedLSB4', 'SignedLSB8', 'SignedMSB2',
2075 : // 'SignedMSB4', 'SignedMSB8', 'UnsignedBitString', 'UnsignedByte',
2076 : // 'UnsignedLSB2', 'UnsignedLSB4', 'UnsignedLSB8', 'UnsignedMSB2',
2077 : // 'UnsignedMSB4', 'UnsignedMSB8'
2078 236 : if (EQUAL(pszDataType, "ComplexLSB16") ||
2079 232 : EQUAL(pszDataType, "ComplexMSB16"))
2080 : {
2081 6 : eDT = GDT_CFloat64;
2082 : }
2083 230 : else if (EQUAL(pszDataType, "ComplexLSB8") ||
2084 226 : EQUAL(pszDataType, "ComplexMSB8"))
2085 : {
2086 4 : eDT = GDT_CFloat32;
2087 : }
2088 226 : else if (EQUAL(pszDataType, "IEEE754LSBDouble") ||
2089 219 : EQUAL(pszDataType, "IEEE754MSBDouble"))
2090 : {
2091 7 : eDT = GDT_Float64;
2092 : }
2093 219 : else if (EQUAL(pszDataType, "IEEE754LSBSingle") ||
2094 208 : EQUAL(pszDataType, "IEEE754MSBSingle"))
2095 : {
2096 12 : eDT = GDT_Float32;
2097 : }
2098 : // SignedBitString unhandled
2099 207 : else if (EQUAL(pszDataType, "SignedByte"))
2100 : {
2101 6 : eDT = GDT_Int8;
2102 : }
2103 201 : else if (EQUAL(pszDataType, "SignedLSB2") ||
2104 197 : EQUAL(pszDataType, "SignedMSB2"))
2105 : {
2106 6 : eDT = GDT_Int16;
2107 : }
2108 195 : else if (EQUAL(pszDataType, "SignedLSB4") ||
2109 191 : EQUAL(pszDataType, "SignedMSB4"))
2110 : {
2111 4 : eDT = GDT_Int32;
2112 : }
2113 191 : else if (EQUAL(pszDataType, "SignedLSB8") ||
2114 187 : EQUAL(pszDataType, "SignedMSB8"))
2115 : {
2116 4 : eDT = GDT_Int64;
2117 : }
2118 187 : else if (EQUAL(pszDataType, "UnsignedByte"))
2119 : {
2120 170 : eDT = GDT_UInt8;
2121 : }
2122 17 : else if (EQUAL(pszDataType, "UnsignedLSB2") ||
2123 9 : EQUAL(pszDataType, "UnsignedMSB2"))
2124 : {
2125 8 : eDT = GDT_UInt16;
2126 : }
2127 9 : else if (EQUAL(pszDataType, "UnsignedLSB4") ||
2128 5 : EQUAL(pszDataType, "UnsignedMSB4"))
2129 : {
2130 4 : eDT = GDT_UInt32;
2131 : }
2132 5 : else if (EQUAL(pszDataType, "UnsignedLSB8") ||
2133 1 : EQUAL(pszDataType, "UnsignedMSB8"))
2134 : {
2135 4 : eDT = GDT_UInt64;
2136 : }
2137 : else
2138 : {
2139 1 : CPLDebug("PDS4", "data_type = '%s' unhandled", pszDataType);
2140 1 : continue;
2141 : }
2142 :
2143 235 : poDS->m_osUnits =
2144 470 : CPLGetXMLValue(psSubIter, "Element_Array.unit", "");
2145 :
2146 235 : double dfValueOffset = CPLAtof(
2147 : CPLGetXMLValue(psSubIter, "Element_Array.value_offset", "0"));
2148 235 : double dfValueScale = CPLAtof(
2149 : CPLGetXMLValue(psSubIter, "Element_Array.scaling_factor", "1"));
2150 :
2151 : // Parse Axis_Array elements
2152 235 : int l_nBands = 1;
2153 235 : int nLines = 0;
2154 235 : int nSamples = 0;
2155 235 : std::vector<int> anElements;
2156 235 : std::vector<std::string> axisNames;
2157 235 : std::vector<char> dimSemantics;
2158 235 : anElements.resize(nDIM);
2159 235 : axisNames.resize(nDIM);
2160 235 : dimSemantics.resize(nDIM);
2161 235 : int iBandIdx = -1;
2162 235 : int iLineIdx = -1;
2163 235 : int iSampleIdx = -1;
2164 235 : int nAxisOKCount = 0;
2165 235 : for (CPLXMLNode *psAxisIter = psSubIter->psChild;
2166 2161 : psAxisIter != nullptr; psAxisIter = psAxisIter->psNext)
2167 : {
2168 1927 : if (psAxisIter->eType != CXT_Element ||
2169 1927 : strcmp(psAxisIter->pszValue, "Axis_Array") != 0)
2170 : {
2171 1224 : continue;
2172 : }
2173 : const char *pszAxisName =
2174 703 : CPLGetXMLValue(psAxisIter, "axis_name", nullptr);
2175 : const char *pszElements =
2176 703 : CPLGetXMLValue(psAxisIter, "elements", nullptr);
2177 : const char *pszSequenceNumber =
2178 703 : CPLGetXMLValue(psAxisIter, "sequence_number", nullptr);
2179 703 : if (pszAxisName == nullptr || pszElements == nullptr ||
2180 : pszSequenceNumber == nullptr)
2181 : {
2182 1 : continue;
2183 : }
2184 702 : int nSeqNumber = atoi(pszSequenceNumber);
2185 702 : if (nSeqNumber < 1 || nSeqNumber > nDIM)
2186 : {
2187 2 : CPLError(CE_Warning, CPLE_AppDefined,
2188 : "Invalid sequence_number = %s", pszSequenceNumber);
2189 2 : continue;
2190 : }
2191 700 : int nElements = atoi(pszElements);
2192 700 : if (nElements <= 0)
2193 : {
2194 1 : CPLError(CE_Warning, CPLE_AppDefined,
2195 : "Invalid elements = %s", pszElements);
2196 1 : continue;
2197 : }
2198 :
2199 699 : const int nIdx = nSeqNumber - 1;
2200 699 : if (STARTS_WITH_CI(pszAxisName, "Line"))
2201 : {
2202 234 : if (iLineIdx < 0)
2203 234 : iLineIdx = nIdx;
2204 : else
2205 : {
2206 0 : CPLError(CE_Warning, CPLE_AppDefined,
2207 : "Several axis with Line identifier");
2208 0 : break;
2209 : }
2210 : }
2211 465 : else if (STARTS_WITH_CI(pszAxisName, "Sample"))
2212 : {
2213 234 : if (iSampleIdx < 0)
2214 234 : iSampleIdx = nIdx;
2215 : else
2216 : {
2217 0 : CPLError(CE_Warning, CPLE_AppDefined,
2218 : "Several axis with Sample identifier");
2219 0 : break;
2220 : }
2221 : }
2222 231 : else if (STARTS_WITH_CI(pszAxisName, "Band"))
2223 : {
2224 229 : if (iBandIdx < 0)
2225 228 : iBandIdx = nIdx;
2226 : else
2227 : {
2228 1 : CPLError(CE_Warning, CPLE_AppDefined,
2229 : "Several axis with Band identifier");
2230 1 : break;
2231 : }
2232 : }
2233 :
2234 698 : anElements[nIdx] = nElements;
2235 698 : axisNames[nIdx] = pszAxisName;
2236 :
2237 698 : ++nAxisOKCount;
2238 : }
2239 :
2240 235 : if (nAxisOKCount != nDIM)
2241 : {
2242 1 : CPLError(CE_Warning, CPLE_AppDefined,
2243 : "Found only %d Axis_Array elements. %d expected",
2244 : nAxisOKCount, nDIM);
2245 1 : continue;
2246 : }
2247 :
2248 234 : if (nDIM == 1)
2249 : {
2250 0 : dimSemantics[0] = 'S';
2251 0 : nLines = 1;
2252 0 : nSamples = anElements[0];
2253 : }
2254 234 : else if (nDIM == 2)
2255 : {
2256 6 : if (iLineIdx < 0 || iSampleIdx < 0)
2257 : {
2258 0 : CPLDebug("PDS4", "Assume that axis %s is Line",
2259 0 : axisNames[0].c_str());
2260 0 : CPLDebug("PDS4", "Assume that axis %s is Sample",
2261 0 : axisNames[1].c_str());
2262 0 : iLineIdx = 0;
2263 0 : iSampleIdx = 1;
2264 : }
2265 6 : CPLAssert(iLineIdx >= 0 && iLineIdx < 2);
2266 6 : CPLAssert(iSampleIdx >= 0 && iSampleIdx < 2);
2267 6 : dimSemantics[iLineIdx] = 'L';
2268 6 : dimSemantics[iSampleIdx] = 'S';
2269 6 : nLines = anElements[iLineIdx];
2270 6 : nSamples = anElements[iSampleIdx];
2271 : }
2272 : else /* if (nDim == 3) */
2273 : {
2274 228 : if (iLineIdx < 0 || iSampleIdx < 0)
2275 : {
2276 0 : CPLDebug("PDS4", "Assume that axis %s is Band",
2277 0 : axisNames[0].c_str());
2278 0 : CPLDebug("PDS4", "Assume that axis %s is Line",
2279 0 : axisNames[1].c_str());
2280 0 : CPLDebug("PDS4", "Assume that axis %s is Sample",
2281 0 : axisNames[2].c_str());
2282 0 : iBandIdx = 0;
2283 0 : iLineIdx = 1;
2284 0 : iSampleIdx = 2;
2285 : }
2286 228 : else if (iBandIdx < 0)
2287 : {
2288 1 : CPLAssert(iLineIdx >= 0 && iLineIdx < 3);
2289 1 : CPLAssert(iSampleIdx >= 0 && iSampleIdx < 3);
2290 1 : bool abUsedIndices[3] = {false, false, false};
2291 1 : abUsedIndices[iLineIdx] = true;
2292 1 : abUsedIndices[iSampleIdx] = true;
2293 4 : for (int i = 0; i < 3; ++i)
2294 : {
2295 3 : if (!abUsedIndices[i])
2296 1 : iBandIdx = i;
2297 : }
2298 : }
2299 228 : CPLAssert(iLineIdx >= 0 && iLineIdx < 3);
2300 228 : CPLAssert(iSampleIdx >= 0 && iSampleIdx < 3);
2301 228 : CPLAssert(iBandIdx >= 0 && iSampleIdx < 3);
2302 228 : dimSemantics[iBandIdx] = 'B';
2303 228 : dimSemantics[iLineIdx] = 'L';
2304 228 : dimSemantics[iSampleIdx] = 'S';
2305 228 : l_nBands = anElements[iBandIdx];
2306 228 : nLines = anElements[iLineIdx];
2307 228 : nSamples = anElements[iSampleIdx];
2308 : }
2309 :
2310 468 : if (!GDALCheckDatasetDimensions(nSamples, nLines) ||
2311 234 : !GDALCheckBandCount(l_nBands, FALSE))
2312 : {
2313 1 : continue;
2314 : }
2315 :
2316 : // Compute pixel, line and band spacing
2317 233 : vsi_l_offset nSpacing = GDALGetDataTypeSizeBytes(eDT);
2318 233 : int nPixelOffset = 0;
2319 233 : int nLineOffset = 0;
2320 233 : vsi_l_offset nBandOffset = 0;
2321 233 : int nCountPreviousDim = 1;
2322 924 : for (int i = nDIM - 1; i >= 0; i--)
2323 : {
2324 693 : if (dimSemantics[i] == 'S')
2325 : {
2326 233 : if (nCountPreviousDim > 0 &&
2327 233 : nSpacing > static_cast<vsi_l_offset>(INT_MAX /
2328 : nCountPreviousDim))
2329 : {
2330 1 : CPLError(CE_Failure, CPLE_NotSupported,
2331 : "Integer overflow");
2332 1 : return nullptr;
2333 : }
2334 232 : nPixelOffset =
2335 232 : static_cast<int>(nSpacing * nCountPreviousDim);
2336 232 : nSpacing = nPixelOffset;
2337 : }
2338 460 : else if (dimSemantics[i] == 'L')
2339 : {
2340 233 : if (nCountPreviousDim > 0 &&
2341 233 : nSpacing > static_cast<vsi_l_offset>(INT_MAX /
2342 : nCountPreviousDim))
2343 : {
2344 1 : CPLError(CE_Failure, CPLE_NotSupported,
2345 : "Integer overflow");
2346 1 : return nullptr;
2347 : }
2348 232 : nLineOffset =
2349 232 : static_cast<int>(nSpacing * nCountPreviousDim);
2350 232 : nSpacing = nLineOffset;
2351 : }
2352 : else
2353 : {
2354 227 : nBandOffset = nSpacing * nCountPreviousDim;
2355 227 : nSpacing = nBandOffset;
2356 : }
2357 691 : nCountPreviousDim = anElements[i];
2358 : }
2359 :
2360 : // Retrieve no data value
2361 231 : bool bNoDataSet = false;
2362 231 : double dfNoData = 0.0;
2363 231 : bool bNoDataSetInt64 = false;
2364 231 : int64_t nNoDataInt64 = 0;
2365 231 : bool bNoDataSetUInt64 = false;
2366 231 : int64_t nNoDataUInt64 = 0;
2367 231 : std::vector<double> adfConstants;
2368 231 : CPLXMLNode *psSC = CPLGetXMLNode(psSubIter, "Special_Constants");
2369 231 : if (psSC)
2370 : {
2371 : const auto GetNoDataFromConstant =
2372 76 : [psSC, eDT, &bNoDataSetInt64, &nNoDataInt64,
2373 : &bNoDataSetUInt64, &nNoDataUInt64, &bNoDataSet,
2374 372 : &dfNoData](const char *pszConstantName)
2375 : {
2376 76 : if (const char *pszVal =
2377 76 : CPLGetXMLValue(psSC, pszConstantName, nullptr))
2378 : {
2379 75 : if (eDT == GDT_Int64)
2380 : {
2381 2 : bNoDataSetInt64 = true;
2382 2 : nNoDataInt64 = std::strtoll(pszVal, nullptr, 10);
2383 : }
2384 73 : else if (eDT == GDT_UInt64)
2385 : {
2386 2 : bNoDataSetUInt64 = true;
2387 2 : nNoDataUInt64 = std::strtoull(pszVal, nullptr, 10);
2388 : }
2389 : else
2390 : {
2391 : auto val =
2392 71 : ConstantToDouble(pszConstantName, pszVal);
2393 71 : if (val)
2394 : {
2395 71 : bNoDataSet = true;
2396 71 : dfNoData = *val;
2397 : }
2398 : }
2399 75 : return true;
2400 : }
2401 1 : return false;
2402 75 : };
2403 :
2404 75 : if (!GetNoDataFromConstant("missing_constant"))
2405 : {
2406 : // For example in https://d34uoeqxvp7znu.cloudfront.net/data/l2/200908/20090811/m3g20090811t012112_l2.xml
2407 1 : GetNoDataFromConstant("invalid_constant");
2408 : }
2409 :
2410 75 : const char *const apszConstantNames[] = {
2411 : "saturated_constant",
2412 : "missing_constant",
2413 : "error_constant",
2414 : "invalid_constant",
2415 : "unknown_constant",
2416 : "not_applicable_constant",
2417 : "high_instrument_saturation",
2418 : "high_representation_saturation",
2419 : "low_instrument_saturation",
2420 : "low_representation_saturation"};
2421 825 : for (const char *pszItem : apszConstantNames)
2422 : {
2423 750 : if (const char *pszConstant =
2424 750 : CPLGetXMLValue(psSC, pszItem, nullptr))
2425 : {
2426 109 : auto val = ConstantToDouble(pszItem, pszConstant);
2427 109 : if (val)
2428 : {
2429 109 : adfConstants.push_back(*val);
2430 : }
2431 : }
2432 : }
2433 : }
2434 :
2435 : // Add subdatasets
2436 231 : const int nSDSIdx = 1 + aosSubdatasets.size() / 2;
2437 : aosSubdatasets.SetNameValue(
2438 : CPLSPrintf("SUBDATASET_%d_NAME", nSDSIdx),
2439 : CPLSPrintf("PDS4:%s:%d:%d", osXMLFilename.c_str(), nFAOIdx,
2440 231 : nArrayIdx));
2441 : aosSubdatasets.SetNameValue(
2442 : CPLSPrintf("SUBDATASET_%d_DESC", nSDSIdx),
2443 : CPLSPrintf("Image file %s, array %s", pszFilename,
2444 : pszArrayName ? pszArrayName
2445 231 : : pszArrayId ? pszArrayId
2446 462 : : CPLSPrintf("%d", nArrayIdx)));
2447 :
2448 231 : if (poDS->nBands != 0)
2449 14 : continue;
2450 :
2451 : const std::string osImageFullFilename = CPLFormFilenameSafe(
2452 217 : CPLGetPathSafe(osXMLFilename.c_str()).c_str(), pszFilename,
2453 217 : nullptr);
2454 217 : if (CPLHasPathTraversal(pszFilename))
2455 : {
2456 0 : CPLError(CE_Failure, CPLE_NotSupported,
2457 : "Path traversal detected in %s", pszFilename);
2458 0 : return nullptr;
2459 : }
2460 217 : VSILFILE *fp = VSIFOpenExL(
2461 : osImageFullFilename.c_str(),
2462 217 : (poOpenInfo->eAccess == GA_Update) ? "rb+" : "rb", true);
2463 217 : if (fp == nullptr)
2464 : {
2465 1 : CPLError(CE_Warning, CPLE_FileIO, "Cannot open %s: %s",
2466 : osImageFullFilename.c_str(), VSIGetLastErrorMsg());
2467 1 : continue;
2468 : }
2469 216 : poDS->nRasterXSize = nSamples;
2470 216 : poDS->nRasterYSize = nLines;
2471 216 : poDS->m_osImageFilename = osImageFullFilename;
2472 216 : poDS->m_fpImage = fp;
2473 216 : poDS->m_bIsLSB = bLSBOrder;
2474 :
2475 35 : if (l_nBands > 1 && dimSemantics[0] == 'B' &&
2476 251 : dimSemantics[1] == 'L' && dimSemantics[2] == 'S')
2477 : {
2478 27 : poDS->GDALDataset::SetMetadataItem(GDALMD_INTERLEAVE, "BAND",
2479 : GDAL_MDD_IMAGE_STRUCTURE);
2480 : }
2481 35 : if (l_nBands > 1 && dimSemantics[0] == 'L' &&
2482 251 : dimSemantics[1] == 'S' && dimSemantics[2] == 'B')
2483 : {
2484 6 : poDS->GDALDataset::SetMetadataItem(GDALMD_INTERLEAVE, "PIXEL",
2485 : GDAL_MDD_IMAGE_STRUCTURE);
2486 : }
2487 :
2488 216 : CPLXMLNode *psOS = CPLGetXMLNode(psSubIter, "Object_Statistics");
2489 216 : const char *pszMin = nullptr;
2490 216 : const char *pszMax = nullptr;
2491 216 : const char *pszMean = nullptr;
2492 216 : const char *pszStdDev = nullptr;
2493 216 : if (psOS)
2494 : {
2495 17 : pszMin = CPLGetXMLValue(psOS, "minimum", nullptr);
2496 17 : pszMax = CPLGetXMLValue(psOS, "maximum", nullptr);
2497 17 : pszMean = CPLGetXMLValue(psOS, "mean", nullptr);
2498 17 : pszStdDev = CPLGetXMLValue(psOS, "standard_deviation", nullptr);
2499 : }
2500 :
2501 501 : for (int i = 0; i < l_nBands; i++)
2502 : {
2503 285 : vsi_l_offset nThisBandOffset = nOffset + nBandOffset * i;
2504 285 : if (bBottomToTop)
2505 : {
2506 2 : nThisBandOffset +=
2507 2 : static_cast<vsi_l_offset>(nLines - 1) * nLineOffset;
2508 : }
2509 285 : if (bRightToLeft)
2510 : {
2511 1 : nThisBandOffset +=
2512 1 : static_cast<vsi_l_offset>(nSamples - 1) * nPixelOffset;
2513 : }
2514 : auto poBand = std::make_unique<PDS4RawRasterBand>(
2515 285 : poDS.get(), i + 1, poDS->m_fpImage, nThisBandOffset,
2516 285 : bRightToLeft ? -nPixelOffset : nPixelOffset,
2517 285 : bBottomToTop ? -nLineOffset : nLineOffset, eDT,
2518 285 : bLSBOrder ? RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN
2519 285 : : RawRasterBand::ByteOrder::ORDER_BIG_ENDIAN);
2520 285 : if (!poBand->IsValid())
2521 0 : return nullptr;
2522 285 : if (bNoDataSet)
2523 : {
2524 75 : poBand->SetNoDataValue(dfNoData);
2525 : }
2526 210 : else if (bNoDataSetInt64)
2527 : {
2528 2 : poBand->SetNoDataValueAsInt64(nNoDataInt64);
2529 : }
2530 208 : else if (bNoDataSetUInt64)
2531 : {
2532 2 : poBand->SetNoDataValueAsUInt64(nNoDataUInt64);
2533 : }
2534 285 : poBand->SetOffset(dfValueOffset);
2535 285 : poBand->SetScale(dfValueScale);
2536 :
2537 285 : if (l_nBands == 1)
2538 : {
2539 181 : if (pszMin)
2540 : {
2541 17 : poBand->GDALRasterBand::SetMetadataItem(
2542 : "STATISTICS_MINIMUM", pszMin);
2543 : }
2544 181 : if (pszMax)
2545 : {
2546 17 : poBand->GDALRasterBand::SetMetadataItem(
2547 : "STATISTICS_MAXIMUM", pszMax);
2548 : }
2549 181 : if (pszMean)
2550 : {
2551 17 : poBand->GDALRasterBand::SetMetadataItem(
2552 : "STATISTICS_MEAN", pszMean);
2553 : }
2554 181 : if (pszStdDev)
2555 : {
2556 17 : poBand->GDALRasterBand::SetMetadataItem(
2557 : "STATISTICS_STDDEV", pszStdDev);
2558 : }
2559 : }
2560 :
2561 : // Only instantiate explicit mask band if we have at least one
2562 : // special constant (that is not the missing_constant,
2563 : // already exposed as nodata value)
2564 558 : if (!GDALDataTypeIsComplex(eDT) &&
2565 538 : (CPLTestBool(CPLGetConfigOption("PDS4_FORCE_MASK", "NO")) ||
2566 499 : adfConstants.size() >= 2 ||
2567 282 : (adfConstants.size() == 1 && !bNoDataSet)))
2568 : {
2569 86 : poBand->SetMaskBand(std::make_unique<PDS4MaskBand>(
2570 86 : poBand.get(), adfConstants));
2571 : }
2572 :
2573 285 : poDS->SetBand(i + 1, std::move(poBand));
2574 : }
2575 : }
2576 : }
2577 :
2578 300 : if (nFAOIdxLookup < 0 && aosSubdatasets.size() > 2)
2579 : {
2580 10 : poDS->GDALDataset::SetMetadata(aosSubdatasets.List(),
2581 : GDAL_MDD_SUBDATASETS);
2582 : }
2583 290 : else if (poDS->nBands == 0 &&
2584 300 : (poOpenInfo->nOpenFlags & GDAL_OF_RASTER) != 0 &&
2585 10 : (poOpenInfo->nOpenFlags & GDAL_OF_VECTOR) == 0)
2586 : {
2587 5 : return nullptr;
2588 : }
2589 285 : else if (poDS->m_apoLayers.empty() &&
2590 395 : (poOpenInfo->nOpenFlags & GDAL_OF_VECTOR) != 0 &&
2591 110 : (poOpenInfo->nOpenFlags & GDAL_OF_RASTER) == 0)
2592 : {
2593 0 : return nullptr;
2594 : }
2595 :
2596 : // Expose XML content in xml:PDS4 metadata domain
2597 295 : GByte *pabyRet = nullptr;
2598 295 : CPL_IGNORE_RET_VAL(VSIIngestFile(nullptr, osXMLFilename, &pabyRet, nullptr,
2599 : 10 * 1024 * 1024));
2600 295 : if (pabyRet)
2601 : {
2602 295 : char *apszXML[2] = {reinterpret_cast<char *>(pabyRet), nullptr};
2603 295 : poDS->GDALDataset::SetMetadata(apszXML, "xml:PDS4");
2604 : }
2605 295 : VSIFree(pabyRet);
2606 :
2607 : /*--------------------------------------------------------------------------*/
2608 : /* Parse georeferencing info */
2609 : /*--------------------------------------------------------------------------*/
2610 295 : poDS->ReadGeoreferencing(psProduct);
2611 :
2612 : /*--------------------------------------------------------------------------*/
2613 : /* Check for overviews */
2614 : /*--------------------------------------------------------------------------*/
2615 295 : poDS->oOvManager.Initialize(poDS.get(), poOpenInfo->pszFilename);
2616 :
2617 : /*--------------------------------------------------------------------------*/
2618 : /* Initialize any PAM information */
2619 : /*--------------------------------------------------------------------------*/
2620 295 : poDS->SetDescription(poOpenInfo->pszFilename);
2621 295 : poDS->TryLoadXML();
2622 :
2623 295 : return poDS;
2624 : }
2625 :
2626 : /************************************************************************/
2627 : /* IsCARTVersionGTE() */
2628 : /************************************************************************/
2629 :
2630 : // Returns true is pszCur >= pszRef
2631 : // Must be things like 1900, 1B00, 1D00_1933 ...
2632 489 : static bool IsCARTVersionGTE(const char *pszCur, const char *pszRef)
2633 : {
2634 489 : return strcmp(pszCur, pszRef) >= 0;
2635 : }
2636 :
2637 : /************************************************************************/
2638 : /* WriteGeoreferencing() */
2639 : /************************************************************************/
2640 :
2641 98 : void PDS4Dataset::WriteGeoreferencing(CPLXMLNode *psCart,
2642 : const char *pszCARTVersion)
2643 : {
2644 98 : bool bHasBoundingBox = false;
2645 98 : double adfX[4] = {0};
2646 98 : double adfY[4] = {0};
2647 98 : CPLString osPrefix;
2648 98 : const char *pszColon = strchr(psCart->pszValue, ':');
2649 98 : if (pszColon)
2650 98 : osPrefix.assign(psCart->pszValue, pszColon - psCart->pszValue + 1);
2651 :
2652 98 : if (m_bGotTransform)
2653 : {
2654 96 : bHasBoundingBox = true;
2655 :
2656 : // upper left
2657 96 : adfX[0] = m_gt.xorig;
2658 96 : adfY[0] = m_gt.yorig;
2659 :
2660 : // upper right
2661 96 : adfX[1] = m_gt.xorig + m_gt.xscale * nRasterXSize;
2662 96 : adfY[1] = m_gt.yorig;
2663 :
2664 : // lower left
2665 96 : adfX[2] = m_gt.xorig;
2666 96 : adfY[2] = m_gt.yorig + m_gt.yscale * nRasterYSize;
2667 :
2668 : // lower right
2669 96 : adfX[3] = m_gt.xorig + m_gt.xscale * nRasterXSize;
2670 96 : adfY[3] = m_gt.yorig + m_gt.yscale * nRasterYSize;
2671 : }
2672 : else
2673 : {
2674 2 : OGRLayer *poLayer = GetLayer(0);
2675 2 : OGREnvelope sEnvelope;
2676 2 : if (poLayer->GetExtent(&sEnvelope) == OGRERR_NONE)
2677 : {
2678 1 : bHasBoundingBox = true;
2679 :
2680 1 : adfX[0] = sEnvelope.MinX;
2681 1 : adfY[0] = sEnvelope.MaxY;
2682 :
2683 1 : adfX[1] = sEnvelope.MaxX;
2684 1 : adfY[1] = sEnvelope.MaxY;
2685 :
2686 1 : adfX[2] = sEnvelope.MinX;
2687 1 : adfY[2] = sEnvelope.MinY;
2688 :
2689 1 : adfX[3] = sEnvelope.MaxX;
2690 1 : adfY[3] = sEnvelope.MinY;
2691 : }
2692 : }
2693 :
2694 98 : if (bHasBoundingBox && !m_oSRS.IsGeographic())
2695 : {
2696 33 : bHasBoundingBox = false;
2697 33 : OGRSpatialReference *poSRSLongLat = m_oSRS.CloneGeogCS();
2698 33 : if (poSRSLongLat)
2699 : {
2700 33 : poSRSLongLat->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
2701 : OGRCoordinateTransformation *poCT =
2702 33 : OGRCreateCoordinateTransformation(&m_oSRS, poSRSLongLat);
2703 33 : if (poCT)
2704 : {
2705 33 : if (poCT->Transform(4, adfX, adfY))
2706 : {
2707 33 : bHasBoundingBox = true;
2708 : }
2709 33 : delete poCT;
2710 : }
2711 33 : delete poSRSLongLat;
2712 : }
2713 : }
2714 :
2715 98 : if (!bHasBoundingBox)
2716 : {
2717 : // Write dummy values
2718 1 : adfX[0] = -180.0;
2719 1 : adfY[0] = 90.0;
2720 1 : adfX[1] = 180.0;
2721 1 : adfY[1] = 90.0;
2722 1 : adfX[2] = -180.0;
2723 1 : adfY[2] = -90.0;
2724 1 : adfX[3] = 180.0;
2725 1 : adfY[3] = -90.0;
2726 : }
2727 :
2728 196 : const char *pszLongitudeDirection = CSLFetchNameValueDef(
2729 98 : m_papszCreationOptions, "LONGITUDE_DIRECTION", "Positive East");
2730 98 : const double dfLongitudeMultiplier =
2731 98 : EQUAL(pszLongitudeDirection, "Positive West") ? -1 : 1;
2732 212 : const auto FixLong = [dfLongitudeMultiplier](double dfLon)
2733 212 : { return dfLon * dfLongitudeMultiplier; };
2734 :
2735 : // Note: starting with CART 1900, Spatial_Domain is actually optional
2736 98 : CPLXMLNode *psSD = CPLCreateXMLNode(psCart, CXT_Element,
2737 196 : (osPrefix + "Spatial_Domain").c_str());
2738 98 : CPLXMLNode *psBC = CPLCreateXMLNode(
2739 196 : psSD, CXT_Element, (osPrefix + "Bounding_Coordinates").c_str());
2740 :
2741 : const char *pszBoundingDegrees =
2742 98 : CSLFetchNameValue(m_papszCreationOptions, "BOUNDING_DEGREES");
2743 98 : double dfWest = FixLong(
2744 98 : std::min(std::min(adfX[0], adfX[1]), std::min(adfX[2], adfX[3])));
2745 98 : double dfEast = FixLong(
2746 98 : std::max(std::max(adfX[0], adfX[1]), std::max(adfX[2], adfX[3])));
2747 : double dfNorth =
2748 98 : std::max(std::max(adfY[0], adfY[1]), std::max(adfY[2], adfY[3]));
2749 : double dfSouth =
2750 98 : std::min(std::min(adfY[0], adfY[1]), std::min(adfY[2], adfY[3]));
2751 98 : if (pszBoundingDegrees)
2752 : {
2753 1 : char **papszTokens = CSLTokenizeString2(pszBoundingDegrees, ",", 0);
2754 1 : if (CSLCount(papszTokens) == 4)
2755 : {
2756 1 : dfWest = CPLAtof(papszTokens[0]);
2757 1 : dfSouth = CPLAtof(papszTokens[1]);
2758 1 : dfEast = CPLAtof(papszTokens[2]);
2759 1 : dfNorth = CPLAtof(papszTokens[3]);
2760 : }
2761 1 : CSLDestroy(papszTokens);
2762 : }
2763 :
2764 196 : CPLAddXMLAttributeAndValue(
2765 : CPLCreateXMLElementAndValue(
2766 196 : psBC, (osPrefix + "west_bounding_coordinate").c_str(),
2767 : CPLSPrintf("%.17g", dfWest)),
2768 : "unit", "deg");
2769 196 : CPLAddXMLAttributeAndValue(
2770 : CPLCreateXMLElementAndValue(
2771 196 : psBC, (osPrefix + "east_bounding_coordinate").c_str(),
2772 : CPLSPrintf("%.17g", dfEast)),
2773 : "unit", "deg");
2774 196 : CPLAddXMLAttributeAndValue(
2775 : CPLCreateXMLElementAndValue(
2776 196 : psBC, (osPrefix + "north_bounding_coordinate").c_str(),
2777 : CPLSPrintf("%.17g", dfNorth)),
2778 : "unit", "deg");
2779 196 : CPLAddXMLAttributeAndValue(
2780 : CPLCreateXMLElementAndValue(
2781 196 : psBC, (osPrefix + "south_bounding_coordinate").c_str(),
2782 : CPLSPrintf("%.17g", dfSouth)),
2783 : "unit", "deg");
2784 :
2785 : CPLXMLNode *psSRI =
2786 98 : CPLCreateXMLNode(psCart, CXT_Element,
2787 196 : (osPrefix + "Spatial_Reference_Information").c_str());
2788 98 : CPLXMLNode *psHCSD = CPLCreateXMLNode(
2789 : psSRI, CXT_Element,
2790 196 : (osPrefix + "Horizontal_Coordinate_System_Definition").c_str());
2791 :
2792 98 : double dfUnrotatedULX = m_gt.xorig;
2793 98 : double dfUnrotatedULY = m_gt.yorig;
2794 98 : double dfUnrotatedResX = m_gt.xscale;
2795 98 : double dfUnrotatedResY = m_gt.yscale;
2796 98 : double dfMapProjectionRotation = 0.0;
2797 98 : if (m_gt.xscale == 0.0 && m_gt.xrot > 0.0 && m_gt.yrot > 0.0 &&
2798 1 : m_gt.yscale == 0.0)
2799 : {
2800 1 : dfUnrotatedULX = m_gt.yorig;
2801 1 : dfUnrotatedULY = -m_gt.xorig;
2802 1 : dfUnrotatedResX = m_gt.yrot;
2803 1 : dfUnrotatedResY = -m_gt.xrot;
2804 1 : dfMapProjectionRotation = 90.0;
2805 : }
2806 :
2807 98 : if (GetRasterCount() || m_oSRS.IsProjected())
2808 : {
2809 97 : CPLXMLNode *psPlanar = CPLCreateXMLNode(psHCSD, CXT_Element,
2810 194 : (osPrefix + "Planar").c_str());
2811 97 : CPLXMLNode *psMP = CPLCreateXMLNode(
2812 194 : psPlanar, CXT_Element, (osPrefix + "Map_Projection").c_str());
2813 97 : const char *pszProjection = m_oSRS.GetAttrValue("PROJECTION");
2814 194 : CPLString pszPDS4ProjectionName = "";
2815 : typedef std::pair<const char *, double> ProjParam;
2816 194 : std::vector<ProjParam> aoProjParams;
2817 :
2818 : const bool bUse_CART_1933_Or_Later =
2819 97 : IsCARTVersionGTE(pszCARTVersion, "1D00_1933");
2820 :
2821 : const bool bUse_CART_1950_Or_Later =
2822 97 : IsCARTVersionGTE(pszCARTVersion, "1G00_1950");
2823 :
2824 : const bool bUse_CART_1970_Or_Later =
2825 97 : IsCARTVersionGTE(pszCARTVersion, "1O00_1970");
2826 :
2827 97 : if (pszProjection == nullptr)
2828 : {
2829 63 : pszPDS4ProjectionName = "Equirectangular";
2830 63 : if (bUse_CART_1933_Or_Later)
2831 : {
2832 61 : aoProjParams.push_back(
2833 61 : ProjParam("latitude_of_projection_origin", 0.0));
2834 61 : aoProjParams.push_back(ProjParam("standard_parallel_1", 0.0));
2835 61 : aoProjParams.push_back(
2836 122 : ProjParam("longitude_of_central_meridian", 0.0));
2837 : }
2838 : else
2839 : {
2840 2 : aoProjParams.push_back(ProjParam("standard_parallel_1", 0.0));
2841 2 : aoProjParams.push_back(
2842 2 : ProjParam("longitude_of_central_meridian", 0.0));
2843 2 : aoProjParams.push_back(
2844 4 : ProjParam("latitude_of_projection_origin", 0.0));
2845 : }
2846 : }
2847 :
2848 34 : else if (EQUAL(pszProjection, SRS_PT_EQUIRECTANGULAR))
2849 : {
2850 2 : pszPDS4ProjectionName = "Equirectangular";
2851 2 : if (bUse_CART_1933_Or_Later)
2852 : {
2853 2 : aoProjParams.push_back(ProjParam(
2854 0 : "latitude_of_projection_origin",
2855 2 : m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2856 2 : aoProjParams.push_back(ProjParam(
2857 0 : "standard_parallel_1",
2858 2 : m_oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 1.0)));
2859 2 : aoProjParams.push_back(
2860 0 : ProjParam("longitude_of_central_meridian",
2861 4 : FixLong(m_oSRS.GetNormProjParm(
2862 : SRS_PP_CENTRAL_MERIDIAN, 0.0))));
2863 : }
2864 : else
2865 : {
2866 0 : aoProjParams.push_back(ProjParam(
2867 0 : "standard_parallel_1",
2868 0 : m_oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 1.0)));
2869 0 : aoProjParams.push_back(
2870 0 : ProjParam("longitude_of_central_meridian",
2871 0 : FixLong(m_oSRS.GetNormProjParm(
2872 : SRS_PP_CENTRAL_MERIDIAN, 0.0))));
2873 0 : aoProjParams.push_back(ProjParam(
2874 0 : "latitude_of_projection_origin",
2875 0 : m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2876 : }
2877 : }
2878 :
2879 32 : else if (EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP))
2880 : {
2881 1 : pszPDS4ProjectionName = "Lambert Conformal Conic";
2882 1 : if (bUse_CART_1933_Or_Later)
2883 : {
2884 1 : if (bUse_CART_1970_Or_Later)
2885 : {
2886 : // Note: in EPSG (and PROJ "metadata" part), there is
2887 : // no standard_parallel_1 parameter for LCC_1SP. But
2888 : // for historical reason we do map the latitude of origin
2889 : // to +lat_1 (in addition to +lat_0). So do the same here.
2890 1 : aoProjParams.push_back(
2891 0 : ProjParam("standard_parallel_1",
2892 2 : m_oSRS.GetNormProjParm(
2893 : SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2894 : }
2895 1 : aoProjParams.push_back(
2896 0 : ProjParam("longitude_of_central_meridian",
2897 1 : FixLong(m_oSRS.GetNormProjParm(
2898 : SRS_PP_CENTRAL_MERIDIAN, 0.0))));
2899 1 : aoProjParams.push_back(ProjParam(
2900 0 : "latitude_of_projection_origin",
2901 1 : m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2902 1 : aoProjParams.push_back(ProjParam(
2903 0 : "scale_factor_at_projection_origin",
2904 2 : m_oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0)));
2905 : }
2906 : else
2907 : {
2908 0 : aoProjParams.push_back(ProjParam(
2909 0 : "scale_factor_at_projection_origin",
2910 0 : m_oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0)));
2911 0 : aoProjParams.push_back(
2912 0 : ProjParam("longitude_of_central_meridian",
2913 0 : FixLong(m_oSRS.GetNormProjParm(
2914 : SRS_PP_CENTRAL_MERIDIAN, 0.0))));
2915 0 : aoProjParams.push_back(ProjParam(
2916 0 : "latitude_of_projection_origin",
2917 0 : m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2918 : }
2919 : }
2920 :
2921 31 : else if (EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP))
2922 : {
2923 1 : pszPDS4ProjectionName = "Lambert Conformal Conic";
2924 1 : aoProjParams.push_back(ProjParam(
2925 0 : "standard_parallel_1",
2926 1 : m_oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0)));
2927 1 : aoProjParams.push_back(ProjParam(
2928 0 : "standard_parallel_2",
2929 1 : m_oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_2, 0.0)));
2930 1 : aoProjParams.push_back(ProjParam(
2931 0 : "longitude_of_central_meridian",
2932 1 : FixLong(m_oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0))));
2933 1 : aoProjParams.push_back(ProjParam(
2934 0 : "latitude_of_projection_origin",
2935 2 : m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2936 : }
2937 :
2938 30 : else if (EQUAL(pszProjection,
2939 : SRS_PT_HOTINE_OBLIQUE_MERCATOR_AZIMUTH_CENTER))
2940 : {
2941 1 : pszPDS4ProjectionName = "Oblique Mercator";
2942 : // Proj params defined later
2943 : }
2944 :
2945 29 : else if (EQUAL(pszProjection,
2946 : SRS_PT_HOTINE_OBLIQUE_MERCATOR_TWO_POINT_NATURAL_ORIGIN))
2947 : {
2948 1 : pszPDS4ProjectionName = "Oblique Mercator";
2949 : // Proj params defined later
2950 : }
2951 :
2952 28 : else if (EQUAL(pszProjection, SRS_PT_POLAR_STEREOGRAPHIC))
2953 : {
2954 1 : pszPDS4ProjectionName = "Polar Stereographic";
2955 1 : if (bUse_CART_1950_Or_Later)
2956 : {
2957 1 : aoProjParams.push_back(
2958 0 : ProjParam("longitude_of_central_meridian",
2959 1 : FixLong(m_oSRS.GetNormProjParm(
2960 : SRS_PP_CENTRAL_MERIDIAN, 0.0))));
2961 1 : aoProjParams.push_back(ProjParam(
2962 0 : "latitude_of_projection_origin",
2963 1 : m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2964 1 : aoProjParams.push_back(ProjParam(
2965 0 : "scale_factor_at_projection_origin",
2966 2 : m_oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0)));
2967 : }
2968 : else
2969 : {
2970 0 : aoProjParams.push_back(
2971 0 : ProjParam(bUse_CART_1933_Or_Later
2972 0 : ? "longitude_of_central_meridian"
2973 : : "straight_vertical_longitude_from_pole",
2974 0 : FixLong(m_oSRS.GetNormProjParm(
2975 : SRS_PP_CENTRAL_MERIDIAN, 0.0))));
2976 0 : aoProjParams.push_back(ProjParam(
2977 0 : "scale_factor_at_projection_origin",
2978 0 : m_oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0)));
2979 0 : aoProjParams.push_back(ProjParam(
2980 0 : "latitude_of_projection_origin",
2981 0 : m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2982 : }
2983 : }
2984 :
2985 27 : else if (EQUAL(pszProjection, SRS_PT_POLYCONIC))
2986 : {
2987 1 : pszPDS4ProjectionName = "Polyconic";
2988 1 : aoProjParams.push_back(ProjParam(
2989 0 : "longitude_of_central_meridian",
2990 1 : m_oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2991 1 : aoProjParams.push_back(ProjParam(
2992 0 : "latitude_of_projection_origin",
2993 2 : m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2994 : }
2995 26 : else if (EQUAL(pszProjection, SRS_PT_SINUSOIDAL))
2996 : {
2997 2 : pszPDS4ProjectionName = "Sinusoidal";
2998 2 : aoProjParams.push_back(ProjParam(
2999 0 : "longitude_of_central_meridian",
3000 2 : FixLong(m_oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0))));
3001 2 : aoProjParams.push_back(ProjParam(
3002 0 : "latitude_of_projection_origin",
3003 4 : m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
3004 : }
3005 :
3006 24 : else if (EQUAL(pszProjection, SRS_PT_TRANSVERSE_MERCATOR))
3007 : {
3008 18 : pszPDS4ProjectionName = "Transverse Mercator";
3009 18 : aoProjParams.push_back(
3010 0 : ProjParam("scale_factor_at_central_meridian",
3011 18 : m_oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0)));
3012 18 : aoProjParams.push_back(ProjParam(
3013 0 : "longitude_of_central_meridian",
3014 18 : m_oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0)));
3015 18 : aoProjParams.push_back(ProjParam(
3016 0 : "latitude_of_projection_origin",
3017 36 : m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
3018 : }
3019 6 : else if (EQUAL(pszProjection, SRS_PT_ORTHOGRAPHIC))
3020 : {
3021 1 : pszPDS4ProjectionName = "Orthographic";
3022 1 : aoProjParams.push_back(ProjParam(
3023 0 : "longitude_of_central_meridian",
3024 1 : FixLong(m_oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0))));
3025 1 : aoProjParams.push_back(ProjParam(
3026 0 : "latitude_of_projection_origin",
3027 2 : m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
3028 : }
3029 :
3030 5 : else if (EQUAL(pszProjection, SRS_PT_MERCATOR_1SP))
3031 : {
3032 2 : pszPDS4ProjectionName = "Mercator";
3033 2 : aoProjParams.push_back(ProjParam(
3034 0 : "longitude_of_central_meridian",
3035 2 : FixLong(m_oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0))));
3036 2 : aoProjParams.push_back(ProjParam(
3037 0 : "latitude_of_projection_origin",
3038 2 : m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
3039 2 : aoProjParams.push_back(
3040 0 : ProjParam("scale_factor_at_projection_origin",
3041 4 : m_oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0)));
3042 : }
3043 :
3044 3 : else if (EQUAL(pszProjection, SRS_PT_MERCATOR_2SP))
3045 : {
3046 1 : pszPDS4ProjectionName = "Mercator";
3047 1 : aoProjParams.push_back(ProjParam(
3048 0 : "standard_parallel_1",
3049 1 : m_oSRS.GetNormProjParm(SRS_PP_STANDARD_PARALLEL_1, 0.0)));
3050 1 : aoProjParams.push_back(ProjParam(
3051 0 : "longitude_of_central_meridian",
3052 1 : FixLong(m_oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0))));
3053 1 : aoProjParams.push_back(ProjParam(
3054 0 : "latitude_of_projection_origin",
3055 2 : m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
3056 : }
3057 :
3058 2 : else if (EQUAL(pszProjection, SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA))
3059 : {
3060 1 : pszPDS4ProjectionName = "Lambert Azimuthal Equal Area";
3061 1 : aoProjParams.push_back(ProjParam(
3062 0 : "longitude_of_central_meridian",
3063 1 : FixLong(m_oSRS.GetNormProjParm(SRS_PP_CENTRAL_MERIDIAN, 0.0))));
3064 1 : aoProjParams.push_back(ProjParam(
3065 0 : "latitude_of_projection_origin",
3066 2 : m_oSRS.GetNormProjParm(SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
3067 : }
3068 :
3069 1 : else if (EQUAL(pszProjection, "custom_proj4"))
3070 : {
3071 : const char *pszProj4 =
3072 1 : m_oSRS.GetExtension("PROJCS", "PROJ4", nullptr);
3073 1 : if (pszProj4 && strstr(pszProj4, "+proj=ob_tran") &&
3074 1 : strstr(pszProj4, "+o_proj=eqc"))
3075 : {
3076 1 : pszPDS4ProjectionName = "Oblique Cylindrical";
3077 : const auto FetchParam =
3078 3 : [](const char *pszProj4Str, const char *pszKey)
3079 : {
3080 6 : CPLString needle;
3081 3 : needle.Printf("+%s=", pszKey);
3082 3 : const char *pszVal = strstr(pszProj4Str, needle.c_str());
3083 3 : if (pszVal)
3084 3 : return CPLAtof(pszVal + needle.size());
3085 0 : return 0.0;
3086 : };
3087 :
3088 1 : double dfLonP = FetchParam(pszProj4, "o_lon_p");
3089 1 : double dfLatP = FetchParam(pszProj4, "o_lat_p");
3090 1 : double dfLon0 = FetchParam(pszProj4, "lon_0");
3091 1 : double dfPoleRotation = -dfLonP;
3092 1 : double dfPoleLatitude = 180 - dfLatP;
3093 1 : double dfPoleLongitude = dfLon0;
3094 :
3095 1 : aoProjParams.push_back(ProjParam("map_projection_rotation",
3096 : dfMapProjectionRotation));
3097 1 : aoProjParams.push_back(
3098 1 : ProjParam("oblique_proj_pole_latitude", dfPoleLatitude));
3099 1 : aoProjParams.push_back(ProjParam("oblique_proj_pole_longitude",
3100 1 : FixLong(dfPoleLongitude)));
3101 1 : aoProjParams.push_back(
3102 2 : ProjParam("oblique_proj_pole_rotation", dfPoleRotation));
3103 : }
3104 : else
3105 : {
3106 0 : CPLError(CE_Warning, CPLE_NotSupported,
3107 : "Projection %s not supported", pszProjection);
3108 : }
3109 : }
3110 : else
3111 : {
3112 0 : CPLError(CE_Warning, CPLE_NotSupported,
3113 : "Projection %s not supported", pszProjection);
3114 : }
3115 194 : CPLCreateXMLElementAndValue(psMP,
3116 194 : (osPrefix + "map_projection_name").c_str(),
3117 : pszPDS4ProjectionName);
3118 97 : CPLXMLNode *psProj = CPLCreateXMLNode(
3119 : psMP, CXT_Element,
3120 194 : CPLString(osPrefix + pszPDS4ProjectionName).replaceAll(' ', '_'));
3121 380 : for (size_t i = 0; i < aoProjParams.size(); i++)
3122 : {
3123 566 : CPLXMLNode *psParam = CPLCreateXMLElementAndValue(
3124 566 : psProj, (osPrefix + aoProjParams[i].first).c_str(),
3125 283 : CPLSPrintf("%.17g", aoProjParams[i].second));
3126 283 : if (!STARTS_WITH(aoProjParams[i].first, "scale_factor"))
3127 : {
3128 261 : CPLAddXMLAttributeAndValue(psParam, "unit", "deg");
3129 : }
3130 : }
3131 :
3132 97 : if (pszProjection &&
3133 34 : EQUAL(pszProjection, SRS_PT_HOTINE_OBLIQUE_MERCATOR_AZIMUTH_CENTER))
3134 : {
3135 : CPLXMLNode *psOLA =
3136 1 : CPLCreateXMLNode(nullptr, CXT_Element,
3137 2 : (osPrefix + "Oblique_Line_Azimuth").c_str());
3138 2 : CPLAddXMLAttributeAndValue(
3139 : CPLCreateXMLElementAndValue(
3140 2 : psOLA, (osPrefix + "azimuthal_angle").c_str(),
3141 : CPLSPrintf("%.17g",
3142 : m_oSRS.GetNormProjParm(SRS_PP_AZIMUTH, 0.0))),
3143 : "unit", "deg");
3144 : ;
3145 : // Not completely sure of this
3146 2 : CPLAddXMLAttributeAndValue(
3147 : CPLCreateXMLElementAndValue(
3148 : psOLA,
3149 2 : (osPrefix + "azimuth_measure_point_longitude").c_str(),
3150 : CPLSPrintf("%.17g", FixLong(m_oSRS.GetNormProjParm(
3151 : SRS_PP_CENTRAL_MERIDIAN, 0.0)))),
3152 : "unit", "deg");
3153 :
3154 1 : if (bUse_CART_1933_Or_Later)
3155 : {
3156 1 : CPLAddXMLChild(psProj, psOLA);
3157 :
3158 1 : CPLAddXMLAttributeAndValue(
3159 : CPLCreateXMLElementAndValue(
3160 : psProj,
3161 2 : (osPrefix + "longitude_of_central_meridian").c_str(),
3162 : "0"),
3163 : "unit", "deg");
3164 :
3165 : const double dfScaleFactor =
3166 1 : m_oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 0.0);
3167 1 : if (dfScaleFactor != 1.0)
3168 : {
3169 0 : CPLError(CE_Warning, CPLE_NotSupported,
3170 : "Scale factor on initial support = %.17g cannot "
3171 : "be encoded in PDS4",
3172 : dfScaleFactor);
3173 : }
3174 : }
3175 : else
3176 : {
3177 0 : CPLCreateXMLElementAndValue(
3178 : psProj,
3179 0 : (osPrefix + "scale_factor_at_projection_origin").c_str(),
3180 : CPLSPrintf("%.17g", m_oSRS.GetNormProjParm(
3181 : SRS_PP_SCALE_FACTOR, 0.0)));
3182 :
3183 0 : CPLAddXMLChild(psProj, psOLA);
3184 : }
3185 :
3186 2 : CPLAddXMLAttributeAndValue(
3187 : CPLCreateXMLElementAndValue(
3188 : psProj,
3189 2 : (osPrefix + "latitude_of_projection_origin").c_str(),
3190 : CPLSPrintf("%.17g", m_oSRS.GetNormProjParm(
3191 : SRS_PP_LATITUDE_OF_ORIGIN, 0.0))),
3192 1 : "unit", "deg");
3193 : }
3194 96 : else if (pszProjection &&
3195 33 : EQUAL(pszProjection,
3196 : SRS_PT_HOTINE_OBLIQUE_MERCATOR_TWO_POINT_NATURAL_ORIGIN))
3197 : {
3198 1 : if (bUse_CART_1933_Or_Later)
3199 : {
3200 : const double dfScaleFactor =
3201 1 : m_oSRS.GetNormProjParm(SRS_PP_SCALE_FACTOR, 0.0);
3202 1 : if (dfScaleFactor != 1.0)
3203 : {
3204 0 : CPLError(CE_Warning, CPLE_NotSupported,
3205 : "Scale factor on initial support = %.17g cannot "
3206 : "be encoded in PDS4",
3207 : dfScaleFactor);
3208 : }
3209 : }
3210 : else
3211 : {
3212 0 : CPLCreateXMLElementAndValue(
3213 : psProj,
3214 0 : (osPrefix + "scale_factor_at_projection_origin").c_str(),
3215 : CPLSPrintf("%.17g", m_oSRS.GetNormProjParm(
3216 : SRS_PP_SCALE_FACTOR, 0.0)));
3217 : }
3218 :
3219 1 : CPLXMLNode *psOLP = CPLCreateXMLNode(
3220 2 : psProj, CXT_Element, (osPrefix + "Oblique_Line_Point").c_str());
3221 1 : CPLXMLNode *psOLPG1 = CPLCreateXMLNode(
3222 : psOLP, CXT_Element,
3223 2 : (osPrefix + "Oblique_Line_Point_Group").c_str());
3224 2 : CPLAddXMLAttributeAndValue(
3225 : CPLCreateXMLElementAndValue(
3226 2 : psOLPG1, (osPrefix + "oblique_line_latitude").c_str(),
3227 : CPLSPrintf("%.17g", m_oSRS.GetNormProjParm(
3228 : SRS_PP_LATITUDE_OF_POINT_1, 0.0))),
3229 : "unit", "deg");
3230 2 : CPLAddXMLAttributeAndValue(
3231 : CPLCreateXMLElementAndValue(
3232 2 : psOLPG1, (osPrefix + "oblique_line_longitude").c_str(),
3233 : CPLSPrintf("%.17g",
3234 : FixLong(m_oSRS.GetNormProjParm(
3235 : SRS_PP_LONGITUDE_OF_POINT_1, 0.0)))),
3236 : "unit", "deg");
3237 1 : CPLXMLNode *psOLPG2 = CPLCreateXMLNode(
3238 : psOLP, CXT_Element,
3239 2 : (osPrefix + "Oblique_Line_Point_Group").c_str());
3240 2 : CPLAddXMLAttributeAndValue(
3241 : CPLCreateXMLElementAndValue(
3242 2 : psOLPG2, (osPrefix + "oblique_line_latitude").c_str(),
3243 : CPLSPrintf("%.17g", m_oSRS.GetNormProjParm(
3244 : SRS_PP_LATITUDE_OF_POINT_2, 0.0))),
3245 : "unit", "deg");
3246 2 : CPLAddXMLAttributeAndValue(
3247 : CPLCreateXMLElementAndValue(
3248 2 : psOLPG2, (osPrefix + "oblique_line_longitude").c_str(),
3249 : CPLSPrintf("%.17g", m_oSRS.GetNormProjParm(
3250 : SRS_PP_LONGITUDE_OF_POINT_2, 0.0))),
3251 : "unit", "deg");
3252 :
3253 1 : if (bUse_CART_1933_Or_Later)
3254 : {
3255 1 : CPLAddXMLAttributeAndValue(
3256 : CPLCreateXMLElementAndValue(
3257 : psProj,
3258 2 : (osPrefix + "longitude_of_central_meridian").c_str(),
3259 : "0"),
3260 : "unit", "deg");
3261 : }
3262 :
3263 2 : CPLAddXMLAttributeAndValue(
3264 : CPLCreateXMLElementAndValue(
3265 : psProj,
3266 2 : (osPrefix + "latitude_of_projection_origin").c_str(),
3267 : CPLSPrintf("%.17g", FixLong(m_oSRS.GetNormProjParm(
3268 : SRS_PP_LATITUDE_OF_ORIGIN, 0.0)))),
3269 : "unit", "deg");
3270 : }
3271 :
3272 97 : CPLXMLNode *psCR = nullptr;
3273 97 : if (m_bGotTransform || !IsCARTVersionGTE(pszCARTVersion, "1B00"))
3274 : {
3275 96 : CPLXMLNode *psPCI = CPLCreateXMLNode(
3276 : psPlanar, CXT_Element,
3277 192 : (osPrefix + "Planar_Coordinate_Information").c_str());
3278 96 : CPLCreateXMLElementAndValue(
3279 192 : psPCI, (osPrefix + "planar_coordinate_encoding_method").c_str(),
3280 : "Coordinate Pair");
3281 96 : psCR = CPLCreateXMLNode(
3282 : psPCI, CXT_Element,
3283 192 : (osPrefix + "Coordinate_Representation").c_str());
3284 : }
3285 97 : const double dfLinearUnits = m_oSRS.GetLinearUnits();
3286 97 : const double dfDegToMeter = m_oSRS.GetSemiMajor() * M_PI / 180.0;
3287 :
3288 97 : if (psCR == nullptr)
3289 : {
3290 : // do nothing
3291 : }
3292 96 : else if (!m_bGotTransform)
3293 : {
3294 0 : CPLAddXMLAttributeAndValue(
3295 : CPLCreateXMLElementAndValue(
3296 0 : psCR, (osPrefix + "pixel_resolution_x").c_str(), "0"),
3297 : "unit", "m/pixel");
3298 0 : CPLAddXMLAttributeAndValue(
3299 : CPLCreateXMLElementAndValue(
3300 0 : psCR, (osPrefix + "pixel_resolution_y").c_str(), "0"),
3301 : "unit", "m/pixel");
3302 0 : CPLAddXMLAttributeAndValue(
3303 : CPLCreateXMLElementAndValue(
3304 0 : psCR, (osPrefix + "pixel_scale_x").c_str(), "0"),
3305 : "unit", "pixel/deg");
3306 0 : CPLAddXMLAttributeAndValue(
3307 : CPLCreateXMLElementAndValue(
3308 0 : psCR, (osPrefix + "pixel_scale_y").c_str(), "0"),
3309 : "unit", "pixel/deg");
3310 : }
3311 96 : else if (m_oSRS.IsGeographic())
3312 : {
3313 126 : CPLAddXMLAttributeAndValue(
3314 : CPLCreateXMLElementAndValue(
3315 126 : psCR, (osPrefix + "pixel_resolution_x").c_str(),
3316 : CPLSPrintf("%.17g", dfUnrotatedResX * dfDegToMeter)),
3317 : "unit", "m/pixel");
3318 63 : CPLAddXMLAttributeAndValue(
3319 : CPLCreateXMLElementAndValue(
3320 126 : psCR, (osPrefix + "pixel_resolution_y").c_str(),
3321 63 : CPLSPrintf("%.17g", -dfUnrotatedResY * dfDegToMeter)),
3322 : "unit", "m/pixel");
3323 126 : CPLAddXMLAttributeAndValue(
3324 : CPLCreateXMLElementAndValue(
3325 126 : psCR, (osPrefix + "pixel_scale_x").c_str(),
3326 : CPLSPrintf("%.17g", 1.0 / (dfUnrotatedResX))),
3327 : "unit", "pixel/deg");
3328 126 : CPLAddXMLAttributeAndValue(
3329 : CPLCreateXMLElementAndValue(
3330 126 : psCR, (osPrefix + "pixel_scale_y").c_str(),
3331 : CPLSPrintf("%.17g", 1.0 / (-dfUnrotatedResY))),
3332 : "unit", "pixel/deg");
3333 : }
3334 33 : else if (m_oSRS.IsProjected())
3335 : {
3336 66 : CPLAddXMLAttributeAndValue(
3337 : CPLCreateXMLElementAndValue(
3338 66 : psCR, (osPrefix + "pixel_resolution_x").c_str(),
3339 : CPLSPrintf("%.17g", dfUnrotatedResX * dfLinearUnits)),
3340 : "unit", "m/pixel");
3341 33 : CPLAddXMLAttributeAndValue(
3342 : CPLCreateXMLElementAndValue(
3343 66 : psCR, (osPrefix + "pixel_resolution_y").c_str(),
3344 33 : CPLSPrintf("%.17g", -dfUnrotatedResY * dfLinearUnits)),
3345 : "unit", "m/pixel");
3346 33 : CPLAddXMLAttributeAndValue(
3347 : CPLCreateXMLElementAndValue(
3348 66 : psCR, (osPrefix + "pixel_scale_x").c_str(),
3349 : CPLSPrintf("%.17g", dfDegToMeter /
3350 33 : (dfUnrotatedResX * dfLinearUnits))),
3351 : "unit", "pixel/deg");
3352 33 : CPLAddXMLAttributeAndValue(
3353 : CPLCreateXMLElementAndValue(
3354 66 : psCR, (osPrefix + "pixel_scale_y").c_str(),
3355 33 : CPLSPrintf("%.17g", dfDegToMeter / (-dfUnrotatedResY *
3356 : dfLinearUnits))),
3357 : "unit", "pixel/deg");
3358 : }
3359 :
3360 97 : if (m_bGotTransform)
3361 : {
3362 : CPLXMLNode *psGT =
3363 96 : CPLCreateXMLNode(psPlanar, CXT_Element,
3364 192 : (osPrefix + "Geo_Transformation").c_str());
3365 : const double dfFalseEasting =
3366 96 : m_oSRS.GetNormProjParm(SRS_PP_FALSE_EASTING, 0.0);
3367 : const double dfFalseNorthing =
3368 96 : m_oSRS.GetNormProjParm(SRS_PP_FALSE_NORTHING, 0.0);
3369 96 : const double dfULX = -dfFalseEasting + dfUnrotatedULX;
3370 96 : const double dfULY = -dfFalseNorthing + dfUnrotatedULY;
3371 96 : if (m_oSRS.IsGeographic())
3372 : {
3373 126 : CPLAddXMLAttributeAndValue(
3374 : CPLCreateXMLElementAndValue(
3375 126 : psGT, (osPrefix + "upperleft_corner_x").c_str(),
3376 : CPLSPrintf("%.17g", dfULX * dfDegToMeter)),
3377 : "unit", "m");
3378 126 : CPLAddXMLAttributeAndValue(
3379 : CPLCreateXMLElementAndValue(
3380 126 : psGT, (osPrefix + "upperleft_corner_y").c_str(),
3381 : CPLSPrintf("%.17g", dfULY * dfDegToMeter)),
3382 : "unit", "m");
3383 : }
3384 33 : else if (m_oSRS.IsProjected())
3385 : {
3386 66 : CPLAddXMLAttributeAndValue(
3387 : CPLCreateXMLElementAndValue(
3388 66 : psGT, (osPrefix + "upperleft_corner_x").c_str(),
3389 : CPLSPrintf("%.17g", dfULX * dfLinearUnits)),
3390 : "unit", "m");
3391 66 : CPLAddXMLAttributeAndValue(
3392 : CPLCreateXMLElementAndValue(
3393 66 : psGT, (osPrefix + "upperleft_corner_y").c_str(),
3394 : CPLSPrintf("%.17g", dfULY * dfLinearUnits)),
3395 : "unit", "m");
3396 : }
3397 : }
3398 : }
3399 : else
3400 : {
3401 1 : CPLXMLNode *psGeographic = CPLCreateXMLNode(
3402 2 : psHCSD, CXT_Element, (osPrefix + "Geographic").c_str());
3403 1 : if (!IsCARTVersionGTE(pszCARTVersion, "1B00"))
3404 : {
3405 0 : CPLAddXMLAttributeAndValue(
3406 : CPLCreateXMLElementAndValue(
3407 0 : psGeographic, (osPrefix + "latitude_resolution").c_str(),
3408 : "0"),
3409 : "unit", "deg");
3410 0 : CPLAddXMLAttributeAndValue(
3411 : CPLCreateXMLElementAndValue(
3412 0 : psGeographic, (osPrefix + "longitude_resolution").c_str(),
3413 : "0"),
3414 : "unit", "deg");
3415 : }
3416 : }
3417 :
3418 98 : CPLXMLNode *psGM = CPLCreateXMLNode(psHCSD, CXT_Element,
3419 196 : (osPrefix + "Geodetic_Model").c_str());
3420 196 : const char *pszLatitudeType = CSLFetchNameValueDef(
3421 98 : m_papszCreationOptions, "LATITUDE_TYPE", "Planetocentric");
3422 : // Fix case
3423 98 : if (EQUAL(pszLatitudeType, "Planetocentric"))
3424 97 : pszLatitudeType = "Planetocentric";
3425 1 : else if (EQUAL(pszLatitudeType, "Planetographic"))
3426 1 : pszLatitudeType = "Planetographic";
3427 98 : CPLCreateXMLElementAndValue(psGM, (osPrefix + "latitude_type").c_str(),
3428 : pszLatitudeType);
3429 :
3430 98 : const char *pszDatum = m_oSRS.GetAttrValue("DATUM");
3431 98 : if (pszDatum && STARTS_WITH(pszDatum, "D_"))
3432 : {
3433 4 : CPLCreateXMLElementAndValue(psGM, (osPrefix + "spheroid_name").c_str(),
3434 : pszDatum + 2);
3435 : }
3436 94 : else if (pszDatum)
3437 : {
3438 94 : CPLCreateXMLElementAndValue(psGM, (osPrefix + "spheroid_name").c_str(),
3439 : pszDatum);
3440 : }
3441 :
3442 98 : double dfSemiMajor = m_oSRS.GetSemiMajor();
3443 98 : double dfSemiMinor = m_oSRS.GetSemiMinor();
3444 98 : const char *pszRadii = CSLFetchNameValue(m_papszCreationOptions, "RADII");
3445 98 : if (pszRadii)
3446 : {
3447 1 : char **papszTokens = CSLTokenizeString2(pszRadii, " ,", 0);
3448 1 : if (CSLCount(papszTokens) == 2)
3449 : {
3450 1 : dfSemiMajor = CPLAtof(papszTokens[0]);
3451 1 : dfSemiMinor = CPLAtof(papszTokens[1]);
3452 : }
3453 1 : CSLDestroy(papszTokens);
3454 : }
3455 :
3456 : const bool bUseLDD1930RadiusNames =
3457 98 : IsCARTVersionGTE(pszCARTVersion, "1B10_1930");
3458 :
3459 196 : CPLAddXMLAttributeAndValue(
3460 : CPLCreateXMLElementAndValue(
3461 : psGM,
3462 196 : (osPrefix +
3463 : (bUseLDD1930RadiusNames ? "a_axis_radius" : "semi_major_radius"))
3464 : .c_str(),
3465 : CPLSPrintf("%.17g", dfSemiMajor)),
3466 : "unit", "m");
3467 : // No, this is not a bug. The PDS4 b_axis_radius/semi_minor_radius is the
3468 : // minor radius on the equatorial plane. Which in WKT doesn't really exist,
3469 : // so reuse the WKT semi major
3470 196 : CPLAddXMLAttributeAndValue(
3471 : CPLCreateXMLElementAndValue(
3472 : psGM,
3473 196 : (osPrefix +
3474 : (bUseLDD1930RadiusNames ? "b_axis_radius" : "semi_minor_radius"))
3475 : .c_str(),
3476 : CPLSPrintf("%.17g", dfSemiMajor)),
3477 : "unit", "m");
3478 196 : CPLAddXMLAttributeAndValue(
3479 : CPLCreateXMLElementAndValue(
3480 : psGM,
3481 196 : (osPrefix +
3482 : (bUseLDD1930RadiusNames ? "c_axis_radius" : "polar_radius"))
3483 : .c_str(),
3484 : CPLSPrintf("%.17g", dfSemiMinor)),
3485 : "unit", "m");
3486 :
3487 : // Fix case
3488 98 : if (EQUAL(pszLongitudeDirection, "Positive East"))
3489 97 : pszLongitudeDirection = "Positive East";
3490 1 : else if (EQUAL(pszLongitudeDirection, "Positive West"))
3491 1 : pszLongitudeDirection = "Positive West";
3492 98 : CPLCreateXMLElementAndValue(psGM,
3493 196 : (osPrefix + "longitude_direction").c_str(),
3494 : pszLongitudeDirection);
3495 98 : }
3496 :
3497 : /************************************************************************/
3498 : /* SubstituteVariables() */
3499 : /************************************************************************/
3500 :
3501 16490 : void PDS4Dataset::SubstituteVariables(CPLXMLNode *psNode, char **papszDict)
3502 : {
3503 16490 : if (psNode->eType == CXT_Text && psNode->pszValue &&
3504 6668 : strstr(psNode->pszValue, "${"))
3505 : {
3506 1382 : CPLString osVal(psNode->pszValue);
3507 :
3508 1556 : if (strstr(psNode->pszValue, "${TITLE}") != nullptr &&
3509 174 : CSLFetchNameValue(papszDict, "VAR_TITLE") == nullptr)
3510 : {
3511 160 : const CPLString osTitle(CPLGetFilename(GetDescription()));
3512 160 : CPLError(CE_Warning, CPLE_AppDefined,
3513 : "VAR_TITLE not defined. Using %s by default",
3514 : osTitle.c_str());
3515 160 : osVal.replaceAll("${TITLE}", osTitle);
3516 : }
3517 :
3518 4599 : for (char **papszIter = papszDict; papszIter && *papszIter; papszIter++)
3519 : {
3520 3217 : if (STARTS_WITH_CI(*papszIter, "VAR_"))
3521 : {
3522 2746 : char *pszKey = nullptr;
3523 2746 : const char *pszValue = CPLParseNameValue(*papszIter, &pszKey);
3524 2746 : if (pszKey && pszValue)
3525 : {
3526 2746 : const char *pszVarName = pszKey + strlen("VAR_");
3527 : osVal.replaceAll(
3528 2746 : (CPLString("${") + pszVarName + "}").c_str(), pszValue);
3529 : osVal.replaceAll(
3530 5492 : CPLString(CPLString("${") + pszVarName + "}")
3531 2746 : .tolower()
3532 : .c_str(),
3533 8238 : CPLString(pszValue).tolower());
3534 2746 : CPLFree(pszKey);
3535 : }
3536 : }
3537 : }
3538 1382 : if (osVal.find("${") != std::string::npos)
3539 : {
3540 778 : CPLError(CE_Warning, CPLE_AppDefined, "%s could not be substituted",
3541 : osVal.c_str());
3542 : }
3543 1382 : CPLFree(psNode->pszValue);
3544 1382 : psNode->pszValue = CPLStrdup(osVal);
3545 : }
3546 :
3547 32799 : for (CPLXMLNode *psIter = psNode->psChild; psIter; psIter = psIter->psNext)
3548 : {
3549 16309 : SubstituteVariables(psIter, papszDict);
3550 : }
3551 16490 : }
3552 :
3553 : /************************************************************************/
3554 : /* InitImageFile() */
3555 : /************************************************************************/
3556 :
3557 93 : bool PDS4Dataset::InitImageFile()
3558 : {
3559 93 : m_bMustInitImageFile = false;
3560 :
3561 93 : int bHasNoData = FALSE;
3562 93 : double dfNoData = 0;
3563 93 : int bHasNoDataAsInt64 = FALSE;
3564 93 : int64_t nNoDataInt64 = 0;
3565 93 : int bHasNoDataAsUInt64 = FALSE;
3566 93 : uint64_t nNoDataUInt64 = 0;
3567 93 : const GDALDataType eDT = GetRasterBand(1)->GetRasterDataType();
3568 93 : if (eDT == GDT_Int64)
3569 : {
3570 4 : nNoDataInt64 =
3571 4 : GetRasterBand(1)->GetNoDataValueAsInt64(&bHasNoDataAsInt64);
3572 4 : if (!bHasNoDataAsInt64)
3573 2 : nNoDataInt64 = 0;
3574 : }
3575 89 : else if (eDT == GDT_UInt64)
3576 : {
3577 4 : nNoDataUInt64 =
3578 4 : GetRasterBand(1)->GetNoDataValueAsUInt64(&bHasNoDataAsUInt64);
3579 4 : if (!bHasNoDataAsUInt64)
3580 2 : nNoDataUInt64 = 0;
3581 : }
3582 : else
3583 : {
3584 85 : dfNoData = GetRasterBand(1)->GetNoDataValue(&bHasNoData);
3585 85 : if (!bHasNoData)
3586 69 : dfNoData = 0;
3587 : }
3588 :
3589 93 : if (m_poExternalDS)
3590 : {
3591 : int nBlockXSize, nBlockYSize;
3592 18 : GetRasterBand(1)->GetBlockSize(&nBlockXSize, &nBlockYSize);
3593 18 : const int nDTSize = GDALGetDataTypeSizeBytes(eDT);
3594 18 : const int nBlockSizeBytes = nBlockXSize * nBlockYSize * nDTSize;
3595 18 : const int l_nBlocksPerColumn = DIV_ROUND_UP(nRasterYSize, nBlockYSize);
3596 :
3597 18 : if (nBands == 1 || EQUAL(m_osInterleave, "BSQ"))
3598 : {
3599 : // We need to make sure that blocks are written in the right order
3600 38 : for (int i = 0; i < nBands; i++)
3601 : {
3602 21 : if (eDT == GDT_Int64)
3603 : {
3604 1 : if (m_poExternalDS->GetRasterBand(i + 1)->RasterIO(
3605 : GF_Write, 0, 0, nRasterXSize, nRasterYSize,
3606 1 : &nNoDataInt64, 1, 1, eDT, 0, 0, nullptr) != CE_None)
3607 : {
3608 0 : return false;
3609 : }
3610 : }
3611 20 : else if (eDT == GDT_UInt64)
3612 : {
3613 1 : if (m_poExternalDS->GetRasterBand(i + 1)->RasterIO(
3614 : GF_Write, 0, 0, nRasterXSize, nRasterYSize,
3615 : &nNoDataUInt64, 1, 1, eDT, 0, 0,
3616 1 : nullptr) != CE_None)
3617 : {
3618 0 : return false;
3619 : }
3620 : }
3621 : else
3622 : {
3623 19 : if (m_poExternalDS->GetRasterBand(i + 1)->Fill(dfNoData) !=
3624 : CE_None)
3625 : {
3626 0 : return false;
3627 : }
3628 : }
3629 : }
3630 17 : m_poExternalDS->FlushCache(false);
3631 :
3632 : // Check that blocks are effectively written in expected order.
3633 17 : GIntBig nLastOffset = 0;
3634 38 : for (int i = 0; i < nBands; i++)
3635 : {
3636 333 : for (int y = 0; y < l_nBlocksPerColumn; y++)
3637 : {
3638 : const char *pszBlockOffset =
3639 624 : m_poExternalDS->GetRasterBand(i + 1)->GetMetadataItem(
3640 312 : CPLSPrintf("BLOCK_OFFSET_%d_%d", 0, y), "TIFF");
3641 312 : if (pszBlockOffset)
3642 : {
3643 312 : GIntBig nOffset = CPLAtoGIntBig(pszBlockOffset);
3644 312 : if (i != 0 || y != 0)
3645 : {
3646 295 : if (nOffset != nLastOffset + nBlockSizeBytes)
3647 : {
3648 0 : CPLError(CE_Warning, CPLE_AppDefined,
3649 : "Block %d,%d band %d not at expected "
3650 : "offset",
3651 : 0, y, i + 1);
3652 0 : return false;
3653 : }
3654 : }
3655 312 : nLastOffset = nOffset;
3656 : }
3657 : else
3658 : {
3659 0 : CPLError(CE_Warning, CPLE_AppDefined,
3660 : "Block %d,%d band %d not at expected "
3661 : "offset",
3662 : 0, y, i + 1);
3663 0 : return false;
3664 : }
3665 : }
3666 : }
3667 : }
3668 : else
3669 : {
3670 1 : void *pBlockData = VSI_MALLOC_VERBOSE(nBlockSizeBytes);
3671 1 : if (pBlockData == nullptr)
3672 0 : return false;
3673 1 : if (eDT == GDT_Int64)
3674 : {
3675 0 : GDALCopyWords(&nNoDataInt64, eDT, 0, pBlockData, eDT, nDTSize,
3676 : nBlockXSize * nBlockYSize);
3677 : }
3678 1 : else if (eDT == GDT_UInt64)
3679 : {
3680 0 : GDALCopyWords(&nNoDataUInt64, eDT, 0, pBlockData, eDT, nDTSize,
3681 : nBlockXSize * nBlockYSize);
3682 : }
3683 : else
3684 : {
3685 1 : GDALCopyWords(&dfNoData, GDT_Float64, 0, pBlockData, eDT,
3686 : nDTSize, nBlockXSize * nBlockYSize);
3687 : }
3688 2 : for (int y = 0; y < l_nBlocksPerColumn; y++)
3689 : {
3690 4 : for (int i = 0; i < nBands; i++)
3691 : {
3692 3 : if (m_poExternalDS->GetRasterBand(i + 1)->WriteBlock(
3693 3 : 0, y, pBlockData) != CE_None)
3694 : {
3695 0 : VSIFree(pBlockData);
3696 0 : return false;
3697 : }
3698 : }
3699 : }
3700 1 : VSIFree(pBlockData);
3701 1 : m_poExternalDS->FlushCache(false);
3702 :
3703 : // Check that blocks are effectively written in expected order.
3704 1 : GIntBig nLastOffset = 0;
3705 2 : for (int y = 0; y < l_nBlocksPerColumn; y++)
3706 : {
3707 : const char *pszBlockOffset =
3708 2 : m_poExternalDS->GetRasterBand(1)->GetMetadataItem(
3709 1 : CPLSPrintf("BLOCK_OFFSET_%d_%d", 0, y), "TIFF");
3710 1 : if (pszBlockOffset)
3711 : {
3712 1 : GIntBig nOffset = CPLAtoGIntBig(pszBlockOffset);
3713 1 : if (y != 0)
3714 : {
3715 0 : if (nOffset !=
3716 0 : nLastOffset +
3717 0 : static_cast<GIntBig>(nBlockSizeBytes) * nBands)
3718 : {
3719 0 : CPLError(CE_Warning, CPLE_AppDefined,
3720 : "Block %d,%d not at expected "
3721 : "offset",
3722 : 0, y);
3723 0 : return false;
3724 : }
3725 : }
3726 1 : nLastOffset = nOffset;
3727 : }
3728 : else
3729 : {
3730 0 : CPLError(CE_Warning, CPLE_AppDefined,
3731 : "Block %d,%d not at expected "
3732 : "offset",
3733 : 0, y);
3734 0 : return false;
3735 : }
3736 : }
3737 : }
3738 :
3739 18 : return true;
3740 : }
3741 :
3742 75 : const int nDTSize = GDALGetDataTypeSizeBytes(eDT);
3743 75 : const vsi_l_offset nFileSize = static_cast<vsi_l_offset>(nRasterXSize) *
3744 75 : nRasterYSize * nBands * nDTSize;
3745 75 : if ((eDT == GDT_Int64 && (nNoDataInt64 == 0 || !bHasNoDataAsInt64)) ||
3746 73 : (eDT == GDT_UInt64 && (nNoDataUInt64 == 0 || !bHasNoDataAsUInt64)) ||
3747 70 : (eDT != GDT_Int64 && eDT != GDT_UInt64 &&
3748 69 : (dfNoData == 0 || !bHasNoData)))
3749 : {
3750 66 : if (VSIFTruncateL(m_fpImage, nFileSize) != 0)
3751 : {
3752 0 : CPLError(CE_Failure, CPLE_FileIO,
3753 : "Cannot create file of size " CPL_FRMT_GUIB " bytes",
3754 : nFileSize);
3755 0 : return false;
3756 : }
3757 : }
3758 : else
3759 : {
3760 9 : size_t nLineSize = static_cast<size_t>(nRasterXSize) * nDTSize;
3761 9 : void *pData = VSI_MALLOC_VERBOSE(nLineSize);
3762 9 : if (pData == nullptr)
3763 0 : return false;
3764 9 : if (eDT == GDT_Int64)
3765 : {
3766 1 : GDALCopyWords(&nNoDataInt64, eDT, 0, pData, eDT, nDTSize,
3767 : nRasterXSize);
3768 : }
3769 8 : else if (eDT == GDT_UInt64)
3770 : {
3771 1 : GDALCopyWords(&nNoDataUInt64, eDT, 0, pData, eDT, nDTSize,
3772 : nRasterXSize);
3773 : }
3774 : else
3775 : {
3776 7 : GDALCopyWords(&dfNoData, GDT_Float64, 0, pData, eDT, nDTSize,
3777 : nRasterXSize);
3778 : }
3779 : #ifdef CPL_MSB
3780 : if (GDALDataTypeIsComplex(eDT))
3781 : {
3782 : GDALSwapWords(pData, nDTSize / 2, nRasterXSize * 2, nDTSize / 2);
3783 : }
3784 : else
3785 : {
3786 : GDALSwapWords(pData, nDTSize, nRasterXSize, nDTSize);
3787 : }
3788 : #endif
3789 18 : for (vsi_l_offset i = 0;
3790 18 : i < static_cast<vsi_l_offset>(nRasterYSize) * nBands; i++)
3791 : {
3792 9 : size_t nBytesWritten = VSIFWriteL(pData, 1, nLineSize, m_fpImage);
3793 9 : if (nBytesWritten != nLineSize)
3794 : {
3795 0 : CPLError(CE_Failure, CPLE_FileIO,
3796 : "Cannot create file of size " CPL_FRMT_GUIB " bytes",
3797 : nFileSize);
3798 0 : VSIFree(pData);
3799 0 : return false;
3800 : }
3801 : }
3802 9 : VSIFree(pData);
3803 : }
3804 75 : return true;
3805 : }
3806 :
3807 : /************************************************************************/
3808 : /* GetSpecialConstants() */
3809 : /************************************************************************/
3810 :
3811 12 : static CPLXMLNode *GetSpecialConstants(const CPLString &osPrefix,
3812 : CPLXMLNode *psFileAreaObservational)
3813 : {
3814 27 : for (CPLXMLNode *psIter = psFileAreaObservational->psChild; psIter;
3815 15 : psIter = psIter->psNext)
3816 : {
3817 48 : if (psIter->eType == CXT_Element &&
3818 48 : STARTS_WITH(psIter->pszValue, (osPrefix + "Array").c_str()))
3819 : {
3820 : CPLXMLNode *psSC =
3821 11 : CPLGetXMLNode(psIter, (osPrefix + "Special_Constants").c_str());
3822 11 : if (psSC)
3823 : {
3824 9 : CPLXMLNode *psNext = psSC->psNext;
3825 9 : psSC->psNext = nullptr;
3826 9 : CPLXMLNode *psRet = CPLCloneXMLTree(psSC);
3827 9 : psSC->psNext = psNext;
3828 9 : return psRet;
3829 : }
3830 : }
3831 : }
3832 3 : return nullptr;
3833 : }
3834 :
3835 : /************************************************************************/
3836 : /* WriteHeaderAppendCase() */
3837 : /************************************************************************/
3838 :
3839 4 : void PDS4Dataset::WriteHeaderAppendCase()
3840 : {
3841 4 : CPLXMLTreeCloser oCloser(CPLParseXMLFile(GetDescription()));
3842 4 : CPLXMLNode *psRoot = oCloser.get();
3843 4 : if (psRoot == nullptr)
3844 0 : return;
3845 4 : CPLString osPrefix;
3846 4 : CPLXMLNode *psProduct = CPLGetXMLNode(psRoot, "=Product_Observational");
3847 4 : if (psProduct == nullptr)
3848 : {
3849 0 : psProduct = CPLGetXMLNode(psRoot, "=pds:Product_Observational");
3850 0 : if (psProduct)
3851 0 : osPrefix = "pds:";
3852 : }
3853 4 : if (psProduct == nullptr)
3854 : {
3855 0 : CPLError(CE_Failure, CPLE_AppDefined,
3856 : "Cannot find Product_Observational element");
3857 0 : return;
3858 : }
3859 4 : CPLXMLNode *psFAO = CPLGetXMLNode(
3860 8 : psProduct, (osPrefix + "File_Area_Observational").c_str());
3861 4 : if (psFAO == nullptr)
3862 : {
3863 0 : CPLError(CE_Failure, CPLE_AppDefined,
3864 : "Cannot find File_Area_Observational element");
3865 0 : return;
3866 : }
3867 :
3868 4 : WriteArray(osPrefix, psFAO, nullptr, nullptr);
3869 :
3870 4 : CPLSerializeXMLTreeToFile(psRoot, GetDescription());
3871 : }
3872 :
3873 : /************************************************************************/
3874 : /* WriteArray() */
3875 : /************************************************************************/
3876 :
3877 135 : void PDS4Dataset::WriteArray(const CPLString &osPrefix, CPLXMLNode *psFAO,
3878 : const char *pszLocalIdentifierDefault,
3879 : CPLXMLNode *psTemplateSpecialConstants)
3880 : {
3881 270 : const char *pszArrayType = CSLFetchNameValueDef(
3882 135 : m_papszCreationOptions, "ARRAY_TYPE", "Array_3D_Image");
3883 135 : const bool bIsArray2D = STARTS_WITH(pszArrayType, "Array_2D");
3884 : CPLXMLNode *psArray =
3885 135 : CPLCreateXMLNode(psFAO, CXT_Element, (osPrefix + pszArrayType).c_str());
3886 :
3887 270 : const char *pszLocalIdentifier = CSLFetchNameValueDef(
3888 135 : m_papszCreationOptions, "ARRAY_IDENTIFIER", pszLocalIdentifierDefault);
3889 135 : if (pszLocalIdentifier)
3890 : {
3891 131 : CPLCreateXMLElementAndValue(psArray,
3892 262 : (osPrefix + "local_identifier").c_str(),
3893 : pszLocalIdentifier);
3894 : }
3895 :
3896 135 : GUIntBig nOffset = m_nBaseOffset;
3897 135 : if (m_poExternalDS)
3898 : {
3899 : const char *pszOffset =
3900 18 : m_poExternalDS->GetRasterBand(1)->GetMetadataItem(
3901 18 : "BLOCK_OFFSET_0_0", "TIFF");
3902 18 : if (pszOffset)
3903 18 : nOffset = CPLAtoGIntBig(pszOffset);
3904 : }
3905 270 : CPLAddXMLAttributeAndValue(
3906 270 : CPLCreateXMLElementAndValue(psArray, (osPrefix + "offset").c_str(),
3907 : CPLSPrintf(CPL_FRMT_GUIB, nOffset)),
3908 : "unit", "byte");
3909 135 : CPLCreateXMLElementAndValue(psArray, (osPrefix + "axes").c_str(),
3910 : (bIsArray2D) ? "2" : "3");
3911 135 : CPLCreateXMLElementAndValue(
3912 270 : psArray, (osPrefix + "axis_index_order").c_str(), "Last Index Fastest");
3913 135 : CPLXMLNode *psElementArray = CPLCreateXMLNode(
3914 270 : psArray, CXT_Element, (osPrefix + "Element_Array").c_str());
3915 135 : GDALDataType eDT = GetRasterBand(1)->GetRasterDataType();
3916 135 : const char *pszDataType =
3917 187 : (eDT == GDT_UInt8) ? "UnsignedByte"
3918 101 : : (eDT == GDT_Int8) ? "SignedByte"
3919 91 : : (eDT == GDT_UInt16) ? "UnsignedLSB2"
3920 79 : : (eDT == GDT_Int16) ? (m_bIsLSB ? "SignedLSB2" : "SignedMSB2")
3921 70 : : (eDT == GDT_UInt32) ? (m_bIsLSB ? "UnsignedLSB4" : "UnsignedMSB4")
3922 62 : : (eDT == GDT_Int32) ? (m_bIsLSB ? "SignedLSB4" : "SignedMSB4")
3923 54 : : (eDT == GDT_UInt64) ? (m_bIsLSB ? "UnsignedLSB8" : "UnsignedMSB8")
3924 46 : : (eDT == GDT_Int64) ? (m_bIsLSB ? "SignedLSB8" : "SignedMSB8")
3925 : : (eDT == GDT_Float32)
3926 35 : ? (m_bIsLSB ? "IEEE754LSBSingle" : "IEEE754MSBSingle")
3927 : : (eDT == GDT_Float64)
3928 22 : ? (m_bIsLSB ? "IEEE754LSBDouble" : "IEEE754MSBDouble")
3929 12 : : (eDT == GDT_CFloat32) ? (m_bIsLSB ? "ComplexLSB8" : "ComplexMSB8")
3930 4 : : (eDT == GDT_CFloat64) ? (m_bIsLSB ? "ComplexLSB16" : "ComplexMSB16")
3931 : : "should not happen";
3932 135 : CPLCreateXMLElementAndValue(psElementArray,
3933 270 : (osPrefix + "data_type").c_str(), pszDataType);
3934 :
3935 135 : const char *pszUnits = GetRasterBand(1)->GetUnitType();
3936 135 : const char *pszUnitsCO = CSLFetchNameValue(m_papszCreationOptions, "UNIT");
3937 135 : if (pszUnitsCO)
3938 : {
3939 0 : pszUnits = pszUnitsCO;
3940 : }
3941 135 : if (pszUnits && pszUnits[0] != 0)
3942 : {
3943 1 : CPLCreateXMLElementAndValue(psElementArray, (osPrefix + "unit").c_str(),
3944 : pszUnits);
3945 : }
3946 :
3947 135 : int bHasScale = FALSE;
3948 135 : double dfScale = GetRasterBand(1)->GetScale(&bHasScale);
3949 135 : if (bHasScale && dfScale != 1.0)
3950 : {
3951 10 : CPLCreateXMLElementAndValue(psElementArray,
3952 10 : (osPrefix + "scaling_factor").c_str(),
3953 : CPLSPrintf("%.17g", dfScale));
3954 : }
3955 :
3956 135 : int bHasOffset = FALSE;
3957 135 : double dfOffset = GetRasterBand(1)->GetOffset(&bHasOffset);
3958 135 : if (bHasOffset && dfOffset != 1.0)
3959 : {
3960 10 : CPLCreateXMLElementAndValue(psElementArray,
3961 10 : (osPrefix + "value_offset").c_str(),
3962 : CPLSPrintf("%.17g", dfOffset));
3963 : }
3964 :
3965 : // Axis definitions
3966 : {
3967 135 : CPLXMLNode *psAxis = CPLCreateXMLNode(
3968 270 : psArray, CXT_Element, (osPrefix + "Axis_Array").c_str());
3969 270 : CPLCreateXMLElementAndValue(
3970 270 : psAxis, (osPrefix + "axis_name").c_str(),
3971 135 : EQUAL(m_osInterleave, "BSQ")
3972 : ? "Band"
3973 : :
3974 : /* EQUAL(m_osInterleave, "BIL") ? "Line" : */
3975 : "Line");
3976 270 : CPLCreateXMLElementAndValue(
3977 270 : psAxis, (osPrefix + "elements").c_str(),
3978 : CPLSPrintf("%d",
3979 135 : EQUAL(m_osInterleave, "BSQ")
3980 : ? nBands
3981 : :
3982 : /* EQUAL(m_osInterleave, "BIL") ? nRasterYSize : */
3983 : nRasterYSize));
3984 135 : CPLCreateXMLElementAndValue(
3985 270 : psAxis, (osPrefix + "sequence_number").c_str(), "1");
3986 : }
3987 : {
3988 135 : CPLXMLNode *psAxis = CPLCreateXMLNode(
3989 270 : psArray, CXT_Element, (osPrefix + "Axis_Array").c_str());
3990 135 : CPLCreateXMLElementAndValue(psAxis, (osPrefix + "axis_name").c_str(),
3991 135 : EQUAL(m_osInterleave, "BSQ") ? "Line"
3992 11 : : EQUAL(m_osInterleave, "BIL") ? "Band"
3993 : : "Sample");
3994 270 : CPLCreateXMLElementAndValue(
3995 270 : psAxis, (osPrefix + "elements").c_str(),
3996 135 : CPLSPrintf("%d", EQUAL(m_osInterleave, "BSQ") ? nRasterYSize
3997 11 : : EQUAL(m_osInterleave, "BIL") ? nBands
3998 : : nRasterXSize));
3999 135 : CPLCreateXMLElementAndValue(
4000 270 : psAxis, (osPrefix + "sequence_number").c_str(), "2");
4001 : }
4002 135 : if (!bIsArray2D)
4003 : {
4004 134 : CPLXMLNode *psAxis = CPLCreateXMLNode(
4005 268 : psArray, CXT_Element, (osPrefix + "Axis_Array").c_str());
4006 134 : CPLCreateXMLElementAndValue(psAxis, (osPrefix + "axis_name").c_str(),
4007 134 : EQUAL(m_osInterleave, "BSQ") ? "Sample"
4008 10 : : EQUAL(m_osInterleave, "BIL") ? "Sample"
4009 : : "Band");
4010 268 : CPLCreateXMLElementAndValue(
4011 268 : psAxis, (osPrefix + "elements").c_str(),
4012 134 : CPLSPrintf("%d", EQUAL(m_osInterleave, "BSQ") ? nRasterXSize
4013 10 : : EQUAL(m_osInterleave, "BIL") ? nRasterXSize
4014 : : nBands));
4015 134 : CPLCreateXMLElementAndValue(
4016 268 : psAxis, (osPrefix + "sequence_number").c_str(), "3");
4017 : }
4018 :
4019 135 : int bHasNoData = FALSE;
4020 135 : int64_t nNoDataInt64 = 0;
4021 135 : uint64_t nNoDataUInt64 = 0;
4022 135 : double dfNoData = 0;
4023 135 : if (eDT == GDT_Int64)
4024 4 : nNoDataInt64 = GetRasterBand(1)->GetNoDataValueAsInt64(&bHasNoData);
4025 131 : else if (eDT == GDT_UInt64)
4026 4 : nNoDataUInt64 = GetRasterBand(1)->GetNoDataValueAsUInt64(&bHasNoData);
4027 : else
4028 127 : dfNoData = GetRasterBand(1)->GetNoDataValue(&bHasNoData);
4029 :
4030 270 : std::string osNoData;
4031 135 : if (bHasNoData)
4032 : {
4033 29 : if (eDT == GDT_Int64)
4034 2 : osNoData = std::to_string(nNoDataInt64);
4035 27 : else if (eDT == GDT_UInt64)
4036 2 : osNoData = std::to_string(nNoDataUInt64);
4037 25 : else if (GDALDataTypeIsInteger(GetRasterBand(1)->GetRasterDataType()))
4038 20 : osNoData = std::to_string(static_cast<int64_t>(dfNoData));
4039 5 : else if (eDT == GDT_Float32)
4040 : {
4041 : uint32_t nVal;
4042 3 : float fVal = static_cast<float>(dfNoData);
4043 3 : memcpy(&nVal, &fVal, sizeof(nVal));
4044 3 : osNoData = CPLSPrintf("0x%08X", nVal);
4045 : }
4046 : else
4047 : {
4048 : uint64_t nVal;
4049 2 : memcpy(&nVal, &dfNoData, sizeof(nVal));
4050 : osNoData =
4051 2 : CPLSPrintf("0x%16llX", static_cast<unsigned long long>(nVal));
4052 : }
4053 : }
4054 :
4055 135 : if (psTemplateSpecialConstants)
4056 : {
4057 9 : CPLAddXMLChild(psArray, psTemplateSpecialConstants);
4058 9 : if (bHasNoData)
4059 : {
4060 : CPLXMLNode *psMC =
4061 6 : CPLGetXMLNode(psTemplateSpecialConstants,
4062 12 : (osPrefix + "missing_constant").c_str());
4063 6 : if (psMC != nullptr)
4064 : {
4065 4 : if (psMC->psChild && psMC->psChild->eType == CXT_Text)
4066 : {
4067 4 : CPLFree(psMC->psChild->pszValue);
4068 4 : psMC->psChild->pszValue = CPLStrdup(osNoData.c_str());
4069 : }
4070 : }
4071 : else
4072 : {
4073 : CPLXMLNode *psSaturatedConstant =
4074 2 : CPLGetXMLNode(psTemplateSpecialConstants,
4075 4 : (osPrefix + "saturated_constant").c_str());
4076 4 : psMC = CPLCreateXMLElementAndValue(
4077 4 : nullptr, (osPrefix + "missing_constant").c_str(),
4078 : osNoData.c_str());
4079 : CPLXMLNode *psNext;
4080 2 : if (psSaturatedConstant)
4081 : {
4082 1 : psNext = psSaturatedConstant->psNext;
4083 1 : psSaturatedConstant->psNext = psMC;
4084 : }
4085 : else
4086 : {
4087 1 : psNext = psTemplateSpecialConstants->psChild;
4088 1 : psTemplateSpecialConstants->psChild = psMC;
4089 : }
4090 2 : psMC->psNext = psNext;
4091 : }
4092 : }
4093 : }
4094 126 : else if (bHasNoData)
4095 : {
4096 23 : CPLXMLNode *psSC = CPLCreateXMLNode(
4097 46 : psArray, CXT_Element, (osPrefix + "Special_Constants").c_str());
4098 46 : CPLCreateXMLElementAndValue(
4099 46 : psSC, (osPrefix + "missing_constant").c_str(), osNoData.c_str());
4100 : }
4101 135 : }
4102 :
4103 : /************************************************************************/
4104 : /* WriteVectorLayers() */
4105 : /************************************************************************/
4106 :
4107 188 : void PDS4Dataset::WriteVectorLayers(CPLXMLNode *psProduct)
4108 : {
4109 376 : CPLString osPrefix;
4110 188 : if (STARTS_WITH(psProduct->pszValue, "pds:"))
4111 0 : osPrefix = "pds:";
4112 :
4113 260 : for (auto &poLayer : m_apoLayers)
4114 : {
4115 72 : if (!poLayer->IsDirtyHeader())
4116 3 : continue;
4117 :
4118 69 : if (poLayer->GetFeatureCount(false) == 0)
4119 : {
4120 16 : CPLError(CE_Warning, CPLE_AppDefined,
4121 : "Writing header for layer %s which has 0 features. "
4122 : "This is not legal in PDS4",
4123 16 : poLayer->GetName());
4124 : }
4125 :
4126 69 : if (poLayer->GetRawFieldCount() == 0)
4127 : {
4128 16 : CPLError(CE_Warning, CPLE_AppDefined,
4129 : "Writing header for layer %s which has 0 fields. "
4130 : "This is not legal in PDS4",
4131 16 : poLayer->GetName());
4132 : }
4133 :
4134 : const std::string osRelativePath(
4135 138 : CPLExtractRelativePath(CPLGetPathSafe(m_osXMLFilename).c_str(),
4136 207 : poLayer->GetFileName(), nullptr));
4137 :
4138 69 : bool bFound = false;
4139 634 : for (CPLXMLNode *psIter = psProduct->psChild; psIter != nullptr;
4140 565 : psIter = psIter->psNext)
4141 : {
4142 733 : if (psIter->eType == CXT_Element &&
4143 164 : strcmp(psIter->pszValue,
4144 733 : (osPrefix + "File_Area_Observational").c_str()) == 0)
4145 : {
4146 23 : const char *pszFilename = CPLGetXMLValue(
4147 : psIter,
4148 46 : (osPrefix + "File." + osPrefix + "file_name").c_str(), "");
4149 23 : if (strcmp(pszFilename, osRelativePath.c_str()) == 0)
4150 : {
4151 4 : poLayer->RefreshFileAreaObservational(psIter);
4152 4 : bFound = true;
4153 4 : break;
4154 : }
4155 : }
4156 : }
4157 69 : if (!bFound)
4158 : {
4159 65 : CPLXMLNode *psFAO = CPLCreateXMLNode(
4160 : psProduct, CXT_Element,
4161 130 : (osPrefix + "File_Area_Observational").c_str());
4162 65 : CPLXMLNode *psFile = CPLCreateXMLNode(psFAO, CXT_Element,
4163 130 : (osPrefix + "File").c_str());
4164 130 : CPLCreateXMLElementAndValue(psFile,
4165 130 : (osPrefix + "file_name").c_str(),
4166 : osRelativePath.c_str());
4167 65 : poLayer->RefreshFileAreaObservational(psFAO);
4168 : }
4169 : }
4170 188 : }
4171 :
4172 : /************************************************************************/
4173 : /* CreateHeader() */
4174 : /************************************************************************/
4175 :
4176 181 : void PDS4Dataset::CreateHeader(CPLXMLNode *psProduct,
4177 : const char *pszCARTVersion)
4178 : {
4179 181 : CPLString osPrefix;
4180 181 : if (STARTS_WITH(psProduct->pszValue, "pds:"))
4181 0 : osPrefix = "pds:";
4182 :
4183 181 : if (m_oSRS.IsEmpty() && GetLayerCount() >= 1)
4184 : {
4185 46 : if (const auto poSRS = GetLayer(0)->GetSpatialRef())
4186 2 : m_oSRS = *poSRS;
4187 : }
4188 :
4189 280 : if (!m_oSRS.IsEmpty() &&
4190 99 : CSLFetchNameValue(m_papszCreationOptions, "VAR_TARGET") == nullptr)
4191 : {
4192 98 : const char *pszTarget = nullptr;
4193 98 : if (fabs(m_oSRS.GetSemiMajor() - 6378137) < 0.001 * 6378137)
4194 : {
4195 77 : pszTarget = "Earth";
4196 77 : m_papszCreationOptions = CSLSetNameValue(
4197 : m_papszCreationOptions, "VAR_TARGET_TYPE", "Planet");
4198 : }
4199 : else
4200 : {
4201 21 : const char *pszDatum = m_oSRS.GetAttrValue("DATUM");
4202 21 : if (pszDatum && STARTS_WITH(pszDatum, "D_"))
4203 : {
4204 3 : pszTarget = pszDatum + 2;
4205 : }
4206 18 : else if (pszDatum)
4207 : {
4208 18 : pszTarget = pszDatum;
4209 : }
4210 : }
4211 98 : if (pszTarget)
4212 : {
4213 98 : m_papszCreationOptions = CSLSetNameValue(m_papszCreationOptions,
4214 : "VAR_TARGET", pszTarget);
4215 : }
4216 : }
4217 181 : SubstituteVariables(psProduct, m_papszCreationOptions);
4218 :
4219 : // Remove <Discipline_Area>/<disp:Display_Settings> if there is no raster
4220 181 : if (GetRasterCount() == 0)
4221 : {
4222 : CPLXMLNode *psDisciplineArea =
4223 141 : CPLGetXMLNode(psProduct, (osPrefix + "Observation_Area." +
4224 94 : osPrefix + "Discipline_Area")
4225 : .c_str());
4226 47 : if (psDisciplineArea)
4227 : {
4228 : CPLXMLNode *psDisplaySettings =
4229 47 : CPLGetXMLNode(psDisciplineArea, "disp:Display_Settings");
4230 47 : if (psDisplaySettings)
4231 : {
4232 47 : CPLRemoveXMLChild(psDisciplineArea, psDisplaySettings);
4233 47 : CPLDestroyXMLNode(psDisplaySettings);
4234 : }
4235 : }
4236 : }
4237 :
4238 : // Depending on the version of the DISP schema, Local_Internal_Reference
4239 : // may be in the disp: namespace or the default one.
4240 : const auto GetLocalIdentifierReferenceFromDisciplineArea =
4241 224 : [](const CPLXMLNode *psDisciplineArea, const char *pszDefault)
4242 : {
4243 224 : return CPLGetXMLValue(
4244 : psDisciplineArea,
4245 : "disp:Display_Settings.Local_Internal_Reference."
4246 : "local_identifier_reference",
4247 : CPLGetXMLValue(
4248 : psDisciplineArea,
4249 : "disp:Display_Settings.disp:Local_Internal_Reference."
4250 : "local_identifier_reference",
4251 224 : pszDefault));
4252 : };
4253 :
4254 181 : if (GetRasterCount() || !m_oSRS.IsEmpty())
4255 : {
4256 : CPLXMLNode *psDisciplineArea =
4257 408 : CPLGetXMLNode(psProduct, (osPrefix + "Observation_Area." +
4258 272 : osPrefix + "Discipline_Area")
4259 : .c_str());
4260 136 : if (GetRasterCount() && !(m_bGotTransform && !m_oSRS.IsEmpty()))
4261 : {
4262 : // if we have no georeferencing, strip any existing georeferencing
4263 : // from the template
4264 37 : if (psDisciplineArea)
4265 : {
4266 : CPLXMLNode *psCart =
4267 33 : CPLGetXMLNode(psDisciplineArea, "cart:Cartography");
4268 33 : if (psCart == nullptr)
4269 32 : psCart = CPLGetXMLNode(psDisciplineArea, "Cartography");
4270 33 : if (psCart)
4271 : {
4272 1 : CPLRemoveXMLChild(psDisciplineArea, psCart);
4273 1 : CPLDestroyXMLNode(psCart);
4274 : }
4275 :
4276 33 : if (CPLGetXMLNode(psDisciplineArea,
4277 33 : "sp:Spectral_Characteristics"))
4278 : {
4279 : const char *pszArrayType =
4280 1 : CSLFetchNameValue(m_papszCreationOptions, "ARRAY_TYPE");
4281 : // The schematron PDS4_SP_1100.sch requires that
4282 : // sp:local_identifier_reference is used by
4283 : // Array_[2D|3D]_Spectrum/pds:local_identifier
4284 1 : if (pszArrayType == nullptr)
4285 : {
4286 1 : m_papszCreationOptions =
4287 1 : CSLSetNameValue(m_papszCreationOptions,
4288 : "ARRAY_TYPE", "Array_3D_Spectrum");
4289 : }
4290 0 : else if (!EQUAL(pszArrayType, "Array_2D_Spectrum") &&
4291 0 : !EQUAL(pszArrayType, "Array_3D_Spectrum"))
4292 : {
4293 0 : CPLError(CE_Warning, CPLE_AppDefined,
4294 : "PDS4_SP_xxxx.sch schematron requires the "
4295 : "use of ARRAY_TYPE=Array_2D_Spectrum or "
4296 : "Array_3D_Spectrum");
4297 : }
4298 : }
4299 : }
4300 : }
4301 : else
4302 : {
4303 99 : if (psDisciplineArea == nullptr)
4304 : {
4305 2 : CPLXMLNode *psTI = CPLGetXMLNode(
4306 4 : psProduct, (osPrefix + "Observation_Area." + osPrefix +
4307 : "Target_Identification")
4308 : .c_str());
4309 2 : if (psTI == nullptr)
4310 : {
4311 1 : CPLError(CE_Failure, CPLE_AppDefined,
4312 : "Cannot find Target_Identification element in "
4313 : "template");
4314 1 : return;
4315 : }
4316 : psDisciplineArea =
4317 1 : CPLCreateXMLNode(nullptr, CXT_Element,
4318 2 : (osPrefix + "Discipline_Area").c_str());
4319 1 : if (psTI->psNext)
4320 0 : psDisciplineArea->psNext = psTI->psNext;
4321 1 : psTI->psNext = psDisciplineArea;
4322 : }
4323 : CPLXMLNode *psCart =
4324 98 : CPLGetXMLNode(psDisciplineArea, "cart:Cartography");
4325 98 : if (psCart == nullptr)
4326 93 : psCart = CPLGetXMLNode(psDisciplineArea, "Cartography");
4327 98 : if (psCart == nullptr)
4328 : {
4329 93 : psCart = CPLCreateXMLNode(psDisciplineArea, CXT_Element,
4330 : "cart:Cartography");
4331 93 : if (CPLGetXMLNode(psProduct, "xmlns:cart") == nullptr)
4332 : {
4333 : CPLXMLNode *psNS =
4334 1 : CPLCreateXMLNode(nullptr, CXT_Attribute, "xmlns:cart");
4335 1 : CPLCreateXMLNode(psNS, CXT_Text,
4336 : "http://pds.nasa.gov/pds4/cart/v1");
4337 1 : CPLAddXMLChild(psProduct, psNS);
4338 : CPLXMLNode *psSchemaLoc =
4339 1 : CPLGetXMLNode(psProduct, "xsi:schemaLocation");
4340 1 : if (psSchemaLoc != nullptr &&
4341 1 : psSchemaLoc->psChild != nullptr &&
4342 1 : psSchemaLoc->psChild->pszValue != nullptr)
4343 : {
4344 2 : CPLString osCartSchema;
4345 1 : if (strstr(psSchemaLoc->psChild->pszValue,
4346 : "PDS4_PDS_1800.xsd"))
4347 : {
4348 : // GDAL 2.4
4349 : osCartSchema = "https://pds.nasa.gov/pds4/cart/v1/"
4350 1 : "PDS4_CART_1700.xsd";
4351 1 : pszCARTVersion = "1700";
4352 : }
4353 0 : else if (strstr(psSchemaLoc->psChild->pszValue,
4354 : "PDS4_PDS_1B00.xsd"))
4355 : {
4356 : // GDAL 3.0
4357 : osCartSchema =
4358 : "https://raw.githubusercontent.com/"
4359 : "nasa-pds-data-dictionaries/ldd-cart/master/"
4360 0 : "build/1.B.0.0/PDS4_CART_1B00.xsd";
4361 0 : pszCARTVersion = "1B00";
4362 : }
4363 0 : else if (strstr(psSchemaLoc->psChild->pszValue,
4364 : "PDS4_PDS_1D00.xsd"))
4365 : {
4366 : // GDAL 3.1
4367 : osCartSchema = "https://pds.nasa.gov/pds4/cart/v1/"
4368 0 : "PDS4_CART_1D00_1933.xsd";
4369 0 : pszCARTVersion = "1D00_1933";
4370 : }
4371 0 : else if (strstr(psSchemaLoc->psChild->pszValue,
4372 : "PDS4_PDS_1G00_1950.xsd"))
4373 : {
4374 : // GDAL 3.4
4375 : osCartSchema = "https://pds.nasa.gov/pds4/cart/v1/"
4376 0 : "PDS4_CART_1G00_1950.xsd";
4377 0 : pszCARTVersion = "1G00_1950";
4378 : }
4379 : else
4380 : {
4381 : // GDAL 3.12
4382 : osCartSchema =
4383 : "https://pds.nasa.gov/pds4/cart/v1/"
4384 0 : "PDS4_CART_" CURRENT_CART_VERSION ".xsd";
4385 0 : pszCARTVersion = CURRENT_CART_VERSION;
4386 : }
4387 1 : CPLString osNewVal(psSchemaLoc->psChild->pszValue);
4388 : osNewVal +=
4389 1 : " http://pds.nasa.gov/pds4/cart/v1 " + osCartSchema;
4390 1 : CPLFree(psSchemaLoc->psChild->pszValue);
4391 1 : psSchemaLoc->psChild->pszValue = CPLStrdup(osNewVal);
4392 : }
4393 : }
4394 : }
4395 : else
4396 : {
4397 5 : if (psCart->psChild)
4398 : {
4399 5 : CPLDestroyXMLNode(psCart->psChild);
4400 5 : psCart->psChild = nullptr;
4401 : }
4402 : }
4403 :
4404 98 : if (IsCARTVersionGTE(pszCARTVersion, "1900"))
4405 : {
4406 : const char *pszLocalIdentifier =
4407 93 : GetLocalIdentifierReferenceFromDisciplineArea(
4408 : psDisciplineArea,
4409 93 : GetRasterCount() == 0 && GetLayerCount() > 0
4410 2 : ? GetLayer(0)->GetName()
4411 : : "image");
4412 93 : CPLXMLNode *psLIR = CPLCreateXMLNode(
4413 : psCart, CXT_Element,
4414 186 : (osPrefix + "Local_Internal_Reference").c_str());
4415 93 : CPLCreateXMLElementAndValue(
4416 186 : psLIR, (osPrefix + "local_identifier_reference").c_str(),
4417 : pszLocalIdentifier);
4418 93 : CPLCreateXMLElementAndValue(
4419 186 : psLIR, (osPrefix + "local_reference_type").c_str(),
4420 : "cartography_parameters_to_image_object");
4421 : }
4422 :
4423 98 : WriteGeoreferencing(psCart, pszCARTVersion);
4424 : }
4425 :
4426 270 : const char *pszVertDir = CSLFetchNameValue(
4427 135 : m_papszCreationOptions, "VAR_VERTICAL_DISPLAY_DIRECTION");
4428 135 : if (pszVertDir)
4429 : {
4430 : CPLXMLNode *psVertDirNode =
4431 1 : CPLGetXMLNode(psDisciplineArea,
4432 : "disp:Display_Settings.disp:Display_Direction."
4433 : "disp:vertical_display_direction");
4434 1 : if (psVertDirNode == nullptr)
4435 : {
4436 0 : CPLError(
4437 : CE_Warning, CPLE_AppDefined,
4438 : "PDS4 template lacks a disp:vertical_display_direction "
4439 : "element where to write %s",
4440 : pszVertDir);
4441 : }
4442 : else
4443 : {
4444 1 : CPLDestroyXMLNode(psVertDirNode->psChild);
4445 1 : psVertDirNode->psChild =
4446 1 : CPLCreateXMLNode(nullptr, CXT_Text, pszVertDir);
4447 : }
4448 : }
4449 : }
4450 : else
4451 : {
4452 : // Remove Observation_Area.Discipline_Area if it contains only
4453 : // <disp:Display_Settings> or is empty
4454 : CPLXMLNode *psObservationArea =
4455 45 : CPLGetXMLNode(psProduct, (osPrefix + "Observation_Area").c_str());
4456 45 : if (psObservationArea)
4457 : {
4458 45 : CPLXMLNode *psDisciplineArea = CPLGetXMLNode(
4459 90 : psObservationArea, (osPrefix + "Discipline_Area").c_str());
4460 45 : if (psDisciplineArea &&
4461 45 : (psDisciplineArea->psChild == nullptr ||
4462 0 : (psDisciplineArea->psChild->eType == CXT_Element &&
4463 0 : psDisciplineArea->psChild->psNext == nullptr &&
4464 0 : strcmp(psDisciplineArea->psChild->pszValue,
4465 : "disp:Display_Settings") == 0)))
4466 : {
4467 45 : CPLRemoveXMLChild(psObservationArea, psDisciplineArea);
4468 45 : CPLDestroyXMLNode(psDisciplineArea);
4469 : }
4470 : }
4471 : }
4472 :
4473 180 : if (m_bStripFileAreaObservationalFromTemplate)
4474 : {
4475 180 : m_bStripFileAreaObservationalFromTemplate = false;
4476 180 : CPLXMLNode *psObservationArea = nullptr;
4477 180 : CPLXMLNode *psPrev = nullptr;
4478 180 : CPLXMLNode *psTemplateSpecialConstants = nullptr;
4479 1618 : for (CPLXMLNode *psIter = psProduct->psChild; psIter != nullptr;)
4480 : {
4481 1809 : if (psIter->eType == CXT_Element &&
4482 1809 : psIter->pszValue == osPrefix + "Observation_Area")
4483 : {
4484 179 : psObservationArea = psIter;
4485 179 : psPrev = psIter;
4486 179 : psIter = psIter->psNext;
4487 : }
4488 1631 : else if (psIter->eType == CXT_Element &&
4489 192 : (psIter->pszValue ==
4490 1631 : osPrefix + "File_Area_Observational" ||
4491 180 : psIter->pszValue ==
4492 1439 : osPrefix + "File_Area_Observational_Supplemental"))
4493 : {
4494 12 : if (psIter->pszValue == osPrefix + "File_Area_Observational")
4495 : {
4496 : psTemplateSpecialConstants =
4497 12 : GetSpecialConstants(osPrefix, psIter);
4498 : }
4499 12 : if (psPrev)
4500 12 : psPrev->psNext = psIter->psNext;
4501 : else
4502 : {
4503 0 : CPLAssert(psProduct->psChild == psIter);
4504 0 : psProduct->psChild = psIter->psNext;
4505 : }
4506 12 : CPLXMLNode *psNext = psIter->psNext;
4507 12 : psIter->psNext = nullptr;
4508 12 : CPLDestroyXMLNode(psIter);
4509 12 : psIter = psNext;
4510 : }
4511 : else
4512 : {
4513 1247 : psPrev = psIter;
4514 1247 : psIter = psIter->psNext;
4515 : }
4516 : }
4517 180 : if (psObservationArea == nullptr)
4518 : {
4519 1 : CPLError(CE_Failure, CPLE_AppDefined,
4520 : "Cannot find Observation_Area in template");
4521 1 : CPLDestroyXMLNode(psTemplateSpecialConstants);
4522 1 : return;
4523 : }
4524 :
4525 179 : if (GetRasterCount())
4526 : {
4527 132 : CPLXMLNode *psFAOPrev = psObservationArea;
4528 134 : while (psFAOPrev->psNext != nullptr &&
4529 4 : psFAOPrev->psNext->eType == CXT_Comment)
4530 : {
4531 2 : psFAOPrev = psFAOPrev->psNext;
4532 : }
4533 132 : if (psFAOPrev->psNext != nullptr)
4534 : {
4535 : // There may be an optional Reference_List element between
4536 : // Observation_Area and File_Area_Observational
4537 4 : if (!(psFAOPrev->psNext->eType == CXT_Element &&
4538 2 : psFAOPrev->psNext->pszValue ==
4539 4 : osPrefix + "Reference_List"))
4540 : {
4541 1 : CPLError(CE_Failure, CPLE_AppDefined,
4542 : "Unexpected content found after Observation_Area "
4543 : "in template");
4544 1 : CPLDestroyXMLNode(psTemplateSpecialConstants);
4545 1 : return;
4546 : }
4547 1 : psFAOPrev = psFAOPrev->psNext;
4548 2 : while (psFAOPrev->psNext != nullptr &&
4549 1 : psFAOPrev->psNext->eType == CXT_Comment)
4550 : {
4551 1 : psFAOPrev = psFAOPrev->psNext;
4552 : }
4553 1 : if (psFAOPrev->psNext != nullptr)
4554 : {
4555 0 : CPLError(CE_Failure, CPLE_AppDefined,
4556 : "Unexpected content found after Reference_List in "
4557 : "template");
4558 0 : CPLDestroyXMLNode(psTemplateSpecialConstants);
4559 0 : return;
4560 : }
4561 : }
4562 :
4563 131 : CPLXMLNode *psFAO = CPLCreateXMLNode(
4564 : nullptr, CXT_Element,
4565 262 : (osPrefix + "File_Area_Observational").c_str());
4566 131 : psFAOPrev->psNext = psFAO;
4567 :
4568 131 : CPLXMLNode *psFile = CPLCreateXMLNode(psFAO, CXT_Element,
4569 262 : (osPrefix + "File").c_str());
4570 262 : CPLCreateXMLElementAndValue(psFile,
4571 262 : (osPrefix + "file_name").c_str(),
4572 : CPLGetFilename(m_osImageFilename));
4573 131 : if (m_bCreatedFromExistingBinaryFile)
4574 : {
4575 7 : CPLCreateXMLNode(psFile, CXT_Comment, PREEXISTING_BINARY_FILE);
4576 : }
4577 : CPLXMLNode *psDisciplineArea =
4578 393 : CPLGetXMLNode(psProduct, (osPrefix + "Observation_Area." +
4579 262 : osPrefix + "Discipline_Area")
4580 : .c_str());
4581 : const char *pszLocalIdentifier =
4582 131 : GetLocalIdentifierReferenceFromDisciplineArea(psDisciplineArea,
4583 : "image");
4584 :
4585 147 : if (m_poExternalDS && m_poExternalDS->GetDriver() &&
4586 16 : EQUAL(m_poExternalDS->GetDriver()->GetDescription(), "GTiff"))
4587 : {
4588 : VSILFILE *fpTemp =
4589 16 : VSIFOpenL(m_poExternalDS->GetDescription(), "rb");
4590 16 : if (fpTemp)
4591 : {
4592 16 : GByte abySignature[4] = {0};
4593 16 : VSIFReadL(abySignature, 1, 4, fpTemp);
4594 16 : VSIFCloseL(fpTemp);
4595 16 : const bool bBigTIFF =
4596 16 : abySignature[2] == 43 || abySignature[3] == 43;
4597 : m_osHeaderParsingStandard =
4598 16 : bBigTIFF ? BIGTIFF_GEOTIFF_STRING : TIFF_GEOTIFF_STRING;
4599 : const char *pszOffset =
4600 16 : m_poExternalDS->GetRasterBand(1)->GetMetadataItem(
4601 16 : "BLOCK_OFFSET_0_0", "TIFF");
4602 16 : if (pszOffset)
4603 16 : m_nBaseOffset = CPLAtoGIntBig(pszOffset);
4604 : }
4605 : }
4606 :
4607 131 : if (!m_osHeaderParsingStandard.empty() && m_nBaseOffset > 0)
4608 : {
4609 22 : CPLXMLNode *psHeader = CPLCreateXMLNode(
4610 44 : psFAO, CXT_Element, (osPrefix + "Header").c_str());
4611 22 : CPLAddXMLAttributeAndValue(
4612 : CPLCreateXMLElementAndValue(
4613 44 : psHeader, (osPrefix + "offset").c_str(), "0"),
4614 : "unit", "byte");
4615 22 : CPLAddXMLAttributeAndValue(
4616 : CPLCreateXMLElementAndValue(
4617 44 : psHeader, (osPrefix + "object_length").c_str(),
4618 : CPLSPrintf(CPL_FRMT_GUIB,
4619 22 : static_cast<GUIntBig>(m_nBaseOffset))),
4620 : "unit", "byte");
4621 44 : CPLCreateXMLElementAndValue(
4622 44 : psHeader, (osPrefix + "parsing_standard_id").c_str(),
4623 : m_osHeaderParsingStandard.c_str());
4624 22 : if (m_osHeaderParsingStandard == TIFF_GEOTIFF_STRING)
4625 : {
4626 18 : CPLCreateXMLElementAndValue(
4627 36 : psHeader, (osPrefix + "description").c_str(),
4628 : "TIFF/GeoTIFF header. The TIFF/GeoTIFF format is used "
4629 : "throughout the geospatial and science communities "
4630 : "to share geographic image data. ");
4631 : }
4632 4 : else if (m_osHeaderParsingStandard == BIGTIFF_GEOTIFF_STRING)
4633 : {
4634 0 : CPLCreateXMLElementAndValue(
4635 0 : psHeader, (osPrefix + "description").c_str(),
4636 : "BigTIFF/GeoTIFF header. The BigTIFF/GeoTIFF format is "
4637 : "used "
4638 : "throughout the geospatial and science communities "
4639 : "to share geographic image data. ");
4640 : }
4641 : }
4642 :
4643 131 : WriteArray(osPrefix, psFAO, pszLocalIdentifier,
4644 : psTemplateSpecialConstants);
4645 : }
4646 : }
4647 : }
4648 :
4649 : /************************************************************************/
4650 : /* WriteHeader() */
4651 : /************************************************************************/
4652 :
4653 193 : void PDS4Dataset::WriteHeader()
4654 : {
4655 : const bool bAppend =
4656 193 : CPLFetchBool(m_papszCreationOptions, "APPEND_SUBDATASET", false);
4657 193 : if (bAppend)
4658 : {
4659 4 : WriteHeaderAppendCase();
4660 5 : return;
4661 : }
4662 :
4663 : CPLXMLNode *psRoot;
4664 189 : if (m_bCreateHeader)
4665 : {
4666 : CPLString osTemplateFilename =
4667 182 : CSLFetchNameValueDef(m_papszCreationOptions, "TEMPLATE", "");
4668 182 : if (!osTemplateFilename.empty())
4669 : {
4670 20 : if (STARTS_WITH(osTemplateFilename, "http://") ||
4671 10 : STARTS_WITH(osTemplateFilename, "https://"))
4672 : {
4673 0 : osTemplateFilename = "/vsicurl_streaming/" + osTemplateFilename;
4674 : }
4675 10 : psRoot = CPLParseXMLFile(osTemplateFilename);
4676 : }
4677 172 : else if (!m_osXMLPDS4.empty())
4678 6 : psRoot = CPLParseXMLString(m_osXMLPDS4);
4679 : else
4680 : {
4681 : #ifndef USE_ONLY_EMBEDDED_RESOURCE_FILES
4682 : #ifdef EMBED_RESOURCE_FILES
4683 : CPLErrorStateBackuper oErrorStateBackuper(CPLQuietErrorHandler);
4684 : #endif
4685 : const char *pszDefaultTemplateFilename =
4686 166 : CPLFindFile("gdal", "pds4_template.xml");
4687 166 : if (pszDefaultTemplateFilename)
4688 : {
4689 166 : psRoot = CPLParseXMLFile(pszDefaultTemplateFilename);
4690 : }
4691 : else
4692 : #endif
4693 : {
4694 : #ifdef EMBED_RESOURCE_FILES
4695 : static const bool bOnce [[maybe_unused]] = []()
4696 : {
4697 : CPLDebug("PDS4", "Using embedded pds4_template.xml");
4698 : return true;
4699 : }();
4700 : psRoot = CPLParseXMLString(PDS4GetEmbeddedTemplate());
4701 : #else
4702 0 : CPLError(CE_Failure, CPLE_AppDefined,
4703 : "Cannot find pds4_template.xml and TEMPLATE "
4704 : "creation option not specified");
4705 0 : return;
4706 : #endif
4707 : }
4708 : }
4709 : }
4710 : else
4711 : {
4712 7 : psRoot = CPLParseXMLFile(m_osXMLFilename);
4713 : }
4714 189 : CPLXMLTreeCloser oCloser(psRoot);
4715 189 : psRoot = oCloser.get();
4716 189 : if (psRoot == nullptr)
4717 0 : return;
4718 189 : CPLXMLNode *psProduct = CPLGetXMLNode(psRoot, "=Product_Observational");
4719 189 : if (psProduct == nullptr)
4720 : {
4721 1 : psProduct = CPLGetXMLNode(psRoot, "=pds:Product_Observational");
4722 : }
4723 189 : if (psProduct == nullptr)
4724 : {
4725 1 : CPLError(CE_Failure, CPLE_AppDefined,
4726 : "Cannot find Product_Observational element in template");
4727 1 : return;
4728 : }
4729 :
4730 188 : if (m_bCreateHeader)
4731 : {
4732 362 : CPLString osCARTVersion(CURRENT_CART_VERSION);
4733 181 : char *pszXML = CPLSerializeXMLTree(psRoot);
4734 181 : if (pszXML)
4735 : {
4736 181 : const char *pszIter = pszXML;
4737 : while (true)
4738 : {
4739 355 : const char *pszCartSchema = strstr(pszIter, "PDS4_CART_");
4740 355 : if (pszCartSchema)
4741 : {
4742 348 : const char *pszXSDExtension = strstr(pszCartSchema, ".xsd");
4743 348 : if (pszXSDExtension &&
4744 348 : pszXSDExtension - pszCartSchema <= 20)
4745 : {
4746 174 : osCARTVersion = pszCartSchema + strlen("PDS4_CART_");
4747 174 : osCARTVersion.resize(pszXSDExtension - pszCartSchema -
4748 : strlen("PDS4_CART_"));
4749 174 : break;
4750 : }
4751 : else
4752 : {
4753 174 : pszIter = pszCartSchema + 1;
4754 : }
4755 : }
4756 : else
4757 : {
4758 7 : break;
4759 : }
4760 174 : }
4761 :
4762 181 : CPLFree(pszXML);
4763 : }
4764 :
4765 181 : CreateHeader(psProduct, osCARTVersion.c_str());
4766 : }
4767 :
4768 188 : WriteVectorLayers(psProduct);
4769 :
4770 188 : CPLSerializeXMLTreeToFile(psRoot, GetDescription());
4771 : }
4772 :
4773 : /************************************************************************/
4774 : /* ICreateLayer() */
4775 : /************************************************************************/
4776 :
4777 66 : OGRLayer *PDS4Dataset::ICreateLayer(const char *pszName,
4778 : const OGRGeomFieldDefn *poGeomFieldDefn,
4779 : CSLConstList papszOptions)
4780 : {
4781 : const char *pszTableType =
4782 66 : CSLFetchNameValueDef(papszOptions, "TABLE_TYPE", "DELIMITED");
4783 66 : if (!EQUAL(pszTableType, "CHARACTER") && !EQUAL(pszTableType, "BINARY") &&
4784 55 : !EQUAL(pszTableType, "DELIMITED"))
4785 : {
4786 0 : return nullptr;
4787 : }
4788 :
4789 66 : const auto eGType = poGeomFieldDefn ? poGeomFieldDefn->GetType() : wkbNone;
4790 : const auto poSpatialRef =
4791 66 : poGeomFieldDefn ? poGeomFieldDefn->GetSpatialRef() : nullptr;
4792 :
4793 126 : const char *pszExt = EQUAL(pszTableType, "CHARACTER") ? "dat"
4794 60 : : EQUAL(pszTableType, "BINARY") ? "bin"
4795 : : "csv";
4796 132 : std::string osBasename(pszName);
4797 629 : for (char &ch : osBasename)
4798 : {
4799 563 : if (!isalnum(static_cast<unsigned char>(ch)) &&
4800 53 : static_cast<unsigned>(ch) <= 127)
4801 53 : ch = '_';
4802 : }
4803 :
4804 198 : CPLString osFullFilename(CPLFormFilenameSafe(
4805 132 : CPLGetPathSafe(m_osXMLFilename.c_str()).c_str(),
4806 198 : CPLGetBasenameSafe(m_osXMLFilename.c_str()).c_str(), nullptr));
4807 66 : osFullFilename += '_';
4808 66 : osFullFilename += osBasename.c_str();
4809 66 : osFullFilename += '.';
4810 66 : osFullFilename += pszExt;
4811 : VSIStatBufL sStat;
4812 66 : if (VSIStatL(osFullFilename, &sStat) == 0)
4813 : {
4814 0 : CPLError(CE_Failure, CPLE_AppDefined,
4815 : "%s already exists. Please delete it before, or "
4816 : "rename the layer",
4817 : osFullFilename.c_str());
4818 0 : return nullptr;
4819 : }
4820 :
4821 66 : if (EQUAL(pszTableType, "DELIMITED"))
4822 : {
4823 : auto poLayer = std::make_unique<PDS4DelimitedTable>(
4824 55 : this, pszName, osFullFilename, true);
4825 55 : if (!poLayer->InitializeNewLayer(poSpatialRef, false, eGType,
4826 : papszOptions))
4827 : {
4828 1 : return nullptr;
4829 : }
4830 54 : m_apoLayers.push_back(
4831 108 : std::make_unique<PDS4EditableLayer>(std::move(poLayer)));
4832 : }
4833 : else
4834 : {
4835 0 : std::unique_ptr<PDS4FixedWidthTable> poLayer;
4836 11 : if (EQUAL(pszTableType, "CHARACTER"))
4837 6 : poLayer = std::make_unique<PDS4TableCharacter>(
4838 12 : this, pszName, osFullFilename, true);
4839 : else
4840 5 : poLayer = std::make_unique<PDS4TableBinary>(this, pszName,
4841 10 : osFullFilename, true);
4842 11 : if (!poLayer->InitializeNewLayer(poSpatialRef, false, eGType,
4843 : papszOptions))
4844 : {
4845 0 : return nullptr;
4846 : }
4847 11 : m_apoLayers.push_back(
4848 22 : std::make_unique<PDS4EditableLayer>(std::move(poLayer)));
4849 : }
4850 65 : return dynamic_cast<OGRLayer *>(m_apoLayers.back().get());
4851 : }
4852 :
4853 : /************************************************************************/
4854 : /* TestCapability() */
4855 : /************************************************************************/
4856 :
4857 69 : bool PDS4Dataset::TestCapability(const char *pszCap) const
4858 : {
4859 69 : if (EQUAL(pszCap, ODsCCreateLayer))
4860 35 : return eAccess == GA_Update;
4861 34 : else if (EQUAL(pszCap, ODsCZGeometries))
4862 6 : return TRUE;
4863 : else
4864 28 : return FALSE;
4865 : }
4866 :
4867 : /************************************************************************/
4868 : /* Create() */
4869 : /************************************************************************/
4870 :
4871 140 : GDALDataset *PDS4Dataset::Create(const char *pszFilename, int nXSize,
4872 : int nYSize, int nBandsIn, GDALDataType eType,
4873 : CSLConstList papszOptions)
4874 : {
4875 280 : return CreateInternal(pszFilename, nullptr, nXSize, nYSize, nBandsIn, eType,
4876 : papszOptions)
4877 140 : .release();
4878 : }
4879 :
4880 : /************************************************************************/
4881 : /* CreateInternal() */
4882 : /************************************************************************/
4883 :
4884 200 : std::unique_ptr<PDS4Dataset> PDS4Dataset::CreateInternal(
4885 : const char *pszFilename, GDALDataset *poSrcDS, int nXSize, int nYSize,
4886 : int nBandsIn, GDALDataType eType, const char *const *papszOptionsIn)
4887 : {
4888 400 : CPLStringList aosOptions(papszOptionsIn);
4889 :
4890 200 : if (nXSize == 0 && nYSize == 0 && nBandsIn == 0 && eType == GDT_Unknown)
4891 : {
4892 : // Vector file creation
4893 94 : auto poDS = std::make_unique<PDS4Dataset>();
4894 47 : poDS->SetDescription(pszFilename);
4895 47 : poDS->nRasterXSize = 0;
4896 47 : poDS->nRasterYSize = 0;
4897 47 : poDS->eAccess = GA_Update;
4898 47 : poDS->m_osXMLFilename = pszFilename;
4899 47 : poDS->m_bCreateHeader = true;
4900 47 : poDS->m_bStripFileAreaObservationalFromTemplate = true;
4901 47 : poDS->m_papszCreationOptions = CSLDuplicate(aosOptions.List());
4902 47 : poDS->m_bUseSrcLabel = aosOptions.FetchBool("USE_SRC_LABEL", true);
4903 47 : return poDS;
4904 : }
4905 :
4906 153 : if (nXSize == 0)
4907 0 : return nullptr;
4908 :
4909 153 : if (!(eType == GDT_UInt8 || eType == GDT_Int8 || eType == GDT_Int16 ||
4910 50 : eType == GDT_UInt16 || eType == GDT_Int32 || eType == GDT_UInt32 ||
4911 35 : eType == GDT_Int64 || eType == GDT_UInt64 || eType == GDT_Float32 ||
4912 20 : eType == GDT_Float64 || eType == GDT_CFloat32 ||
4913 10 : eType == GDT_CFloat64))
4914 : {
4915 6 : CPLError(
4916 : CE_Failure, CPLE_NotSupported,
4917 : "The PDS4 driver does not supporting creating files of type %s.",
4918 : GDALGetDataTypeName(eType));
4919 6 : return nullptr;
4920 : }
4921 :
4922 147 : if (nBandsIn == 0)
4923 : {
4924 1 : CPLError(CE_Failure, CPLE_NotSupported, "Invalid number of bands");
4925 1 : return nullptr;
4926 : }
4927 :
4928 : const char *pszArrayType =
4929 146 : aosOptions.FetchNameValueDef("ARRAY_TYPE", "Array_3D_Image");
4930 146 : const bool bIsArray2D = STARTS_WITH(pszArrayType, "Array_2D");
4931 146 : if (nBandsIn > 1 && bIsArray2D)
4932 : {
4933 1 : CPLError(CE_Failure, CPLE_NotSupported,
4934 : "ARRAY_TYPE=%s is not supported for a multi-band raster",
4935 : pszArrayType);
4936 1 : return nullptr;
4937 : }
4938 :
4939 : /* -------------------------------------------------------------------- */
4940 : /* Compute pixel, line and band offsets */
4941 : /* -------------------------------------------------------------------- */
4942 145 : const int nItemSize = GDALGetDataTypeSizeBytes(eType);
4943 : int nLineOffset, nPixelOffset;
4944 : vsi_l_offset nBandOffset;
4945 :
4946 : const char *pszInterleave =
4947 145 : aosOptions.FetchNameValueDef(GDALMD_INTERLEAVE, "BSQ");
4948 145 : if (bIsArray2D)
4949 1 : pszInterleave = "BIP";
4950 :
4951 145 : if (EQUAL(pszInterleave, "BIP"))
4952 : {
4953 4 : nPixelOffset = nItemSize * nBandsIn;
4954 4 : if (nPixelOffset > INT_MAX / nBandsIn)
4955 : {
4956 0 : return nullptr;
4957 : }
4958 4 : nLineOffset = nPixelOffset * nXSize;
4959 4 : nBandOffset = nItemSize;
4960 : }
4961 141 : else if (EQUAL(pszInterleave, "BSQ"))
4962 : {
4963 137 : nPixelOffset = nItemSize;
4964 137 : if (nPixelOffset > INT_MAX / nXSize)
4965 : {
4966 0 : return nullptr;
4967 : }
4968 137 : nLineOffset = nPixelOffset * nXSize;
4969 137 : nBandOffset = static_cast<vsi_l_offset>(nLineOffset) * nYSize;
4970 : }
4971 4 : else if (EQUAL(pszInterleave, "BIL"))
4972 : {
4973 3 : nPixelOffset = nItemSize;
4974 3 : if (nPixelOffset > INT_MAX / nBandsIn ||
4975 3 : nPixelOffset * nBandsIn > INT_MAX / nXSize)
4976 : {
4977 0 : return nullptr;
4978 : }
4979 3 : nLineOffset = nItemSize * nBandsIn * nXSize;
4980 3 : nBandOffset = static_cast<vsi_l_offset>(nItemSize) * nXSize;
4981 : }
4982 : else
4983 : {
4984 1 : CPLError(CE_Failure, CPLE_NotSupported, "Invalid value for INTERLEAVE");
4985 1 : return nullptr;
4986 : }
4987 :
4988 : const char *pszImageFormat =
4989 144 : aosOptions.FetchNameValueDef("IMAGE_FORMAT", "RAW");
4990 144 : const char *pszImageExtension = aosOptions.FetchNameValueDef(
4991 144 : "IMAGE_EXTENSION", EQUAL(pszImageFormat, "RAW") ? "img" : "tif");
4992 : CPLString osImageFilename(aosOptions.FetchNameValueDef(
4993 : "IMAGE_FILENAME",
4994 288 : CPLResetExtensionSafe(pszFilename, pszImageExtension).c_str()));
4995 :
4996 144 : const bool bAppend = aosOptions.FetchBool("APPEND_SUBDATASET", false);
4997 144 : if (bAppend)
4998 : {
4999 4 : GDALOpenInfo oOpenInfo(pszFilename, GA_ReadOnly);
5000 4 : auto poExistingPDS4 = OpenInternal(&oOpenInfo);
5001 4 : if (!poExistingPDS4)
5002 : {
5003 0 : return nullptr;
5004 : }
5005 4 : osImageFilename = poExistingPDS4->m_osImageFilename;
5006 4 : poExistingPDS4.reset();
5007 :
5008 : auto poImageDS = std::unique_ptr<GDALDataset>(
5009 8 : GDALDataset::Open(osImageFilename, GDAL_OF_RASTER));
5010 6 : if (poImageDS && poImageDS->GetDriver() &&
5011 2 : EQUAL(poImageDS->GetDriver()->GetDescription(), "GTiff"))
5012 : {
5013 2 : pszImageFormat = "GEOTIFF";
5014 : }
5015 : }
5016 :
5017 144 : GDALDataset *poExternalDS = nullptr;
5018 144 : VSILFILE *fpImage = nullptr;
5019 144 : vsi_l_offset nBaseOffset = 0;
5020 144 : bool bIsLSB = true;
5021 288 : CPLString osHeaderParsingStandard;
5022 : const bool bCreateLabelOnly =
5023 144 : aosOptions.FetchBool("CREATE_LABEL_ONLY", false);
5024 144 : if (bCreateLabelOnly)
5025 : {
5026 8 : if (poSrcDS == nullptr)
5027 : {
5028 0 : CPLError(
5029 : CE_Failure, CPLE_AppDefined,
5030 : "CREATE_LABEL_ONLY is only compatible with CreateCopy() mode");
5031 1 : return nullptr;
5032 : }
5033 8 : RawBinaryLayout sLayout;
5034 8 : if (!poSrcDS->GetRawBinaryLayout(sLayout))
5035 : {
5036 1 : CPLError(CE_Failure, CPLE_AppDefined,
5037 : "Source dataset is not compatible with raw binary format");
5038 1 : return nullptr;
5039 : }
5040 7 : if ((nBandsIn > 1 &&
5041 7 : sLayout.eInterleaving == RawBinaryLayout::Interleaving::UNKNOWN) ||
5042 6 : (nBandsIn == 1 &&
5043 6 : !(sLayout.nPixelOffset == nItemSize &&
5044 6 : sLayout.nLineOffset == sLayout.nPixelOffset * nXSize)))
5045 : {
5046 0 : CPLError(CE_Failure, CPLE_AppDefined,
5047 : "Source dataset has an interleaving not handled in PDS4");
5048 0 : return nullptr;
5049 : }
5050 7 : fpImage = VSIFOpenL(sLayout.osRawFilename.c_str(), "rb");
5051 7 : if (fpImage == nullptr)
5052 : {
5053 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot open raw image %s",
5054 : sLayout.osRawFilename.c_str());
5055 0 : return nullptr;
5056 : }
5057 7 : osImageFilename = sLayout.osRawFilename;
5058 7 : if (nBandsIn == 1 ||
5059 1 : sLayout.eInterleaving == RawBinaryLayout::Interleaving::BIP)
5060 7 : pszInterleave = "BIP";
5061 0 : else if (sLayout.eInterleaving == RawBinaryLayout::Interleaving::BIL)
5062 0 : pszInterleave = "BIL";
5063 : else
5064 0 : pszInterleave = "BSQ";
5065 7 : nBaseOffset = sLayout.nImageOffset;
5066 7 : nPixelOffset = static_cast<int>(sLayout.nPixelOffset);
5067 7 : nLineOffset = static_cast<int>(sLayout.nLineOffset);
5068 7 : nBandOffset = static_cast<vsi_l_offset>(sLayout.nBandOffset);
5069 7 : bIsLSB = sLayout.bLittleEndianOrder;
5070 7 : auto poSrcDriver = poSrcDS->GetDriver();
5071 7 : if (poSrcDriver)
5072 : {
5073 7 : auto pszDriverName = poSrcDriver->GetDescription();
5074 7 : if (EQUAL(pszDriverName, "GTiff"))
5075 : {
5076 2 : GByte abySignature[4] = {0};
5077 2 : VSIFReadL(abySignature, 1, 4, fpImage);
5078 2 : const bool bBigTIFF =
5079 2 : abySignature[2] == 43 || abySignature[3] == 43;
5080 : osHeaderParsingStandard =
5081 2 : bBigTIFF ? BIGTIFF_GEOTIFF_STRING : TIFF_GEOTIFF_STRING;
5082 : }
5083 5 : else if (EQUAL(pszDriverName, "ISIS3"))
5084 : {
5085 1 : osHeaderParsingStandard = "ISIS3";
5086 : }
5087 4 : else if (EQUAL(pszDriverName, "VICAR"))
5088 : {
5089 1 : osHeaderParsingStandard = "VICAR2";
5090 : }
5091 3 : else if (EQUAL(pszDriverName, "PDS"))
5092 : {
5093 1 : osHeaderParsingStandard = "PDS3";
5094 : }
5095 2 : else if (EQUAL(pszDriverName, "FITS"))
5096 : {
5097 1 : osHeaderParsingStandard = "FITS 3.0";
5098 : aosOptions.SetNameValue("VAR_VERTICAL_DISPLAY_DIRECTION",
5099 1 : "Bottom to Top");
5100 : }
5101 : }
5102 : }
5103 136 : else if (EQUAL(pszImageFormat, "GEOTIFF"))
5104 : {
5105 20 : if (EQUAL(pszInterleave, "BIL"))
5106 : {
5107 2 : if (aosOptions.FetchBool("@INTERLEAVE_ADDED_AUTOMATICALLY", false))
5108 : {
5109 1 : pszInterleave = "BSQ";
5110 : }
5111 : else
5112 : {
5113 1 : CPLError(CE_Failure, CPLE_AppDefined,
5114 : "INTERLEAVE=BIL not supported for GeoTIFF in PDS4");
5115 1 : return nullptr;
5116 : }
5117 : }
5118 : GDALDriver *poDrv =
5119 19 : static_cast<GDALDriver *>(GDALGetDriverByName("GTiff"));
5120 19 : if (poDrv == nullptr)
5121 : {
5122 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find GTiff driver");
5123 0 : return nullptr;
5124 : }
5125 19 : char **papszGTiffOptions = nullptr;
5126 : #ifdef notdef
5127 : // In practice I can't see which option we can really use
5128 : const char *pszGTiffOptions =
5129 : CSLFetchNameValueDef(papszOptions, "GEOTIFF_OPTIONS", "");
5130 : char **papszTokens = CSLTokenizeString2(pszGTiffOptions, ",", 0);
5131 : if (CPLFetchBool(papszTokens, "TILED", false))
5132 : {
5133 : CSLDestroy(papszTokens);
5134 : CPLError(CE_Failure, CPLE_AppDefined,
5135 : "Tiled GeoTIFF is not supported for PDS4");
5136 : return NULL;
5137 : }
5138 : if (!EQUAL(CSLFetchNameValueDef(papszTokens, "COMPRESS", "NONE"),
5139 : "NONE"))
5140 : {
5141 : CSLDestroy(papszTokens);
5142 : CPLError(CE_Failure, CPLE_AppDefined,
5143 : "Compressed GeoTIFF is not supported for PDS4");
5144 : return NULL;
5145 : }
5146 : papszGTiffOptions =
5147 : CSLSetNameValue(papszGTiffOptions, "ENDIANNESS", "LITTLE");
5148 : for (int i = 0; papszTokens[i] != NULL; i++)
5149 : {
5150 : papszGTiffOptions = CSLAddString(papszGTiffOptions, papszTokens[i]);
5151 : }
5152 : CSLDestroy(papszTokens);
5153 : #endif
5154 :
5155 : papszGTiffOptions =
5156 19 : CSLSetNameValue(papszGTiffOptions, GDALMD_INTERLEAVE,
5157 19 : EQUAL(pszInterleave, "BSQ") ? "BAND" : "PIXEL");
5158 : // Will make sure that our blocks at nodata are not optimized
5159 : // away but indeed well written
5160 19 : papszGTiffOptions = CSLSetNameValue(
5161 : papszGTiffOptions, "@WRITE_EMPTY_TILES_SYNCHRONOUSLY", "YES");
5162 19 : if (nBandsIn > 1 && EQUAL(pszInterleave, "BSQ"))
5163 : {
5164 : papszGTiffOptions =
5165 2 : CSLSetNameValue(papszGTiffOptions, "BLOCKYSIZE", "1");
5166 : }
5167 :
5168 19 : if (bAppend)
5169 : {
5170 : papszGTiffOptions =
5171 2 : CSLAddString(papszGTiffOptions, "APPEND_SUBDATASET=YES");
5172 : }
5173 :
5174 19 : poExternalDS = poDrv->Create(osImageFilename, nXSize, nYSize, nBandsIn,
5175 : eType, papszGTiffOptions);
5176 19 : CSLDestroy(papszGTiffOptions);
5177 19 : if (poExternalDS == nullptr)
5178 : {
5179 1 : CPLError(CE_Failure, CPLE_FileIO, "Cannot create %s",
5180 : osImageFilename.c_str());
5181 1 : return nullptr;
5182 : }
5183 : }
5184 : else
5185 : {
5186 230 : fpImage = VSIFOpenL(
5187 : osImageFilename,
5188 : bAppend ? "rb+"
5189 114 : : VSISupportsRandomWrite(osImageFilename.c_str(), true) ? "wb+"
5190 : : "wb");
5191 116 : if (fpImage == nullptr)
5192 : {
5193 2 : CPLError(CE_Failure, CPLE_FileIO, "Cannot create %s",
5194 : osImageFilename.c_str());
5195 2 : return nullptr;
5196 : }
5197 114 : if (bAppend)
5198 : {
5199 2 : VSIFSeekL(fpImage, 0, SEEK_END);
5200 2 : nBaseOffset = VSIFTellL(fpImage);
5201 : }
5202 : }
5203 :
5204 278 : auto poDS = std::make_unique<PDS4Dataset>();
5205 139 : poDS->SetDescription(pszFilename);
5206 139 : poDS->m_bMustInitImageFile = true;
5207 139 : poDS->m_fpImage = fpImage;
5208 139 : poDS->m_nBaseOffset = nBaseOffset;
5209 139 : poDS->m_poExternalDS = poExternalDS;
5210 139 : poDS->nRasterXSize = nXSize;
5211 139 : poDS->nRasterYSize = nYSize;
5212 139 : poDS->eAccess = GA_Update;
5213 139 : poDS->m_osImageFilename = std::move(osImageFilename);
5214 139 : poDS->m_bCreateHeader = true;
5215 139 : poDS->m_bStripFileAreaObservationalFromTemplate = true;
5216 139 : poDS->m_osInterleave = pszInterleave;
5217 139 : poDS->m_papszCreationOptions = CSLDuplicate(aosOptions.List());
5218 139 : poDS->m_bUseSrcLabel = aosOptions.FetchBool("USE_SRC_LABEL", true);
5219 139 : poDS->m_bIsLSB = bIsLSB;
5220 139 : poDS->m_osHeaderParsingStandard = std::move(osHeaderParsingStandard);
5221 139 : poDS->m_bCreatedFromExistingBinaryFile = bCreateLabelOnly;
5222 :
5223 139 : if (EQUAL(pszInterleave, "BIP"))
5224 : {
5225 10 : poDS->GDALDataset::SetMetadataItem(GDALMD_INTERLEAVE, "PIXEL",
5226 : GDAL_MDD_IMAGE_STRUCTURE);
5227 : }
5228 129 : else if (EQUAL(pszInterleave, "BSQ"))
5229 : {
5230 128 : poDS->GDALDataset::SetMetadataItem(GDALMD_INTERLEAVE, "BAND",
5231 : GDAL_MDD_IMAGE_STRUCTURE);
5232 : }
5233 :
5234 338 : for (int i = 0; i < nBandsIn; i++)
5235 : {
5236 199 : if (poDS->m_poExternalDS != nullptr)
5237 : {
5238 : auto poBand = std::make_unique<PDS4WrapperRasterBand>(
5239 24 : poDS->m_poExternalDS->GetRasterBand(i + 1));
5240 24 : poDS->SetBand(i + 1, std::move(poBand));
5241 : }
5242 : else
5243 : {
5244 : auto poBand = std::make_unique<PDS4RawRasterBand>(
5245 175 : poDS.get(), i + 1, poDS->m_fpImage,
5246 175 : poDS->m_nBaseOffset + nBandOffset * i, nPixelOffset,
5247 : nLineOffset, eType,
5248 175 : bIsLSB ? RawRasterBand::ByteOrder::ORDER_LITTLE_ENDIAN
5249 175 : : RawRasterBand::ByteOrder::ORDER_BIG_ENDIAN);
5250 175 : poDS->SetBand(i + 1, std::move(poBand));
5251 : }
5252 : }
5253 :
5254 139 : return poDS;
5255 : }
5256 :
5257 : /************************************************************************/
5258 : /* PDS4GetUnderlyingDataset() */
5259 : /************************************************************************/
5260 :
5261 64 : static GDALDataset *PDS4GetUnderlyingDataset(GDALDataset *poSrcDS)
5262 : {
5263 128 : if (poSrcDS->GetDriver() != nullptr &&
5264 64 : poSrcDS->GetDriver() == GDALGetDriverByName("VRT"))
5265 : {
5266 3 : VRTDataset *poVRTDS = cpl::down_cast<VRTDataset *>(poSrcDS);
5267 3 : poSrcDS = poVRTDS->GetSingleSimpleSource();
5268 : }
5269 :
5270 64 : return poSrcDS;
5271 : }
5272 :
5273 : /************************************************************************/
5274 : /* CreateCopy() */
5275 : /************************************************************************/
5276 :
5277 64 : GDALDataset *PDS4Dataset::CreateCopy(const char *pszFilename,
5278 : GDALDataset *poSrcDS, int bStrict,
5279 : CSLConstList papszOptions,
5280 : GDALProgressFunc pfnProgress,
5281 : void *pProgressData)
5282 : {
5283 : const char *pszImageFormat =
5284 64 : CSLFetchNameValueDef(papszOptions, "IMAGE_FORMAT", "RAW");
5285 64 : GDALDataset *poSrcUnderlyingDS = PDS4GetUnderlyingDataset(poSrcDS);
5286 64 : if (poSrcUnderlyingDS == nullptr)
5287 0 : poSrcUnderlyingDS = poSrcDS;
5288 72 : if (EQUAL(pszImageFormat, "GEOTIFF") &&
5289 8 : strcmp(poSrcUnderlyingDS->GetDescription(),
5290 : CSLFetchNameValueDef(
5291 : papszOptions, "IMAGE_FILENAME",
5292 72 : CPLResetExtensionSafe(pszFilename, "tif").c_str())) == 0)
5293 : {
5294 1 : CPLError(CE_Failure, CPLE_NotSupported,
5295 : "Output file has same name as input file");
5296 1 : return nullptr;
5297 : }
5298 63 : if (poSrcDS->GetRasterCount() == 0)
5299 : {
5300 1 : CPLError(CE_Failure, CPLE_NotSupported, "Unsupported band count");
5301 1 : return nullptr;
5302 : }
5303 :
5304 62 : const bool bAppend = CPLFetchBool(papszOptions, "APPEND_SUBDATASET", false);
5305 62 : if (bAppend)
5306 : {
5307 6 : GDALOpenInfo oOpenInfo(pszFilename, GA_ReadOnly);
5308 6 : auto poExistingDS = OpenInternal(&oOpenInfo);
5309 6 : if (poExistingDS)
5310 : {
5311 6 : GDALGeoTransform existingGT;
5312 : const bool bExistingHasGT =
5313 6 : poExistingDS->GetGeoTransform(existingGT) == CE_None;
5314 6 : GDALGeoTransform gt;
5315 6 : const bool bSrcHasGT = poSrcDS->GetGeoTransform(gt) == CE_None;
5316 :
5317 6 : CPLString osExistingProj4;
5318 6 : if (const auto poExistingSRS = poExistingDS->GetSpatialRef())
5319 : {
5320 6 : char *pszExistingProj4 = nullptr;
5321 6 : poExistingSRS->exportToProj4(&pszExistingProj4);
5322 6 : if (pszExistingProj4)
5323 6 : osExistingProj4 = pszExistingProj4;
5324 6 : CPLFree(pszExistingProj4);
5325 : }
5326 6 : CPLString osSrcProj4;
5327 6 : if (const auto poSrcSRS = poSrcDS->GetSpatialRef())
5328 : {
5329 6 : char *pszSrcProj4 = nullptr;
5330 6 : poSrcSRS->exportToProj4(&pszSrcProj4);
5331 6 : if (pszSrcProj4)
5332 6 : osSrcProj4 = pszSrcProj4;
5333 6 : CPLFree(pszSrcProj4);
5334 : }
5335 :
5336 6 : poExistingDS.reset();
5337 :
5338 : const auto maxRelErrorGT =
5339 6 : [](const GDALGeoTransform >1, const GDALGeoTransform >2)
5340 : {
5341 6 : double maxRelError = 0.0;
5342 42 : for (int i = 0; i < 6; i++)
5343 : {
5344 36 : if (gt1[i] == 0.0)
5345 : {
5346 12 : maxRelError = std::max(maxRelError, std::abs(gt2[i]));
5347 : }
5348 : else
5349 : {
5350 24 : maxRelError =
5351 48 : std::max(maxRelError, std::abs(gt2[i] - gt1[i]) /
5352 24 : std::abs(gt1[i]));
5353 : }
5354 : }
5355 6 : return maxRelError;
5356 : };
5357 :
5358 6 : if ((bExistingHasGT && !bSrcHasGT) ||
5359 18 : (!bExistingHasGT && bSrcHasGT) ||
5360 6 : (bExistingHasGT && bSrcHasGT &&
5361 6 : maxRelErrorGT(existingGT, gt) > 1e-10))
5362 : {
5363 1 : CPLError(bStrict ? CE_Failure : CE_Warning, CPLE_NotSupported,
5364 : "Appending to a dataset with a different "
5365 : "geotransform is not supported");
5366 1 : if (bStrict)
5367 1 : return nullptr;
5368 : }
5369 : // Do proj string comparison, as it is unlikely that
5370 : // OGRSpatialReference::IsSame() will lead to identical reasons due
5371 : // to PDS changing CRS names, etc...
5372 5 : if (osExistingProj4 != osSrcProj4)
5373 : {
5374 1 : CPLError(bStrict ? CE_Failure : CE_Warning, CPLE_NotSupported,
5375 : "Appending to a dataset with a different "
5376 : "coordinate reference system is not supported");
5377 1 : if (bStrict)
5378 1 : return nullptr;
5379 : }
5380 : }
5381 : }
5382 :
5383 60 : const int nXSize = poSrcDS->GetRasterXSize();
5384 60 : const int nYSize = poSrcDS->GetRasterYSize();
5385 60 : const int nBands = poSrcDS->GetRasterCount();
5386 60 : GDALDataType eType = poSrcDS->GetRasterBand(1)->GetRasterDataType();
5387 : auto poDS = CreateInternal(pszFilename, poSrcDS, nXSize, nYSize, nBands,
5388 120 : eType, papszOptions);
5389 60 : if (poDS == nullptr)
5390 5 : return nullptr;
5391 :
5392 55 : GDALGeoTransform gt;
5393 55 : if (poSrcDS->GetGeoTransform(gt) == CE_None && gt != GDALGeoTransform())
5394 : {
5395 52 : poDS->SetGeoTransform(gt);
5396 : }
5397 :
5398 110 : if (poSrcDS->GetProjectionRef() != nullptr &&
5399 55 : strlen(poSrcDS->GetProjectionRef()) > 0)
5400 : {
5401 52 : poDS->SetProjection(poSrcDS->GetProjectionRef());
5402 : }
5403 :
5404 137 : for (int i = 1; i <= nBands; i++)
5405 : {
5406 82 : int bHasNoData = false;
5407 :
5408 82 : if (poSrcDS->GetRasterBand(i)->GetRasterDataType() == GDT_Int64)
5409 : {
5410 : const auto nNoData =
5411 0 : poSrcDS->GetRasterBand(i)->GetNoDataValueAsInt64(&bHasNoData);
5412 0 : if (bHasNoData)
5413 0 : poDS->GetRasterBand(i)->SetNoDataValueAsInt64(nNoData);
5414 : }
5415 82 : else if (poSrcDS->GetRasterBand(i)->GetRasterDataType() == GDT_UInt64)
5416 : {
5417 : const auto nNoData =
5418 0 : poSrcDS->GetRasterBand(i)->GetNoDataValueAsUInt64(&bHasNoData);
5419 0 : if (bHasNoData)
5420 0 : poDS->GetRasterBand(i)->SetNoDataValueAsUInt64(nNoData);
5421 : }
5422 : else
5423 : {
5424 : const double dfNoData =
5425 82 : poSrcDS->GetRasterBand(i)->GetNoDataValue(&bHasNoData);
5426 82 : if (bHasNoData)
5427 14 : poDS->GetRasterBand(i)->SetNoDataValue(dfNoData);
5428 : }
5429 :
5430 82 : const double dfOffset = poSrcDS->GetRasterBand(i)->GetOffset();
5431 82 : if (dfOffset != 0.0)
5432 3 : poDS->GetRasterBand(i)->SetOffset(dfOffset);
5433 :
5434 82 : const double dfScale = poSrcDS->GetRasterBand(i)->GetScale();
5435 82 : if (dfScale != 1.0)
5436 3 : poDS->GetRasterBand(i)->SetScale(dfScale);
5437 :
5438 164 : poDS->GetRasterBand(i)->SetUnitType(
5439 82 : poSrcDS->GetRasterBand(i)->GetUnitType());
5440 : }
5441 :
5442 55 : if (poDS->m_bUseSrcLabel)
5443 : {
5444 53 : CSLConstList papszMD_PDS4 = poSrcDS->GetMetadata("xml:PDS4");
5445 53 : if (papszMD_PDS4 != nullptr)
5446 : {
5447 6 : poDS->SetMetadata(papszMD_PDS4, "xml:PDS4");
5448 : }
5449 : }
5450 :
5451 55 : if (poDS->m_poExternalDS == nullptr)
5452 : {
5453 : // We don't need to initialize the imagery as we are going to copy it
5454 : // completely
5455 46 : poDS->m_bMustInitImageFile = false;
5456 : }
5457 :
5458 55 : if (!CPLFetchBool(papszOptions, "CREATE_LABEL_ONLY", false))
5459 : {
5460 48 : CPLErr eErr = GDALDatasetCopyWholeRaster(poSrcDS, poDS.get(), nullptr,
5461 : pfnProgress, pProgressData);
5462 48 : poDS->FlushCache(false);
5463 48 : if (eErr != CE_None)
5464 : {
5465 1 : return nullptr;
5466 : }
5467 :
5468 47 : if (CPLFetchBool(papszOptions, "PROPAGATE_SRC_METADATA", true))
5469 : {
5470 46 : CSLConstList papszISIS3MD = poSrcDS->GetMetadata("json:ISIS3");
5471 46 : if (papszISIS3MD)
5472 : {
5473 2 : poDS->SetMetadata(papszISIS3MD, "json:ISIS3");
5474 :
5475 2 : if (poDS->m_poExternalDS)
5476 1 : poDS->m_poExternalDS->SetMetadata(papszISIS3MD,
5477 1 : "json:ISIS3");
5478 : }
5479 : }
5480 : }
5481 :
5482 54 : return poDS.release();
5483 : }
5484 :
5485 : /************************************************************************/
5486 : /* Delete() */
5487 : /************************************************************************/
5488 :
5489 109 : CPLErr PDS4Dataset::Delete(const char *pszFilename)
5490 :
5491 : {
5492 : /* -------------------------------------------------------------------- */
5493 : /* Collect file list. */
5494 : /* -------------------------------------------------------------------- */
5495 218 : GDALOpenInfo oOpenInfo(pszFilename, GA_ReadOnly);
5496 218 : auto poDS = PDS4Dataset::OpenInternal(&oOpenInfo);
5497 109 : if (poDS == nullptr)
5498 : {
5499 0 : if (CPLGetLastErrorNo() == 0)
5500 0 : CPLError(CE_Failure, CPLE_OpenFailed,
5501 : "Unable to open %s to obtain file list.", pszFilename);
5502 :
5503 0 : return CE_Failure;
5504 : }
5505 :
5506 109 : char **papszFileList = poDS->GetFileList();
5507 218 : CPLString osImageFilename = poDS->m_osImageFilename;
5508 : bool bCreatedFromExistingBinaryFile =
5509 109 : poDS->m_bCreatedFromExistingBinaryFile;
5510 :
5511 109 : poDS.reset();
5512 :
5513 109 : if (CSLCount(papszFileList) == 0)
5514 : {
5515 0 : CPLError(CE_Failure, CPLE_NotSupported,
5516 : "Unable to determine files associated with %s, "
5517 : "delete fails.",
5518 : pszFilename);
5519 0 : CSLDestroy(papszFileList);
5520 0 : return CE_Failure;
5521 : }
5522 :
5523 : /* -------------------------------------------------------------------- */
5524 : /* Delete all files. */
5525 : /* -------------------------------------------------------------------- */
5526 109 : CPLErr eErr = CE_None;
5527 392 : for (int i = 0; papszFileList[i] != nullptr; ++i)
5528 : {
5529 297 : if (bCreatedFromExistingBinaryFile &&
5530 14 : EQUAL(papszFileList[i], osImageFilename))
5531 : {
5532 7 : continue;
5533 : }
5534 276 : if (VSIUnlink(papszFileList[i]) != 0)
5535 : {
5536 0 : CPLError(CE_Failure, CPLE_AppDefined, "Deleting %s failed:\n%s",
5537 0 : papszFileList[i], VSIStrerror(errno));
5538 0 : eErr = CE_Failure;
5539 : }
5540 : }
5541 :
5542 109 : CSLDestroy(papszFileList);
5543 :
5544 109 : return eErr;
5545 : }
5546 :
5547 : /************************************************************************/
5548 : /* GDALRegister_PDS4() */
5549 : /************************************************************************/
5550 :
5551 2062 : void GDALRegister_PDS4()
5552 :
5553 : {
5554 2062 : if (GDALGetDriverByName(PDS4_DRIVER_NAME) != nullptr)
5555 263 : return;
5556 :
5557 1799 : GDALDriver *poDriver = new GDALDriver();
5558 1799 : PDS4DriverSetCommonMetadata(poDriver);
5559 :
5560 1799 : poDriver->pfnOpen = PDS4Dataset::Open;
5561 1799 : poDriver->pfnCreate = PDS4Dataset::Create;
5562 1799 : poDriver->pfnCreateCopy = PDS4Dataset::CreateCopy;
5563 1799 : poDriver->pfnDelete = PDS4Dataset::Delete;
5564 :
5565 1799 : GetGDALDriverManager()->RegisterDriver(poDriver);
5566 : }
|
