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