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