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