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