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