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