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