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