1 : /******************************************************************************
2 : * $Id $
3 : *
4 : * Project: SPOT Dimap Driver
5 : * Purpose: Implementation of SPOT Dimap driver.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : * Docs: http://www.spotimage.fr/dimap/spec/documentation/refdoc.htm
9 : *
10 : ******************************************************************************
11 : * Copyright (c) 2007, Frank Warmerdam <warmerdam@pobox.com>
12 : *
13 : * Permission is hereby granted, free of charge, to any person obtaining a
14 : * copy of this software and associated documentation files (the "Software"),
15 : * to deal in the Software without restriction, including without limitation
16 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
17 : * and/or sell copies of the Software, and to permit persons to whom the
18 : * Software is furnished to do so, subject to the following conditions:
19 : *
20 : * The above copyright notice and this permission notice shall be included
21 : * in all copies or substantial portions of the Software.
22 : *
23 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
24 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
25 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
26 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
27 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
28 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
29 : * DEALINGS IN THE SOFTWARE.
30 : ****************************************************************************/
31 :
32 : #include "gdal_pam.h"
33 : #include "cpl_minixml.h"
34 : #include "ogr_spatialref.h"
35 :
36 : CPL_CVSID("$Id: dimapdataset.cpp 17664 2009-09-21 21:16:45Z rouault $");
37 :
38 : CPL_C_START
39 : void GDALRegister_DIMAP(void);
40 : CPL_C_END
41 :
42 : /************************************************************************/
43 : /* ==================================================================== */
44 : /* DIMAPDataset */
45 : /* ==================================================================== */
46 : /************************************************************************/
47 :
48 : class DIMAPDataset : public GDALPamDataset
49 : {
50 : CPLXMLNode *psProduct;
51 :
52 : GDALDataset *poImageDS;
53 :
54 : int nGCPCount;
55 : GDAL_GCP *pasGCPList;
56 : char *pszGCPProjection;
57 :
58 : CPLString osProjection;
59 :
60 : int bHaveGeoTransform;
61 : double adfGeoTransform[6];
62 :
63 : char **papszXMLDimapMetadata;
64 :
65 : public:
66 : DIMAPDataset();
67 : ~DIMAPDataset();
68 :
69 : virtual const char *GetProjectionRef(void);
70 : virtual CPLErr GetGeoTransform( double * );
71 : virtual int GetGCPCount();
72 : virtual const char *GetGCPProjection();
73 : virtual const GDAL_GCP *GetGCPs();
74 : virtual char **GetMetadata( const char *pszDomain );
75 : virtual char **GetFileList(void);
76 :
77 : static int Identify( GDALOpenInfo * );
78 : static GDALDataset *Open( GDALOpenInfo * );
79 :
80 : CPLXMLNode *GetProduct() { return psProduct; }
81 : };
82 :
83 : /************************************************************************/
84 : /* ==================================================================== */
85 : /* DIMAPDataset */
86 : /* ==================================================================== */
87 : /************************************************************************/
88 :
89 : /************************************************************************/
90 : /* DIMAPDataset() */
91 : /************************************************************************/
92 :
93 2 : DIMAPDataset::DIMAPDataset()
94 : {
95 2 : psProduct = NULL;
96 :
97 2 : nGCPCount = 0;
98 2 : pasGCPList = NULL;
99 2 : pszGCPProjection = CPLStrdup("");
100 :
101 2 : poImageDS = NULL;
102 2 : bHaveGeoTransform = FALSE;
103 :
104 2 : papszXMLDimapMetadata = NULL;
105 2 : }
106 :
107 : /************************************************************************/
108 : /* ~DIMAPDataset() */
109 : /************************************************************************/
110 :
111 4 : DIMAPDataset::~DIMAPDataset()
112 :
113 : {
114 2 : FlushCache();
115 :
116 2 : CPLDestroyXMLNode( psProduct );
117 :
118 2 : CPLFree( pszGCPProjection );
119 2 : if( nGCPCount > 0 )
120 : {
121 2 : GDALDeinitGCPs( nGCPCount, pasGCPList );
122 2 : CPLFree( pasGCPList );
123 : }
124 :
125 2 : if( poImageDS != NULL )
126 2 : delete poImageDS;
127 :
128 2 : CSLDestroy(papszXMLDimapMetadata);
129 :
130 : /* -------------------------------------------------------------------- */
131 : /* Disconnect the bands so our destructor doesn't try and */
132 : /* delete them since they really belonged to poImageDS. */
133 : /* -------------------------------------------------------------------- */
134 : int iBand;
135 4 : for( iBand = 0; iBand < GetRasterCount(); iBand++ )
136 2 : papoBands[iBand] = NULL;
137 4 : }
138 :
139 : /************************************************************************/
140 : /* GetMetadata() */
141 : /* */
142 : /* We implement special support for fetching the full product */
143 : /* metadata as xml. */
144 : /************************************************************************/
145 :
146 2 : char **DIMAPDataset::GetMetadata( const char *pszDomain )
147 :
148 : {
149 2 : if( pszDomain && EQUAL(pszDomain,"xml:dimap") )
150 : {
151 0 : if (papszXMLDimapMetadata == NULL)
152 : {
153 0 : papszXMLDimapMetadata = (char **) CPLCalloc(sizeof(char*),2);
154 0 : papszXMLDimapMetadata[0] = CPLSerializeXMLTree( psProduct );
155 : }
156 0 : return papszXMLDimapMetadata;
157 : }
158 : else
159 2 : return GDALPamDataset::GetMetadata( pszDomain );
160 : }
161 :
162 : /************************************************************************/
163 : /* GetProjectionRef() */
164 : /************************************************************************/
165 :
166 0 : const char *DIMAPDataset::GetProjectionRef()
167 :
168 : {
169 0 : if( strlen(osProjection) > 0 )
170 0 : return osProjection;
171 : else
172 0 : return GDALPamDataset::GetProjectionRef();
173 : }
174 :
175 : /************************************************************************/
176 : /* GetGeoTransform() */
177 : /************************************************************************/
178 :
179 0 : CPLErr DIMAPDataset::GetGeoTransform( double *padfGeoTransform )
180 :
181 : {
182 0 : if( bHaveGeoTransform )
183 : {
184 0 : memcpy( padfGeoTransform, adfGeoTransform, sizeof(double)*6 );
185 0 : return CE_None;
186 : }
187 : else
188 0 : return GDALPamDataset::GetGeoTransform( padfGeoTransform );
189 : }
190 :
191 : /************************************************************************/
192 : /* GetFileList() */
193 : /************************************************************************/
194 :
195 0 : char **DIMAPDataset::GetFileList()
196 :
197 : {
198 0 : char **papszFileList = GDALPamDataset::GetFileList();
199 0 : char **papszImageFiles = poImageDS->GetFileList();
200 :
201 0 : papszFileList = CSLInsertStrings( papszFileList, -1, papszImageFiles );
202 :
203 0 : CSLDestroy( papszImageFiles );
204 :
205 0 : return papszFileList;
206 : }
207 :
208 : /************************************************************************/
209 : /* Identify() */
210 : /************************************************************************/
211 :
212 9592 : int DIMAPDataset::Identify( GDALOpenInfo * poOpenInfo )
213 :
214 : {
215 9592 : if( poOpenInfo->nHeaderBytes >= 100 )
216 : {
217 901 : if( strstr((const char *) poOpenInfo->pabyHeader,
218 : "<Dimap_Document" ) == NULL )
219 899 : return FALSE;
220 : else
221 2 : return TRUE;
222 : }
223 8691 : else if( poOpenInfo->bIsDirectory )
224 : {
225 : VSIStatBufL sStat;
226 :
227 : CPLString osMDFilename =
228 22 : CPLFormCIFilename( poOpenInfo->pszFilename, "METADATA.DIM", NULL );
229 :
230 22 : if( VSIStatL( osMDFilename, &sStat ) == 0 )
231 0 : return TRUE;
232 : else
233 22 : return FALSE;
234 : }
235 :
236 8669 : return FALSE;
237 : }
238 :
239 : /************************************************************************/
240 : /* Open() */
241 : /************************************************************************/
242 :
243 1863 : GDALDataset *DIMAPDataset::Open( GDALOpenInfo * poOpenInfo )
244 :
245 : {
246 1863 : if( !Identify( poOpenInfo ) )
247 1861 : return NULL;
248 :
249 : /* -------------------------------------------------------------------- */
250 : /* Confirm the requested access is supported. */
251 : /* -------------------------------------------------------------------- */
252 2 : if( poOpenInfo->eAccess == GA_Update )
253 : {
254 : CPLError( CE_Failure, CPLE_NotSupported,
255 : "The DIMAP driver does not support update access to existing"
256 0 : " datasets.\n" );
257 0 : return NULL;
258 : }
259 : /* -------------------------------------------------------------------- */
260 : /* Get the metadata filename. */
261 : /* -------------------------------------------------------------------- */
262 2 : CPLString osMDFilename;
263 :
264 2 : if( poOpenInfo->bIsDirectory )
265 : {
266 : osMDFilename =
267 0 : CPLFormCIFilename( poOpenInfo->pszFilename, "METADATA.DIM", NULL );
268 : }
269 : else
270 2 : osMDFilename = poOpenInfo->pszFilename;
271 :
272 : /* -------------------------------------------------------------------- */
273 : /* Ingest the xml file. */
274 : /* -------------------------------------------------------------------- */
275 : CPLXMLNode *psProduct, *psImageAttributes;
276 :
277 2 : psProduct = CPLParseXMLFile( osMDFilename );
278 2 : if( psProduct == NULL )
279 0 : return NULL;
280 :
281 2 : CPLXMLNode *psDoc = CPLGetXMLNode( psProduct, "=Dimap_Document" );
282 2 : psImageAttributes = CPLGetXMLNode( psDoc, "Raster_Dimensions" );
283 2 : if( psImageAttributes == NULL )
284 : {
285 : CPLError( CE_Failure, CPLE_OpenFailed,
286 0 : "Failed to find <Raster_Dimensions> in document." );
287 0 : return NULL;
288 : }
289 :
290 :
291 : /* -------------------------------------------------------------------- */
292 : /* Create the dataset. */
293 : /* -------------------------------------------------------------------- */
294 2 : DIMAPDataset *poDS = new DIMAPDataset();
295 :
296 2 : poDS->psProduct = psProduct;
297 :
298 : /* -------------------------------------------------------------------- */
299 : /* Get overall image information. */
300 : /* -------------------------------------------------------------------- */
301 : #ifdef DEBUG
302 : int nBands =
303 : atoi(CPLGetXMLValue( psImageAttributes, "NBANDS", "-1" ));
304 : #endif
305 :
306 : poDS->nRasterXSize =
307 4 : atoi(CPLGetXMLValue( psImageAttributes, "NCOLS", "-1" ));
308 : poDS->nRasterYSize =
309 2 : atoi(CPLGetXMLValue( psImageAttributes, "NROWS", "-1" ));
310 :
311 : /* -------------------------------------------------------------------- */
312 : /* Get the name of the underlying file. */
313 : /* -------------------------------------------------------------------- */
314 : const char *pszHref = CPLGetXMLValue(
315 2 : psDoc, "Data_Access.Data_File.DATA_FILE_PATH.href", "" );
316 2 : CPLString osPath = CPLGetPath(osMDFilename);
317 : CPLString osImageFilename =
318 2 : CPLFormFilename( osPath, pszHref, NULL );
319 :
320 : /* -------------------------------------------------------------------- */
321 : /* Try and open the file. */
322 : /* -------------------------------------------------------------------- */
323 2 : poDS->poImageDS = (GDALDataset *) GDALOpen( osImageFilename, GA_ReadOnly );
324 2 : if( poDS->poImageDS == NULL )
325 : {
326 0 : delete poDS;
327 0 : return NULL;
328 : }
329 :
330 : /* -------------------------------------------------------------------- */
331 : /* Attach the bands. */
332 : /* -------------------------------------------------------------------- */
333 : int iBand;
334 : CPLAssert( nBands == poDS->poImageDS->GetRasterCount() );
335 :
336 4 : for( iBand = 1; iBand <= poDS->poImageDS->GetRasterCount(); iBand++ )
337 2 : poDS->SetBand( iBand, poDS->poImageDS->GetRasterBand( iBand ) );
338 :
339 : /* -------------------------------------------------------------------- */
340 : /* Try to collect simple insertion point. */
341 : /* -------------------------------------------------------------------- */
342 : CPLXMLNode *psGeoLoc =
343 2 : CPLGetXMLNode( psDoc, "Geoposition.Geoposition_Insert" );
344 :
345 2 : if( psGeoLoc != NULL )
346 : {
347 0 : poDS->bHaveGeoTransform = TRUE;
348 0 : poDS->adfGeoTransform[0] = atof(CPLGetXMLValue(psGeoLoc,"ULXMAP","0"));
349 0 : poDS->adfGeoTransform[1] = atof(CPLGetXMLValue(psGeoLoc,"XDIM","0"));
350 0 : poDS->adfGeoTransform[2] = 0.0;
351 0 : poDS->adfGeoTransform[3] = atof(CPLGetXMLValue(psGeoLoc,"ULYMAP","0"));
352 0 : poDS->adfGeoTransform[4] = 0.0;
353 0 : poDS->adfGeoTransform[5] = -atof(CPLGetXMLValue(psGeoLoc,"YDIM","0"));
354 : }
355 :
356 : /* -------------------------------------------------------------------- */
357 : /* Collect GCPs. */
358 : /* -------------------------------------------------------------------- */
359 2 : psGeoLoc = CPLGetXMLNode( psDoc, "Geoposition.Geoposition_Points" );
360 :
361 2 : if( psGeoLoc != NULL )
362 : {
363 : CPLXMLNode *psNode;
364 :
365 : // count gcps.
366 2 : poDS->nGCPCount = 0;
367 10 : for( psNode = psGeoLoc->psChild; psNode != NULL;
368 : psNode = psNode->psNext )
369 : {
370 8 : if( EQUAL(psNode->pszValue,"Tie_Point") )
371 8 : poDS->nGCPCount++ ;
372 : }
373 :
374 : poDS->pasGCPList = (GDAL_GCP *)
375 2 : CPLCalloc(sizeof(GDAL_GCP),poDS->nGCPCount);
376 :
377 2 : poDS->nGCPCount = 0;
378 :
379 10 : for( psNode = psGeoLoc->psChild; psNode != NULL;
380 : psNode = psNode->psNext )
381 : {
382 : char szID[32];
383 8 : GDAL_GCP *psGCP = poDS->pasGCPList + poDS->nGCPCount;
384 :
385 8 : if( !EQUAL(psNode->pszValue,"Tie_Point") )
386 0 : continue;
387 :
388 8 : poDS->nGCPCount++ ;
389 :
390 8 : sprintf( szID, "%d", poDS->nGCPCount );
391 8 : psGCP->pszId = CPLStrdup( szID );
392 8 : psGCP->pszInfo = CPLStrdup("");
393 : psGCP->dfGCPPixel =
394 8 : atof(CPLGetXMLValue(psNode,"TIE_POINT_DATA_X","0"))-0.5;
395 : psGCP->dfGCPLine =
396 8 : atof(CPLGetXMLValue(psNode,"TIE_POINT_DATA_Y","0"))-0.5;
397 : psGCP->dfGCPX =
398 8 : atof(CPLGetXMLValue(psNode,"TIE_POINT_CRS_X",""));
399 : psGCP->dfGCPY =
400 8 : atof(CPLGetXMLValue(psNode,"TIE_POINT_CRS_Y",""));
401 : psGCP->dfGCPZ =
402 8 : atof(CPLGetXMLValue(psNode,"TIE_POINT_CRS_Z",""));
403 : }
404 : }
405 :
406 : /* -------------------------------------------------------------------- */
407 : /* Collect the CRS. For now we look only for EPSG codes. */
408 : /* -------------------------------------------------------------------- */
409 : const char *pszSRS = CPLGetXMLValue(
410 : psDoc,
411 : "Coordinate_Reference_System.Horizontal_CS.HORIZONTAL_CS_CODE",
412 2 : NULL );
413 :
414 2 : if( pszSRS != NULL )
415 : {
416 2 : OGRSpatialReference oSRS;
417 2 : if( oSRS.SetFromUserInput( pszSRS ) == OGRERR_NONE )
418 : {
419 2 : if( poDS->nGCPCount > 0 )
420 : {
421 2 : CPLFree(poDS->pszGCPProjection);
422 2 : oSRS.exportToWkt( &(poDS->pszGCPProjection) );
423 : }
424 : else
425 : {
426 0 : char *pszProjection = NULL;
427 0 : oSRS.exportToWkt( &pszProjection );
428 0 : poDS->osProjection = pszProjection;
429 0 : CPLFree( pszProjection );
430 : }
431 2 : }
432 : }
433 :
434 : /* -------------------------------------------------------------------- */
435 : /* Translate other metadata of interest. */
436 : /* -------------------------------------------------------------------- */
437 : static const char *apszMetadataTranslation[] =
438 : {
439 : "Production", "",
440 : "Production.Facility", "FACILITY_",
441 : "Dataset_Sources.Source_Information.Scene_Source", "",
442 : "Data_Processing", "",
443 : "Image_Interpretation.Spectral_Band_Info", "SPECTRAL_",
444 : NULL, NULL
445 : };
446 :
447 : int iTrItem;
448 :
449 12 : for( iTrItem = 0; apszMetadataTranslation[iTrItem] != NULL; iTrItem += 2 )
450 : {
451 : CPLXMLNode *psParent =
452 10 : CPLGetXMLNode( psDoc, apszMetadataTranslation[iTrItem] );
453 :
454 10 : if( psParent == NULL )
455 2 : continue;
456 :
457 : // hackey logic to support directly access a name/value entry
458 : // or a parent element with many name/values.
459 :
460 : CPLXMLNode *psTarget;
461 8 : if( psParent->psChild != NULL
462 : && psParent->psChild->eType == CXT_Text )
463 0 : psTarget = psParent;
464 : else
465 8 : psTarget = psParent->psChild;
466 :
467 70 : for( ; psTarget != NULL && psTarget != psParent;
468 : psTarget = psTarget->psNext )
469 : {
470 62 : if( psTarget->eType == CXT_Element
471 : && psTarget->psChild != NULL
472 : && psTarget->psChild->eType == CXT_Text )
473 : {
474 54 : CPLString osName = apszMetadataTranslation[iTrItem+1];
475 :
476 54 : osName += psTarget->pszValue;
477 54 : poDS->SetMetadataItem( osName, psTarget->psChild->pszValue );
478 : }
479 : }
480 : }
481 :
482 : /* -------------------------------------------------------------------- */
483 : /* Set Band metadata from the <Spectral_Band_Info> content */
484 : /* -------------------------------------------------------------------- */
485 :
486 : CPLXMLNode *psImageInterpretationNode =
487 2 : CPLGetXMLNode( psDoc, "Image_Interpretation" );
488 2 : if (psImageInterpretationNode != NULL)
489 : {
490 2 : CPLXMLNode *psSpectralBandInfoNode = psImageInterpretationNode->psChild;
491 6 : while (psSpectralBandInfoNode != NULL)
492 : {
493 2 : if (psSpectralBandInfoNode->eType == CXT_Element &&
494 : EQUAL(psSpectralBandInfoNode->pszValue, "Spectral_Band_Info"))
495 : {
496 2 : CPLXMLNode *psTag = psSpectralBandInfoNode->psChild;
497 2 : int nBandIndex = 0;
498 16 : while(psTag != NULL)
499 : {
500 12 : if (psTag->eType == CXT_Element && psTag->psChild != NULL &&
501 : psTag->psChild->eType == CXT_Text && psTag->pszValue != NULL)
502 : {
503 12 : if (EQUAL(psTag->pszValue, "BAND_INDEX"))
504 : {
505 2 : nBandIndex = atoi(psTag->psChild->pszValue);
506 2 : if (nBandIndex <= 0 ||
507 : nBandIndex > poDS->poImageDS->GetRasterCount())
508 : {
509 : CPLError(CE_Warning, CPLE_AppDefined,
510 0 : "Bad BAND_INDEX value : %s", psTag->psChild->pszValue);
511 0 : nBandIndex = 0;
512 : }
513 : }
514 10 : else if (nBandIndex >= 1)
515 : {
516 : poDS->GetRasterBand(nBandIndex)->SetMetadataItem(
517 10 : psTag->pszValue, psTag->psChild->pszValue);
518 : }
519 : }
520 12 : psTag = psTag->psNext;
521 : }
522 : }
523 2 : psSpectralBandInfoNode = psSpectralBandInfoNode->psNext;
524 : }
525 : }
526 :
527 : /* -------------------------------------------------------------------- */
528 : /* Initialize any PAM information. */
529 : /* -------------------------------------------------------------------- */
530 2 : poDS->SetDescription( osMDFilename );
531 2 : poDS->TryLoadXML();
532 :
533 : /* -------------------------------------------------------------------- */
534 : /* Check for overviews. */
535 : /* -------------------------------------------------------------------- */
536 2 : poDS->oOvManager.Initialize( poDS, osMDFilename );
537 :
538 2 : return( poDS );
539 : }
540 :
541 : /************************************************************************/
542 : /* GetGCPCount() */
543 : /************************************************************************/
544 :
545 1 : int DIMAPDataset::GetGCPCount()
546 :
547 : {
548 1 : return nGCPCount;
549 : }
550 :
551 : /************************************************************************/
552 : /* GetGCPProjection() */
553 : /************************************************************************/
554 :
555 1 : const char *DIMAPDataset::GetGCPProjection()
556 :
557 : {
558 1 : return pszGCPProjection;
559 : }
560 :
561 : /************************************************************************/
562 : /* GetGCPs() */
563 : /************************************************************************/
564 :
565 1 : const GDAL_GCP *DIMAPDataset::GetGCPs()
566 :
567 : {
568 1 : return pasGCPList;
569 : }
570 :
571 : /************************************************************************/
572 : /* GDALRegister_DIMAP() */
573 : /************************************************************************/
574 :
575 338 : void GDALRegister_DIMAP()
576 :
577 : {
578 : GDALDriver *poDriver;
579 :
580 338 : if( GDALGetDriverByName( "DIMAP" ) == NULL )
581 : {
582 336 : poDriver = new GDALDriver();
583 :
584 336 : poDriver->SetDescription( "DIMAP" );
585 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
586 336 : "SPOT DIMAP" );
587 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
588 336 : "frmt_various.html#DIMAP" );
589 :
590 336 : poDriver->pfnOpen = DIMAPDataset::Open;
591 336 : poDriver->pfnIdentify = DIMAPDataset::Identify;
592 :
593 336 : GetGDALDriverManager()->RegisterDriver( poDriver );
594 : }
595 338 : }
596 :
|
