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