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