1 : /******************************************************************************
2 : * $Id: ersdataset.cpp 24927 2012-09-16 12:14:44Z rouault $
3 : *
4 : * Project: ERMapper .ers Driver
5 : * Purpose: Implementation of .ers driver.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2007, Frank Warmerdam <warmerdam@pobox.com>
10 : *
11 : * Permission is hereby granted, free of charge, to any person obtaining a
12 : * copy of this software and associated documentation files (the "Software"),
13 : * to deal in the Software without restriction, including without limitation
14 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
15 : * and/or sell copies of the Software, and to permit persons to whom the
16 : * Software is furnished to do so, subject to the following conditions:
17 : *
18 : * The above copyright notice and this permission notice shall be included
19 : * in all copies or substantial portions of the Software.
20 : *
21 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
22 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
23 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
24 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
25 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
26 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
27 : * DEALINGS IN THE SOFTWARE.
28 : ****************************************************************************/
29 :
30 : #include "rawdataset.h"
31 : #include "ogr_spatialref.h"
32 : #include "cpl_string.h"
33 : #include "ershdrnode.h"
34 :
35 : CPL_CVSID("$Id: ersdataset.cpp 24927 2012-09-16 12:14:44Z rouault $");
36 :
37 : /************************************************************************/
38 : /* ==================================================================== */
39 : /* ERSDataset */
40 : /* ==================================================================== */
41 : /************************************************************************/
42 :
43 : class ERSRasterBand;
44 :
45 : class ERSDataset : public RawDataset
46 : {
47 : friend class ERSRasterBand;
48 :
49 : VSILFILE *fpImage; // image data file.
50 : GDALDataset *poDepFile;
51 :
52 : int bGotTransform;
53 : double adfGeoTransform[6];
54 : char *pszProjection;
55 :
56 : CPLString osRawFilename;
57 :
58 : int bHDRDirty;
59 : ERSHdrNode *poHeader;
60 :
61 : const char *Find( const char *, const char * );
62 :
63 : int nGCPCount;
64 : GDAL_GCP *pasGCPList;
65 : char *pszGCPProjection;
66 :
67 : void ReadGCPs();
68 :
69 : int bHasNoDataValue;
70 : double dfNoDataValue;
71 :
72 : CPLString osProj;
73 : CPLString osDatum;
74 : CPLString osUnits;
75 : void WriteProjectionInfo(const char* pszProj,
76 : const char* pszDatum,
77 : const char* pszUnits);
78 :
79 : CPLStringList oERSMetadataList;
80 :
81 : protected:
82 : virtual int CloseDependentDatasets();
83 :
84 : public:
85 : ERSDataset();
86 : ~ERSDataset();
87 :
88 : virtual void FlushCache(void);
89 : virtual CPLErr GetGeoTransform( double * padfTransform );
90 : virtual CPLErr SetGeoTransform( double *padfTransform );
91 : virtual const char *GetProjectionRef(void);
92 : virtual CPLErr SetProjection( const char * );
93 : virtual char **GetFileList(void);
94 :
95 : virtual int GetGCPCount();
96 : virtual const char *GetGCPProjection();
97 : virtual const GDAL_GCP *GetGCPs();
98 : virtual CPLErr SetGCPs( int nGCPCount, const GDAL_GCP *pasGCPList,
99 : const char *pszGCPProjection );
100 :
101 : virtual const char *GetMetadataItem( const char * pszName,
102 : const char * pszDomain = "" );
103 : virtual char **GetMetadata( const char * pszDomain = "" );
104 :
105 : static GDALDataset *Open( GDALOpenInfo * );
106 : static int Identify( GDALOpenInfo * );
107 : static GDALDataset *Create( const char * pszFilename,
108 : int nXSize, int nYSize, int nBands,
109 : GDALDataType eType, char ** papszParmList );
110 : };
111 :
112 : /************************************************************************/
113 : /* ERSDataset() */
114 : /************************************************************************/
115 :
116 50 : ERSDataset::ERSDataset()
117 : {
118 50 : fpImage = NULL;
119 50 : poDepFile = NULL;
120 50 : pszProjection = CPLStrdup("");
121 50 : bGotTransform = FALSE;
122 50 : adfGeoTransform[0] = 0.0;
123 50 : adfGeoTransform[1] = 1.0;
124 50 : adfGeoTransform[2] = 0.0;
125 50 : adfGeoTransform[3] = 0.0;
126 50 : adfGeoTransform[4] = 0.0;
127 50 : adfGeoTransform[5] = 1.0;
128 50 : poHeader = NULL;
129 50 : bHDRDirty = FALSE;
130 :
131 50 : nGCPCount = 0;
132 50 : pasGCPList = NULL;
133 50 : pszGCPProjection = CPLStrdup("");
134 :
135 50 : bHasNoDataValue = FALSE;
136 50 : dfNoDataValue = 0.0;
137 50 : }
138 :
139 : /************************************************************************/
140 : /* ~ERSDataset() */
141 : /************************************************************************/
142 :
143 50 : ERSDataset::~ERSDataset()
144 :
145 : {
146 50 : FlushCache();
147 :
148 50 : if( fpImage != NULL )
149 : {
150 50 : VSIFCloseL( fpImage );
151 : }
152 :
153 50 : CloseDependentDatasets();
154 :
155 50 : CPLFree( pszProjection );
156 :
157 50 : CPLFree( pszGCPProjection );
158 50 : if( nGCPCount > 0 )
159 : {
160 3 : GDALDeinitGCPs( nGCPCount, pasGCPList );
161 3 : CPLFree( pasGCPList );
162 : }
163 :
164 50 : if( poHeader != NULL )
165 50 : delete poHeader;
166 50 : }
167 :
168 : /************************************************************************/
169 : /* CloseDependentDatasets() */
170 : /************************************************************************/
171 :
172 50 : int ERSDataset::CloseDependentDatasets()
173 : {
174 50 : int bHasDroppedRef = RawDataset::CloseDependentDatasets();
175 :
176 50 : if( poDepFile != NULL )
177 : {
178 : int iBand;
179 :
180 0 : bHasDroppedRef = TRUE;
181 :
182 0 : for( iBand = 0; iBand < nBands; iBand++ )
183 0 : papoBands[iBand] = NULL;
184 0 : nBands = 0;
185 :
186 0 : GDALClose( (GDALDatasetH) poDepFile );
187 0 : poDepFile = NULL;
188 : }
189 :
190 50 : return bHasDroppedRef;
191 : }
192 :
193 : /************************************************************************/
194 : /* FlushCache() */
195 : /************************************************************************/
196 :
197 50 : void ERSDataset::FlushCache()
198 :
199 : {
200 50 : if( bHDRDirty )
201 : {
202 33 : VSILFILE * fpERS = VSIFOpenL( GetDescription(), "w" );
203 33 : if( fpERS == NULL )
204 : {
205 : CPLError( CE_Failure, CPLE_OpenFailed,
206 : "Unable to rewrite %s header.",
207 0 : GetDescription() );
208 : }
209 : else
210 : {
211 33 : VSIFPrintfL( fpERS, "DatasetHeader Begin\n" );
212 33 : poHeader->WriteSelf( fpERS, 1 );
213 33 : VSIFPrintfL( fpERS, "DatasetHeader End\n" );
214 33 : VSIFCloseL( fpERS );
215 : }
216 : }
217 :
218 50 : RawDataset::FlushCache();
219 50 : }
220 :
221 : /************************************************************************/
222 : /* GetMetadataItem() */
223 : /************************************************************************/
224 :
225 47 : const char *ERSDataset::GetMetadataItem( const char * pszName,
226 : const char * pszDomain )
227 : {
228 47 : if (pszDomain != NULL && EQUAL(pszDomain, "ERS") && pszName != NULL)
229 : {
230 3 : if (EQUAL(pszName, "PROJ"))
231 1 : return osProj.size() ? osProj.c_str() : NULL;
232 2 : if (EQUAL(pszName, "DATUM"))
233 1 : return osDatum.size() ? osDatum.c_str() : NULL;
234 1 : if (EQUAL(pszName, "UNITS"))
235 1 : return osUnits.size() ? osUnits.c_str() : NULL;
236 : }
237 44 : return GDALPamDataset::GetMetadataItem(pszName, pszDomain);
238 : }
239 :
240 : /************************************************************************/
241 : /* GetMetadata() */
242 : /************************************************************************/
243 :
244 3 : char **ERSDataset::GetMetadata( const char *pszDomain )
245 :
246 : {
247 3 : if( pszDomain != NULL && EQUAL(pszDomain, "ERS") )
248 : {
249 1 : oERSMetadataList.Clear();
250 1 : if (osProj.size())
251 1 : oERSMetadataList.AddString(CPLSPrintf("%s=%s", "PROJ", osProj.c_str()));
252 1 : if (osDatum.size())
253 1 : oERSMetadataList.AddString(CPLSPrintf("%s=%s", "DATUM", osDatum.c_str()));
254 1 : if (osUnits.size())
255 1 : oERSMetadataList.AddString(CPLSPrintf("%s=%s", "UNITS", osUnits.c_str()));
256 1 : return oERSMetadataList.List();
257 : }
258 : else
259 2 : return GDALPamDataset::GetMetadata( pszDomain );
260 : }
261 :
262 : /************************************************************************/
263 : /* GetGCPCount() */
264 : /************************************************************************/
265 :
266 3 : int ERSDataset::GetGCPCount()
267 :
268 : {
269 3 : return nGCPCount;
270 : }
271 :
272 : /************************************************************************/
273 : /* GetGCPProjection() */
274 : /************************************************************************/
275 :
276 1 : const char *ERSDataset::GetGCPProjection()
277 :
278 : {
279 1 : return pszGCPProjection;
280 : }
281 :
282 : /************************************************************************/
283 : /* GetGCPs() */
284 : /************************************************************************/
285 :
286 1 : const GDAL_GCP *ERSDataset::GetGCPs()
287 :
288 : {
289 1 : return pasGCPList;
290 : }
291 :
292 : /************************************************************************/
293 : /* SetGCPs() */
294 : /************************************************************************/
295 :
296 1 : CPLErr ERSDataset::SetGCPs( int nGCPCountIn, const GDAL_GCP *pasGCPListIn,
297 : const char *pszGCPProjectionIn )
298 :
299 : {
300 : /* -------------------------------------------------------------------- */
301 : /* Clean old gcps. */
302 : /* -------------------------------------------------------------------- */
303 1 : CPLFree( pszGCPProjection );
304 1 : pszGCPProjection = NULL;
305 :
306 1 : if( nGCPCount > 0 )
307 : {
308 0 : GDALDeinitGCPs( nGCPCount, pasGCPList );
309 0 : CPLFree( pasGCPList );
310 :
311 0 : pasGCPList = NULL;
312 0 : nGCPCount = 0;
313 : }
314 :
315 : /* -------------------------------------------------------------------- */
316 : /* Copy new ones. */
317 : /* -------------------------------------------------------------------- */
318 1 : nGCPCount = nGCPCountIn;
319 1 : pasGCPList = GDALDuplicateGCPs( nGCPCount, pasGCPListIn );
320 1 : pszGCPProjection = CPLStrdup( pszGCPProjectionIn );
321 :
322 : /* -------------------------------------------------------------------- */
323 : /* Setup the header contents corresponding to these GCPs. */
324 : /* -------------------------------------------------------------------- */
325 1 : bHDRDirty = TRUE;
326 :
327 1 : poHeader->Set( "RasterInfo.WarpControl.WarpType", "Polynomial" );
328 1 : if( nGCPCount > 6 )
329 0 : poHeader->Set( "RasterInfo.WarpControl.WarpOrder", "2" );
330 : else
331 1 : poHeader->Set( "RasterInfo.WarpControl.WarpOrder", "1" );
332 1 : poHeader->Set( "RasterInfo.WarpControl.WarpSampling", "Nearest" );
333 :
334 : /* -------------------------------------------------------------------- */
335 : /* Translate the projection. */
336 : /* -------------------------------------------------------------------- */
337 1 : OGRSpatialReference oSRS( pszGCPProjection );
338 : char szERSProj[32], szERSDatum[32], szERSUnits[32];
339 :
340 1 : oSRS.exportToERM( szERSProj, szERSDatum, szERSUnits );
341 :
342 : /* Write the above computed values, unless they have been overriden by */
343 : /* the creation options PROJ, DATUM or UNITS */
344 :
345 : poHeader->Set( "RasterInfo.WarpControl.CoordinateSpace.Datum",
346 : CPLString().Printf( "\"%s\"",
347 1 : (osDatum.size()) ? osDatum.c_str() : szERSDatum ) );
348 : poHeader->Set( "RasterInfo.WarpControl.CoordinateSpace.Projection",
349 : CPLString().Printf( "\"%s\"",
350 1 : (osProj.size()) ? osProj.c_str() : szERSProj ) );
351 : poHeader->Set( "RasterInfo.WarpControl.CoordinateSpace.CoordinateType",
352 1 : CPLString().Printf( "EN" ) );
353 : poHeader->Set( "RasterInfo.WarpControl.CoordinateSpace.Units",
354 : CPLString().Printf( "\"%s\"",
355 1 : (osUnits.size()) ? osUnits.c_str() : szERSUnits ) );
356 : poHeader->Set( "RasterInfo.WarpControl.CoordinateSpace.Rotation",
357 1 : "0:0:0.0" );
358 :
359 : /* -------------------------------------------------------------------- */
360 : /* Translate the GCPs. */
361 : /* -------------------------------------------------------------------- */
362 1 : CPLString osControlPoints = "{\n";
363 : int iGCP;
364 :
365 5 : for( iGCP = 0; iGCP < nGCPCount; iGCP++ )
366 : {
367 4 : CPLString osLine;
368 :
369 4 : CPLString osId = pasGCPList[iGCP].pszId;
370 4 : if( strlen(osId) == 0 )
371 4 : osId.Printf( "%d", iGCP + 1 );
372 :
373 : osLine.Printf( "\t\t\t\t\"%s\"\tYes\tYes\t%.6f\t%.6f\t%.15g\t%.15g\t%.15g\n",
374 : osId.c_str(),
375 4 : pasGCPList[iGCP].dfGCPPixel,
376 4 : pasGCPList[iGCP].dfGCPLine,
377 4 : pasGCPList[iGCP].dfGCPX,
378 4 : pasGCPList[iGCP].dfGCPY,
379 20 : pasGCPList[iGCP].dfGCPZ );
380 4 : osControlPoints += osLine;
381 : }
382 1 : osControlPoints += "\t\t}";
383 :
384 1 : poHeader->Set( "RasterInfo.WarpControl.ControlPoints", osControlPoints );
385 :
386 1 : return CE_None;
387 : }
388 :
389 : /************************************************************************/
390 : /* GetProjectionRef() */
391 : /************************************************************************/
392 :
393 5 : const char *ERSDataset::GetProjectionRef()
394 :
395 : {
396 : // try xml first
397 5 : const char* pszPrj = GDALPamDataset::GetProjectionRef();
398 5 : if(pszPrj && strlen(pszPrj) > 0)
399 0 : return pszPrj;
400 :
401 5 : return pszProjection;
402 : }
403 :
404 : /************************************************************************/
405 : /* SetProjection() */
406 : /************************************************************************/
407 :
408 31 : CPLErr ERSDataset::SetProjection( const char *pszSRS )
409 :
410 : {
411 31 : if( pszProjection && EQUAL(pszSRS,pszProjection) )
412 0 : return CE_None;
413 :
414 31 : if( pszSRS == NULL )
415 0 : pszSRS = "";
416 :
417 31 : CPLFree( pszProjection );
418 31 : pszProjection = CPLStrdup(pszSRS);
419 :
420 31 : OGRSpatialReference oSRS( pszSRS );
421 : char szERSProj[32], szERSDatum[32], szERSUnits[32];
422 :
423 31 : oSRS.exportToERM( szERSProj, szERSDatum, szERSUnits );
424 :
425 : /* Write the above computed values, unless they have been overriden by */
426 : /* the creation options PROJ, DATUM or UNITS */
427 : WriteProjectionInfo( (osProj.size()) ? osProj.c_str() : szERSProj,
428 : (osDatum.size()) ? osDatum.c_str() : szERSDatum,
429 31 : (osUnits.size()) ? osUnits.c_str() : szERSUnits );
430 :
431 31 : return CE_None;
432 : }
433 :
434 : /************************************************************************/
435 : /* WriteProjectionInfo() */
436 : /************************************************************************/
437 :
438 32 : void ERSDataset::WriteProjectionInfo(const char* pszProj,
439 : const char* pszDatum,
440 : const char* pszUnits)
441 : {
442 32 : bHDRDirty = TRUE;
443 : poHeader->Set( "CoordinateSpace.Datum",
444 32 : CPLString().Printf( "\"%s\"", pszDatum ) );
445 : poHeader->Set( "CoordinateSpace.Projection",
446 64 : CPLString().Printf( "\"%s\"", pszProj ) );
447 : poHeader->Set( "CoordinateSpace.CoordinateType",
448 64 : CPLString().Printf( "EN" ) );
449 : poHeader->Set( "CoordinateSpace.Units",
450 64 : CPLString().Printf( "\"%s\"", pszUnits ) );
451 : poHeader->Set( "CoordinateSpace.Rotation",
452 32 : "0:0:0.0" );
453 :
454 : /* -------------------------------------------------------------------- */
455 : /* It seems that CoordinateSpace needs to come before */
456 : /* RasterInfo. Try moving it up manually. */
457 : /* -------------------------------------------------------------------- */
458 32 : int iRasterInfo = -1;
459 32 : int iCoordSpace = -1;
460 : int i;
461 :
462 223 : for( i = 0; i < poHeader->nItemCount; i++ )
463 : {
464 223 : if( EQUAL(poHeader->papszItemName[i],"RasterInfo") )
465 31 : iRasterInfo = i;
466 :
467 223 : if( EQUAL(poHeader->papszItemName[i],"CoordinateSpace") )
468 : {
469 32 : iCoordSpace = i;
470 32 : break;
471 : }
472 : }
473 :
474 32 : if( iCoordSpace > iRasterInfo && iRasterInfo != -1 )
475 : {
476 62 : for( i = iCoordSpace; i > 0 && i != iRasterInfo; i-- )
477 : {
478 : char *pszTemp;
479 :
480 31 : ERSHdrNode *poTemp = poHeader->papoItemChild[i];
481 31 : poHeader->papoItemChild[i] = poHeader->papoItemChild[i-1];
482 31 : poHeader->papoItemChild[i-1] = poTemp;
483 :
484 31 : pszTemp = poHeader->papszItemName[i];
485 31 : poHeader->papszItemName[i] = poHeader->papszItemName[i-1];
486 31 : poHeader->papszItemName[i-1] = pszTemp;
487 :
488 31 : pszTemp = poHeader->papszItemValue[i];
489 31 : poHeader->papszItemValue[i] = poHeader->papszItemValue[i-1];
490 31 : poHeader->papszItemValue[i-1] = pszTemp;
491 : }
492 : }
493 32 : }
494 :
495 : /************************************************************************/
496 : /* GetGeoTransform() */
497 : /************************************************************************/
498 :
499 20 : CPLErr ERSDataset::GetGeoTransform( double * padfTransform )
500 :
501 : {
502 20 : if( bGotTransform )
503 : {
504 20 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
505 20 : return CE_None;
506 : }
507 : else
508 : {
509 0 : return GDALPamDataset::GetGeoTransform( padfTransform );
510 : }
511 : }
512 :
513 : /************************************************************************/
514 : /* SetGeoTransform() */
515 : /************************************************************************/
516 :
517 30 : CPLErr ERSDataset::SetGeoTransform( double *padfTransform )
518 :
519 : {
520 30 : if( memcmp( padfTransform, adfGeoTransform, sizeof(double)*6 ) == 0 )
521 0 : return CE_None;
522 :
523 30 : if( adfGeoTransform[2] != 0 || adfGeoTransform[4] != 0 )
524 : {
525 : CPLError( CE_Failure, CPLE_AppDefined,
526 0 : "Rotated and skewed geotransforms not currently supported for ERS driver." );
527 0 : return CE_Failure;
528 : }
529 :
530 30 : bGotTransform = TRUE;
531 30 : memcpy( adfGeoTransform, padfTransform, sizeof(double) * 6 );
532 :
533 30 : bHDRDirty = TRUE;
534 :
535 : poHeader->Set( "RasterInfo.CellInfo.Xdimension",
536 30 : CPLString().Printf( "%.15g", fabs(adfGeoTransform[1]) ) );
537 : poHeader->Set( "RasterInfo.CellInfo.Ydimension",
538 60 : CPLString().Printf( "%.15g", fabs(adfGeoTransform[5]) ) );
539 : poHeader->Set( "RasterInfo.RegistrationCoord.Eastings",
540 60 : CPLString().Printf( "%.15g", adfGeoTransform[0] ) );
541 : poHeader->Set( "RasterInfo.RegistrationCoord.Northings",
542 60 : CPLString().Printf( "%.15g", adfGeoTransform[3] ) );
543 :
544 30 : return CE_None;
545 : }
546 :
547 : /************************************************************************/
548 : /* ERSDMS2Dec() */
549 : /* */
550 : /* Convert ERS DMS format to decimal degrees. Input is like */
551 : /* "-180:00:00". */
552 : /************************************************************************/
553 :
554 4 : static double ERSDMS2Dec( const char *pszDMS )
555 :
556 : {
557 4 : char **papszTokens = CSLTokenizeStringComplex( pszDMS, ":", FALSE, FALSE );
558 :
559 4 : if( CSLCount(papszTokens) != 3 )
560 : {
561 0 : CSLDestroy(papszTokens);
562 0 : return CPLAtof( pszDMS );
563 : }
564 : else
565 : {
566 4 : double dfResult = fabs(CPLAtof(papszTokens[0]))
567 4 : + CPLAtof(papszTokens[1]) / 60.0
568 8 : + CPLAtof(papszTokens[2]) / 3600.0;
569 :
570 4 : if( CPLAtof(papszTokens[0]) < 0 )
571 3 : dfResult *= -1;
572 :
573 4 : CSLDestroy( papszTokens );
574 4 : return dfResult;
575 : }
576 : }
577 :
578 : /************************************************************************/
579 : /* GetFileList() */
580 : /************************************************************************/
581 :
582 6 : char **ERSDataset::GetFileList()
583 :
584 : {
585 6 : char **papszFileList = NULL;
586 :
587 : // Main data file, etc.
588 6 : papszFileList = GDALPamDataset::GetFileList();
589 :
590 : // Add raw data file if we have one.
591 6 : if( strlen(osRawFilename) > 0 )
592 6 : papszFileList = CSLAddString( papszFileList, osRawFilename );
593 :
594 : // If we have a dependent file, merge it's list of files in.
595 6 : if( poDepFile )
596 : {
597 0 : char **papszDepFiles = poDepFile->GetFileList();
598 : papszFileList =
599 0 : CSLInsertStrings( papszFileList, -1, papszDepFiles );
600 0 : CSLDestroy( papszDepFiles );
601 : }
602 :
603 6 : return papszFileList;
604 : }
605 :
606 : /************************************************************************/
607 : /* ReadGCPs() */
608 : /* */
609 : /* Read the GCPs from the header. */
610 : /************************************************************************/
611 :
612 2 : void ERSDataset::ReadGCPs()
613 :
614 : {
615 : const char *pszCP =
616 2 : poHeader->Find( "RasterInfo.WarpControl.ControlPoints", NULL );
617 :
618 2 : if( pszCP == NULL )
619 0 : return;
620 :
621 : /* -------------------------------------------------------------------- */
622 : /* Parse the control points. They will look something like: */
623 : /* */
624 : /* "1035" Yes No 2344.650885 3546.419458 483270.73 3620906.21 3.105 */
625 : /* -------------------------------------------------------------------- */
626 2 : char **papszTokens = CSLTokenizeStringComplex( pszCP, "{ \t}", TRUE,FALSE);
627 : int nItemsPerLine;
628 2 : int nItemCount = CSLCount(papszTokens);
629 :
630 : /* -------------------------------------------------------------------- */
631 : /* Work out if we have elevation values or not. */
632 : /* -------------------------------------------------------------------- */
633 2 : if( nItemCount == 7 )
634 0 : nItemsPerLine = 7;
635 2 : else if( nItemCount == 8 )
636 0 : nItemsPerLine = 8;
637 2 : else if( nItemCount < 14 )
638 : {
639 0 : CPLDebug("ERS", "Invalid item count for ControlPoints");
640 0 : CSLDestroy( papszTokens );
641 0 : return;
642 : }
643 2 : else if( EQUAL(papszTokens[8],"Yes") || EQUAL(papszTokens[8],"No") )
644 0 : nItemsPerLine = 7;
645 4 : else if( EQUAL(papszTokens[9],"Yes") || EQUAL(papszTokens[9],"No") )
646 2 : nItemsPerLine = 8;
647 : else
648 : {
649 0 : CPLDebug("ERS", "Invalid format for ControlPoints");
650 0 : CSLDestroy( papszTokens );
651 0 : return;
652 : }
653 :
654 : /* -------------------------------------------------------------------- */
655 : /* Setup GCPs. */
656 : /* -------------------------------------------------------------------- */
657 : int iGCP;
658 :
659 2 : CPLAssert( nGCPCount == 0 );
660 :
661 2 : nGCPCount = nItemCount / nItemsPerLine;
662 2 : pasGCPList = (GDAL_GCP *) CPLCalloc(nGCPCount,sizeof(GDAL_GCP));
663 2 : GDALInitGCPs( nGCPCount, pasGCPList );
664 :
665 10 : for( iGCP = 0; iGCP < nGCPCount; iGCP++ )
666 : {
667 8 : GDAL_GCP *psGCP = pasGCPList + iGCP;
668 :
669 8 : CPLFree( psGCP->pszId );
670 8 : psGCP->pszId = CPLStrdup(papszTokens[iGCP*nItemsPerLine+0]);
671 8 : psGCP->dfGCPPixel = atof(papszTokens[iGCP*nItemsPerLine+3]);
672 8 : psGCP->dfGCPLine = atof(papszTokens[iGCP*nItemsPerLine+4]);
673 8 : psGCP->dfGCPX = atof(papszTokens[iGCP*nItemsPerLine+5]);
674 8 : psGCP->dfGCPY = atof(papszTokens[iGCP*nItemsPerLine+6]);
675 8 : if( nItemsPerLine == 8 )
676 8 : psGCP->dfGCPZ = atof(papszTokens[iGCP*nItemsPerLine+7]);
677 : }
678 :
679 2 : CSLDestroy( papszTokens );
680 :
681 : /* -------------------------------------------------------------------- */
682 : /* Parse the GCP projection. */
683 : /* -------------------------------------------------------------------- */
684 2 : OGRSpatialReference oSRS;
685 :
686 2 : osProj = poHeader->Find( "RasterInfo.WarpControl.CoordinateSpace.Projection", "" );
687 2 : osDatum = poHeader->Find( "RasterInfo.WarpControl.CoordinateSpace.Datum", "" );
688 2 : osUnits = poHeader->Find( "RasterInfo.WarpControl.CoordinateSpace.Units", "" );
689 :
690 : oSRS.importFromERM( osProj.size() ? osProj.c_str() : "RAW",
691 : osDatum.size() ? osDatum.c_str() : "WGS84",
692 2 : osUnits.size() ? osUnits.c_str() : "METERS" );
693 :
694 2 : CPLFree( pszGCPProjection );
695 2 : oSRS.exportToWkt( &pszGCPProjection );
696 : }
697 :
698 : /************************************************************************/
699 : /* ==================================================================== */
700 : /* ERSRasterBand */
701 : /* ==================================================================== */
702 : /************************************************************************/
703 :
704 : class ERSRasterBand : public RawRasterBand
705 90 : {
706 : public:
707 :
708 : ERSRasterBand( GDALDataset *poDS, int nBand, void * fpRaw,
709 : vsi_l_offset nImgOffset, int nPixelOffset,
710 : int nLineOffset,
711 : GDALDataType eDataType, int bNativeOrder,
712 : int bIsVSIL = FALSE, int bOwnsFP = FALSE );
713 :
714 : virtual double GetNoDataValue( int *pbSuccess = NULL );
715 : virtual CPLErr SetNoDataValue( double );
716 : };
717 :
718 : /************************************************************************/
719 : /* ERSRasterBand() */
720 : /************************************************************************/
721 :
722 90 : ERSRasterBand::ERSRasterBand( GDALDataset *poDS, int nBand, void * fpRaw,
723 : vsi_l_offset nImgOffset, int nPixelOffset,
724 : int nLineOffset,
725 : GDALDataType eDataType, int bNativeOrder,
726 : int bIsVSIL, int bOwnsFP ) :
727 : RawRasterBand(poDS, nBand, fpRaw, nImgOffset, nPixelOffset,
728 90 : nLineOffset, eDataType, bNativeOrder, bIsVSIL, bOwnsFP)
729 : {
730 90 : }
731 :
732 : /************************************************************************/
733 : /* GetNoDataValue() */
734 : /************************************************************************/
735 :
736 11 : double ERSRasterBand::GetNoDataValue( int *pbSuccess )
737 : {
738 11 : ERSDataset* poGDS = (ERSDataset*) poDS;
739 11 : if (poGDS->bHasNoDataValue)
740 : {
741 1 : if (pbSuccess)
742 1 : *pbSuccess = TRUE;
743 1 : return poGDS->dfNoDataValue;
744 : }
745 :
746 10 : return RawRasterBand::GetNoDataValue(pbSuccess);
747 : }
748 :
749 : /************************************************************************/
750 : /* SetNoDataValue() */
751 : /************************************************************************/
752 :
753 1 : CPLErr ERSRasterBand::SetNoDataValue( double dfNoDataValue )
754 : {
755 1 : ERSDataset* poGDS = (ERSDataset*) poDS;
756 1 : if (!poGDS->bHasNoDataValue || poGDS->dfNoDataValue != dfNoDataValue)
757 : {
758 1 : poGDS->bHasNoDataValue = TRUE;
759 1 : poGDS->dfNoDataValue = dfNoDataValue;
760 :
761 1 : poGDS->bHDRDirty = TRUE;
762 : poGDS->poHeader->Set( "RasterInfo.NullCellValue",
763 1 : CPLString().Printf( "%.16g", dfNoDataValue) );
764 : }
765 1 : return CE_None;
766 : }
767 :
768 : /************************************************************************/
769 : /* Identify() */
770 : /************************************************************************/
771 :
772 12752 : int ERSDataset::Identify( GDALOpenInfo * poOpenInfo )
773 :
774 : {
775 : /* -------------------------------------------------------------------- */
776 : /* We assume the user selects the .ers file. */
777 : /* -------------------------------------------------------------------- */
778 12752 : if( poOpenInfo->nHeaderBytes > 15
779 : && EQUALN((const char *) poOpenInfo->pabyHeader,"Algorithm Begin",15) )
780 : {
781 : CPLError( CE_Failure, CPLE_OpenFailed,
782 : "%s appears to be an algorithm ERS file, which is not currently supported.",
783 0 : poOpenInfo->pszFilename );
784 0 : return FALSE;
785 : }
786 :
787 : /* -------------------------------------------------------------------- */
788 : /* We assume the user selects the .ers file. */
789 : /* -------------------------------------------------------------------- */
790 12752 : if( poOpenInfo->nHeaderBytes < 15
791 : || !EQUALN((const char *) poOpenInfo->pabyHeader,"DatasetHeader ",14) )
792 12702 : return FALSE;
793 :
794 50 : return TRUE;
795 : }
796 :
797 : /************************************************************************/
798 : /* Open() */
799 : /************************************************************************/
800 :
801 2659 : GDALDataset *ERSDataset::Open( GDALOpenInfo * poOpenInfo )
802 :
803 : {
804 2659 : if( !Identify( poOpenInfo ) )
805 2609 : return NULL;
806 :
807 : /* -------------------------------------------------------------------- */
808 : /* Open the .ers file, and read the first line. */
809 : /* -------------------------------------------------------------------- */
810 50 : VSILFILE *fpERS = VSIFOpenL( poOpenInfo->pszFilename, "rb" );
811 :
812 50 : if( fpERS == NULL )
813 0 : return NULL;
814 :
815 50 : CPLReadLineL( fpERS );
816 :
817 : /* -------------------------------------------------------------------- */
818 : /* Now ingest the rest of the file as a tree of header nodes. */
819 : /* -------------------------------------------------------------------- */
820 50 : ERSHdrNode *poHeader = new ERSHdrNode();
821 :
822 50 : if( !poHeader->ParseChildren( fpERS ) )
823 : {
824 0 : delete poHeader;
825 0 : VSIFCloseL( fpERS );
826 0 : return NULL;
827 : }
828 :
829 50 : VSIFCloseL( fpERS );
830 :
831 : /* -------------------------------------------------------------------- */
832 : /* Do we have the minimum required information from this header? */
833 : /* -------------------------------------------------------------------- */
834 50 : if( poHeader->Find( "RasterInfo.NrOfLines" ) == NULL
835 : || poHeader->Find( "RasterInfo.NrOfCellsPerLine" ) == NULL
836 : || poHeader->Find( "RasterInfo.NrOfBands" ) == NULL )
837 : {
838 0 : if( poHeader->FindNode( "Algorithm" ) != NULL )
839 : {
840 : CPLError( CE_Failure, CPLE_OpenFailed,
841 : "%s appears to be an algorithm ERS file, which is not currently supported.",
842 0 : poOpenInfo->pszFilename );
843 : }
844 0 : delete poHeader;
845 0 : return NULL;
846 : }
847 :
848 : /* -------------------------------------------------------------------- */
849 : /* Create a corresponding GDALDataset. */
850 : /* -------------------------------------------------------------------- */
851 : ERSDataset *poDS;
852 :
853 50 : poDS = new ERSDataset();
854 50 : poDS->poHeader = poHeader;
855 50 : poDS->eAccess = poOpenInfo->eAccess;
856 :
857 : /* -------------------------------------------------------------------- */
858 : /* Capture some information from the file that is of interest. */
859 : /* -------------------------------------------------------------------- */
860 50 : int nBands = atoi(poHeader->Find( "RasterInfo.NrOfBands" ));
861 50 : poDS->nRasterXSize = atoi(poHeader->Find( "RasterInfo.NrOfCellsPerLine" ));
862 50 : poDS->nRasterYSize = atoi(poHeader->Find( "RasterInfo.NrOfLines" ));
863 :
864 100 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize) ||
865 : !GDALCheckBandCount(nBands, FALSE))
866 : {
867 0 : delete poDS;
868 0 : return NULL;
869 : }
870 :
871 : /* -------------------------------------------------------------------- */
872 : /* Get the HeaderOffset if it exists in the header */
873 : /* -------------------------------------------------------------------- */
874 50 : GIntBig nHeaderOffset = 0;
875 50 : if( poHeader->Find( "HeaderOffset" ) != NULL )
876 : {
877 1 : nHeaderOffset = atoi(poHeader->Find( "HeaderOffset" ));
878 : }
879 :
880 : /* -------------------------------------------------------------------- */
881 : /* Establish the data type. */
882 : /* -------------------------------------------------------------------- */
883 : GDALDataType eType;
884 : CPLString osCellType = poHeader->Find( "RasterInfo.CellType",
885 50 : "Unsigned8BitInteger" );
886 50 : if( EQUAL(osCellType,"Unsigned8BitInteger") )
887 22 : eType = GDT_Byte;
888 28 : else if( EQUAL(osCellType,"Signed8BitInteger") )
889 4 : eType = GDT_Byte;
890 24 : else if( EQUAL(osCellType,"Unsigned16BitInteger") )
891 3 : eType = GDT_UInt16;
892 21 : else if( EQUAL(osCellType,"Signed16BitInteger") )
893 4 : eType = GDT_Int16;
894 17 : else if( EQUAL(osCellType,"Unsigned32BitInteger") )
895 3 : eType = GDT_UInt32;
896 14 : else if( EQUAL(osCellType,"Signed32BitInteger") )
897 3 : eType = GDT_Int32;
898 11 : else if( EQUAL(osCellType,"IEEE4ByteReal") )
899 8 : eType = GDT_Float32;
900 3 : else if( EQUAL(osCellType,"IEEE8ByteReal") )
901 3 : eType = GDT_Float64;
902 : else
903 : {
904 0 : CPLDebug( "ERS", "Unknown CellType '%s'", osCellType.c_str() );
905 0 : eType = GDT_Byte;
906 : }
907 :
908 : /* -------------------------------------------------------------------- */
909 : /* Pick up the word order. */
910 : /* -------------------------------------------------------------------- */
911 : int bNative;
912 :
913 : #ifdef CPL_LSB
914 : bNative = EQUAL(poHeader->Find( "ByteOrder", "LSBFirst" ),
915 50 : "LSBFirst");
916 : #else
917 : bNative = EQUAL(poHeader->Find( "ByteOrder", "MSBFirst" ),
918 : "MSBFirst");
919 : #endif
920 :
921 : /* -------------------------------------------------------------------- */
922 : /* Figure out the name of the target file. */
923 : /* -------------------------------------------------------------------- */
924 50 : CPLString osPath = CPLGetPath( poOpenInfo->pszFilename );
925 50 : CPLString osDataFile = poHeader->Find( "DataFile", "" );
926 50 : CPLString osDataFilePath;
927 :
928 50 : if( osDataFile.length() == 0 ) // just strip off extension.
929 : {
930 50 : osDataFile = CPLGetFilename( poOpenInfo->pszFilename );
931 50 : osDataFile = osDataFile.substr( 0, osDataFile.find_last_of('.') );
932 : }
933 :
934 50 : osDataFilePath = CPLFormFilename( osPath, osDataFile, NULL );
935 :
936 : /* -------------------------------------------------------------------- */
937 : /* DataSetType = Translated files are links to things like ecw */
938 : /* files. */
939 : /* -------------------------------------------------------------------- */
940 50 : if( EQUAL(poHeader->Find("DataSetType",""),"Translated") )
941 : {
942 : poDS->poDepFile = (GDALDataset *)
943 0 : GDALOpenShared( osDataFilePath, poOpenInfo->eAccess );
944 :
945 0 : if( poDS->poDepFile != NULL
946 : && poDS->poDepFile->GetRasterCount() >= nBands )
947 : {
948 : int iBand;
949 :
950 0 : for( iBand = 0; iBand < nBands; iBand++ )
951 : {
952 : // Assume pixel interleaved.
953 : poDS->SetBand( iBand+1,
954 0 : poDS->poDepFile->GetRasterBand( iBand+1 ) );
955 : }
956 : }
957 : }
958 :
959 : /* ==================================================================== */
960 : /* While ERStorage indicates a raw file. */
961 : /* ==================================================================== */
962 50 : else if( EQUAL(poHeader->Find("DataSetType",""),"ERStorage") )
963 : {
964 : // Open data file.
965 50 : if( poOpenInfo->eAccess == GA_Update )
966 34 : poDS->fpImage = VSIFOpenL( osDataFilePath, "r+" );
967 : else
968 16 : poDS->fpImage = VSIFOpenL( osDataFilePath, "r" );
969 :
970 50 : poDS->osRawFilename = osDataFilePath;
971 :
972 50 : if( poDS->fpImage != NULL )
973 : {
974 50 : int iWordSize = GDALGetDataTypeSize(eType) / 8;
975 : int iBand;
976 :
977 140 : for( iBand = 0; iBand < nBands; iBand++ )
978 : {
979 : // Assume pixel interleaved.
980 : poDS->SetBand(
981 : iBand+1,
982 : new ERSRasterBand( poDS, iBand+1, poDS->fpImage,
983 : nHeaderOffset
984 : + iWordSize * iBand * poDS->nRasterXSize,
985 : iWordSize,
986 : iWordSize * nBands * poDS->nRasterXSize,
987 90 : eType, bNative, TRUE ));
988 90 : if( EQUAL(osCellType,"Signed8BitInteger") )
989 : poDS->GetRasterBand(iBand+1)->
990 : SetMetadataItem( "PIXELTYPE", "SIGNEDBYTE",
991 12 : "IMAGE_STRUCTURE" );
992 : }
993 : }
994 : }
995 :
996 : /* -------------------------------------------------------------------- */
997 : /* Otherwise we have an error! */
998 : /* -------------------------------------------------------------------- */
999 50 : if( poDS->nBands == 0 )
1000 : {
1001 0 : delete poDS;
1002 0 : return NULL;
1003 : }
1004 :
1005 : /* -------------------------------------------------------------------- */
1006 : /* Look for band descriptions. */
1007 : /* -------------------------------------------------------------------- */
1008 50 : int iChild, iBand = 0;
1009 50 : ERSHdrNode *poRI = poHeader->FindNode( "RasterInfo" );
1010 :
1011 287 : for( iChild = 0;
1012 : poRI != NULL && iChild < poRI->nItemCount && iBand < poDS->nBands;
1013 : iChild++ )
1014 : {
1015 263 : if( poRI->papoItemChild[iChild] != NULL
1016 26 : && EQUAL(poRI->papszItemName[iChild],"BandId") )
1017 : {
1018 : const char *pszValue =
1019 9 : poRI->papoItemChild[iChild]->Find( "Value", NULL );
1020 :
1021 9 : iBand++;
1022 9 : if( pszValue )
1023 : {
1024 9 : CPLPushErrorHandler( CPLQuietErrorHandler );
1025 9 : poDS->GetRasterBand( iBand )->SetDescription( pszValue );
1026 9 : CPLPopErrorHandler();
1027 : }
1028 :
1029 9 : pszValue = poRI->papoItemChild[iChild]->Find( "Units", NULL );
1030 9 : if ( pszValue )
1031 : {
1032 2 : CPLPushErrorHandler( CPLQuietErrorHandler );
1033 2 : poDS->GetRasterBand( iBand )->SetUnitType( pszValue );
1034 2 : CPLPopErrorHandler();
1035 : }
1036 : }
1037 : }
1038 :
1039 : /* -------------------------------------------------------------------- */
1040 : /* Look for projection. */
1041 : /* -------------------------------------------------------------------- */
1042 50 : OGRSpatialReference oSRS;
1043 :
1044 50 : poDS->osProj = poHeader->Find( "CoordinateSpace.Projection", "" );
1045 50 : poDS->osDatum = poHeader->Find( "CoordinateSpace.Datum", "" );
1046 50 : poDS->osUnits = poHeader->Find( "CoordinateSpace.Units", "" );
1047 :
1048 : oSRS.importFromERM( poDS->osProj.size() ? poDS->osProj.c_str() : "RAW",
1049 : poDS->osDatum.size() ? poDS->osDatum.c_str() : "WGS84",
1050 50 : poDS->osUnits.size() ? poDS->osUnits.c_str() : "METERS" );
1051 :
1052 50 : CPLFree( poDS->pszProjection );
1053 50 : oSRS.exportToWkt( &(poDS->pszProjection) );
1054 :
1055 : /* -------------------------------------------------------------------- */
1056 : /* Look for the geotransform. */
1057 : /* -------------------------------------------------------------------- */
1058 50 : if( poHeader->Find( "RasterInfo.RegistrationCoord.Eastings", NULL ) )
1059 : {
1060 6 : poDS->bGotTransform = TRUE;
1061 : poDS->adfGeoTransform[0] = CPLAtof(
1062 6 : poHeader->Find( "RasterInfo.RegistrationCoord.Eastings", "" ));
1063 : poDS->adfGeoTransform[1] = CPLAtof(
1064 6 : poHeader->Find( "RasterInfo.CellInfo.Xdimension", "1.0" ));
1065 6 : poDS->adfGeoTransform[2] = 0.0;
1066 : poDS->adfGeoTransform[3] = CPLAtof(
1067 6 : poHeader->Find( "RasterInfo.RegistrationCoord.Northings", "" ));
1068 6 : poDS->adfGeoTransform[4] = 0.0;
1069 : poDS->adfGeoTransform[5] = -CPLAtof(
1070 6 : poHeader->Find( "RasterInfo.CellInfo.Ydimension", "1.0" ));
1071 : }
1072 44 : else if( poHeader->Find( "RasterInfo.RegistrationCoord.Latitude", NULL )
1073 : && poHeader->Find( "RasterInfo.CellInfo.Xdimension", NULL ) )
1074 : {
1075 2 : poDS->bGotTransform = TRUE;
1076 : poDS->adfGeoTransform[0] = ERSDMS2Dec(
1077 2 : poHeader->Find( "RasterInfo.RegistrationCoord.Longitude", "" ));
1078 : poDS->adfGeoTransform[1] = CPLAtof(
1079 2 : poHeader->Find( "RasterInfo.CellInfo.Xdimension", "" ));
1080 2 : poDS->adfGeoTransform[2] = 0.0;
1081 : poDS->adfGeoTransform[3] = ERSDMS2Dec(
1082 2 : poHeader->Find( "RasterInfo.RegistrationCoord.Latitude", "" ));
1083 2 : poDS->adfGeoTransform[4] = 0.0;
1084 : poDS->adfGeoTransform[5] = -CPLAtof(
1085 2 : poHeader->Find( "RasterInfo.CellInfo.Ydimension", "" ));
1086 : }
1087 :
1088 : /* -------------------------------------------------------------------- */
1089 : /* Adjust if we have a registration cell. */
1090 : /* -------------------------------------------------------------------- */
1091 50 : int iCellX = atoi(poHeader->Find("RasterInfo.RegistrationCellX", "1"));
1092 50 : int iCellY = atoi(poHeader->Find("RasterInfo.RegistrationCellY", "1"));
1093 :
1094 50 : if( poDS->bGotTransform )
1095 : {
1096 : poDS->adfGeoTransform[0] -=
1097 8 : (iCellX-1) * poDS->adfGeoTransform[1]
1098 16 : + (iCellY-1) * poDS->adfGeoTransform[2];
1099 : poDS->adfGeoTransform[3] -=
1100 8 : (iCellX-1) * poDS->adfGeoTransform[4]
1101 16 : + (iCellY-1) * poDS->adfGeoTransform[5];
1102 : }
1103 :
1104 : /* -------------------------------------------------------------------- */
1105 : /* Check for null values. */
1106 : /* -------------------------------------------------------------------- */
1107 50 : if( poHeader->Find( "RasterInfo.NullCellValue", NULL ) )
1108 : {
1109 5 : poDS->bHasNoDataValue = TRUE;
1110 5 : poDS->dfNoDataValue = CPLAtofM(poHeader->Find( "RasterInfo.NullCellValue" ));
1111 :
1112 5 : if (poDS->poDepFile != NULL)
1113 : {
1114 0 : CPLPushErrorHandler( CPLQuietErrorHandler );
1115 :
1116 0 : for( iBand = 1; iBand <= poDS->nBands; iBand++ )
1117 0 : poDS->GetRasterBand(iBand)->SetNoDataValue(poDS->dfNoDataValue);
1118 :
1119 0 : CPLPopErrorHandler();
1120 : }
1121 : }
1122 :
1123 : /* -------------------------------------------------------------------- */
1124 : /* Do we have an "All" region? */
1125 : /* -------------------------------------------------------------------- */
1126 50 : ERSHdrNode *poAll = NULL;
1127 :
1128 289 : for( iChild = 0;
1129 : poRI != NULL && iChild < poRI->nItemCount;
1130 : iChild++ )
1131 : {
1132 266 : if( poRI->papoItemChild[iChild] != NULL
1133 27 : && EQUAL(poRI->papszItemName[iChild],"RegionInfo") )
1134 : {
1135 0 : if( EQUAL(poRI->papoItemChild[iChild]->Find("RegionName",""),
1136 : "All") )
1137 0 : poAll = poRI->papoItemChild[iChild];
1138 : }
1139 : }
1140 :
1141 : /* -------------------------------------------------------------------- */
1142 : /* Do we have statistics? */
1143 : /* -------------------------------------------------------------------- */
1144 50 : if( poAll && poAll->FindNode( "Stats" ) )
1145 : {
1146 0 : CPLPushErrorHandler( CPLQuietErrorHandler );
1147 :
1148 0 : for( iBand = 1; iBand <= poDS->nBands; iBand++ )
1149 : {
1150 : const char *pszValue =
1151 0 : poAll->FindElem( "Stats.MinimumValue", iBand-1 );
1152 :
1153 0 : if( pszValue )
1154 : poDS->GetRasterBand(iBand)->SetMetadataItem(
1155 0 : "STATISTICS_MINIMUM", pszValue );
1156 :
1157 0 : pszValue = poAll->FindElem( "Stats.MaximumValue", iBand-1 );
1158 :
1159 0 : if( pszValue )
1160 : poDS->GetRasterBand(iBand)->SetMetadataItem(
1161 0 : "STATISTICS_MAXIMUM", pszValue );
1162 :
1163 0 : pszValue = poAll->FindElem( "Stats.MeanValue", iBand-1 );
1164 :
1165 0 : if( pszValue )
1166 : poDS->GetRasterBand(iBand)->SetMetadataItem(
1167 0 : "STATISTICS_MEAN", pszValue );
1168 :
1169 0 : pszValue = poAll->FindElem( "Stats.MedianValue", iBand-1 );
1170 :
1171 0 : if( pszValue )
1172 : poDS->GetRasterBand(iBand)->SetMetadataItem(
1173 0 : "STATISTICS_MEDIAN", pszValue );
1174 : }
1175 :
1176 0 : CPLPopErrorHandler();
1177 :
1178 : }
1179 :
1180 : /* -------------------------------------------------------------------- */
1181 : /* Do we have GCPs. */
1182 : /* -------------------------------------------------------------------- */
1183 50 : if( poHeader->FindNode( "RasterInfo.WarpControl" ) )
1184 2 : poDS->ReadGCPs();
1185 :
1186 : /* -------------------------------------------------------------------- */
1187 : /* Initialize any PAM information. */
1188 : /* -------------------------------------------------------------------- */
1189 50 : poDS->SetDescription( poOpenInfo->pszFilename );
1190 50 : poDS->TryLoadXML();
1191 :
1192 : // if no SR in xml, try aux
1193 50 : const char* pszPrj = poDS->GDALPamDataset::GetProjectionRef();
1194 50 : if( !pszPrj || strlen(pszPrj) == 0 )
1195 : {
1196 : // try aux
1197 50 : GDALDataset* poAuxDS = GDALFindAssociatedAuxFile( poOpenInfo->pszFilename, GA_ReadOnly, poDS );
1198 50 : if( poAuxDS )
1199 : {
1200 0 : pszPrj = poAuxDS->GetProjectionRef();
1201 0 : if( pszPrj && strlen(pszPrj) > 0 )
1202 : {
1203 0 : CPLFree( poDS->pszProjection );
1204 0 : poDS->pszProjection = CPLStrdup(pszPrj);
1205 : }
1206 :
1207 0 : GDALClose( poAuxDS );
1208 : }
1209 : }
1210 : /* -------------------------------------------------------------------- */
1211 : /* Check for overviews. */
1212 : /* -------------------------------------------------------------------- */
1213 50 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
1214 :
1215 50 : return( poDS );
1216 : }
1217 :
1218 : /************************************************************************/
1219 : /* Create() */
1220 : /************************************************************************/
1221 :
1222 60 : GDALDataset *ERSDataset::Create( const char * pszFilename,
1223 : int nXSize, int nYSize, int nBands,
1224 : GDALDataType eType, char ** papszOptions )
1225 :
1226 : {
1227 : /* -------------------------------------------------------------------- */
1228 : /* Verify settings. */
1229 : /* -------------------------------------------------------------------- */
1230 60 : if (nBands <= 0)
1231 : {
1232 : CPLError( CE_Failure, CPLE_NotSupported,
1233 1 : "ERS driver does not support %d bands.\n", nBands);
1234 1 : return NULL;
1235 : }
1236 :
1237 59 : if( eType != GDT_Byte && eType != GDT_Int16 && eType != GDT_UInt16
1238 : && eType != GDT_Int32 && eType != GDT_UInt32
1239 : && eType != GDT_Float32 && eType != GDT_Float64 )
1240 : {
1241 : CPLError( CE_Failure, CPLE_AppDefined,
1242 : "The ERS driver does not supporting creating files of types %s.",
1243 12 : GDALGetDataTypeName( eType ) );
1244 12 : return NULL;
1245 : }
1246 :
1247 : /* -------------------------------------------------------------------- */
1248 : /* Work out the name we want to use for the .ers and binary */
1249 : /* data files. */
1250 : /* -------------------------------------------------------------------- */
1251 47 : CPLString osBinFile, osErsFile;
1252 :
1253 47 : if( EQUAL(CPLGetExtension( pszFilename ), "ers") )
1254 : {
1255 6 : osErsFile = pszFilename;
1256 6 : osBinFile = osErsFile.substr(0,osErsFile.length()-4);
1257 : }
1258 : else
1259 : {
1260 41 : osBinFile = pszFilename;
1261 41 : osErsFile = osBinFile + ".ers";
1262 : }
1263 :
1264 : /* -------------------------------------------------------------------- */
1265 : /* Work out some values we will write. */
1266 : /* -------------------------------------------------------------------- */
1267 47 : const char *pszCellType = "Unsigned8BitInteger";
1268 :
1269 47 : if( eType == GDT_Byte )
1270 22 : pszCellType = "Unsigned8BitInteger";
1271 25 : else if( eType == GDT_Int16 )
1272 4 : pszCellType = "Signed16BitInteger";
1273 21 : else if( eType == GDT_UInt16 )
1274 4 : pszCellType = "Unsigned16BitInteger";
1275 17 : else if( eType == GDT_Int32 )
1276 4 : pszCellType = "Signed32BitInteger";
1277 13 : else if( eType == GDT_UInt32 )
1278 4 : pszCellType = "Unsigned32BitInteger";
1279 9 : else if( eType == GDT_Float32 )
1280 5 : pszCellType = "IEEE4ByteReal";
1281 4 : else if( eType == GDT_Float64 )
1282 4 : pszCellType = "IEEE8ByteReal";
1283 : else
1284 : {
1285 0 : CPLAssert( FALSE );
1286 : }
1287 :
1288 : /* -------------------------------------------------------------------- */
1289 : /* Handling for signed eight bit data. */
1290 : /* -------------------------------------------------------------------- */
1291 47 : const char *pszPixelType = CSLFetchNameValue( papszOptions, "PIXELTYPE" );
1292 47 : if( pszPixelType
1293 : && EQUAL(pszPixelType,"SIGNEDBYTE")
1294 : && eType == GDT_Byte )
1295 1 : pszCellType = "Signed8BitInteger";
1296 :
1297 : /* -------------------------------------------------------------------- */
1298 : /* Write binary file. */
1299 : /* -------------------------------------------------------------------- */
1300 : GUIntBig nSize;
1301 47 : GByte byZero = 0;
1302 :
1303 47 : VSILFILE *fpBin = VSIFOpenL( osBinFile, "w" );
1304 :
1305 47 : if( fpBin == NULL )
1306 : {
1307 : CPLError( CE_Failure, CPLE_FileIO,
1308 : "Failed to create %s:\n%s",
1309 13 : osBinFile.c_str(), VSIStrerror( errno ) );
1310 13 : return NULL;
1311 : }
1312 :
1313 : nSize = nXSize * (GUIntBig) nYSize
1314 34 : * nBands * (GDALGetDataTypeSize(eType) / 8);
1315 34 : if( VSIFSeekL( fpBin, nSize-1, SEEK_SET ) != 0
1316 : || VSIFWriteL( &byZero, 1, 1, fpBin ) != 1 )
1317 : {
1318 : CPLError( CE_Failure, CPLE_FileIO,
1319 : "Failed to write %s:\n%s",
1320 0 : osBinFile.c_str(), VSIStrerror( errno ) );
1321 0 : VSIFCloseL( fpBin );
1322 0 : return NULL;
1323 : }
1324 34 : VSIFCloseL( fpBin );
1325 :
1326 :
1327 : /* -------------------------------------------------------------------- */
1328 : /* Try writing header file. */
1329 : /* -------------------------------------------------------------------- */
1330 34 : VSILFILE *fpERS = VSIFOpenL( osErsFile, "w" );
1331 :
1332 34 : if( fpERS == NULL )
1333 : {
1334 : CPLError( CE_Failure, CPLE_FileIO,
1335 : "Failed to create %s:\n%s",
1336 0 : osErsFile.c_str(), VSIStrerror( errno ) );
1337 0 : return NULL;
1338 : }
1339 :
1340 34 : VSIFPrintfL( fpERS, "DatasetHeader Begin\n" );
1341 34 : VSIFPrintfL( fpERS, "\tVersion\t\t = \"6.0\"\n" );
1342 34 : VSIFPrintfL( fpERS, "\tName\t\t= \"%s\"\n", CPLGetFilename(osErsFile) );
1343 :
1344 : // Last updated requires timezone info which we don't necessarily get
1345 : // get from VSICTime() so perhaps it is better to omit this.
1346 : // VSIFPrintfL( fpERS, "\tLastUpdated\t= %s",
1347 : // VSICTime( VSITime( NULL ) ) );
1348 :
1349 34 : VSIFPrintfL( fpERS, "\tDataSetType\t= ERStorage\n" );
1350 34 : VSIFPrintfL( fpERS, "\tDataType\t= Raster\n" );
1351 34 : VSIFPrintfL( fpERS, "\tByteOrder\t= LSBFirst\n" );
1352 34 : VSIFPrintfL( fpERS, "\tRasterInfo Begin\n" );
1353 34 : VSIFPrintfL( fpERS, "\t\tCellType\t= %s\n", pszCellType );
1354 34 : VSIFPrintfL( fpERS, "\t\tNrOfLines\t= %d\n", nYSize );
1355 34 : VSIFPrintfL( fpERS, "\t\tNrOfCellsPerLine\t= %d\n", nXSize );
1356 34 : VSIFPrintfL( fpERS, "\t\tNrOfBands\t= %d\n", nBands );
1357 34 : VSIFPrintfL( fpERS, "\tRasterInfo End\n" );
1358 34 : if( VSIFPrintfL( fpERS, "DatasetHeader End\n" ) < 17 )
1359 : {
1360 : CPLError( CE_Failure, CPLE_FileIO,
1361 : "Failed to write %s:\n%s",
1362 0 : osErsFile.c_str(), VSIStrerror( errno ) );
1363 0 : return NULL;
1364 : }
1365 :
1366 34 : VSIFCloseL( fpERS );
1367 :
1368 : /* -------------------------------------------------------------------- */
1369 : /* Reopen. */
1370 : /* -------------------------------------------------------------------- */
1371 34 : GDALOpenInfo oOpenInfo( osErsFile, GA_Update );
1372 34 : ERSDataset* poDS = (ERSDataset*) Open( &oOpenInfo );
1373 34 : if (poDS == NULL)
1374 0 : return NULL;
1375 :
1376 : /* -------------------------------------------------------------------- */
1377 : /* Fetch DATUM, PROJ and UNITS creation option */
1378 : /* -------------------------------------------------------------------- */
1379 34 : const char *pszDatum = CSLFetchNameValue( papszOptions, "DATUM" );
1380 34 : if (pszDatum)
1381 1 : poDS->osDatum = pszDatum;
1382 34 : const char *pszProj = CSLFetchNameValue( papszOptions, "PROJ" );
1383 34 : if (pszProj)
1384 1 : poDS->osProj = pszProj;
1385 34 : const char *pszUnits = CSLFetchNameValue( papszOptions, "UNITS" );
1386 34 : if (pszUnits)
1387 1 : poDS->osUnits = pszUnits;
1388 :
1389 34 : if (pszDatum || pszProj || pszUnits)
1390 : {
1391 : poDS->WriteProjectionInfo(pszProj ? pszProj : "RAW",
1392 : pszDatum ? pszDatum : "RAW",
1393 1 : pszUnits ? pszUnits : "METERS");
1394 : }
1395 :
1396 34 : return poDS;
1397 : }
1398 :
1399 : /************************************************************************/
1400 : /* GDALRegister_ERS() */
1401 : /************************************************************************/
1402 :
1403 582 : void GDALRegister_ERS()
1404 :
1405 : {
1406 : GDALDriver *poDriver;
1407 :
1408 582 : if( GDALGetDriverByName( "ERS" ) == NULL )
1409 : {
1410 561 : poDriver = new GDALDriver();
1411 :
1412 561 : poDriver->SetDescription( "ERS" );
1413 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
1414 561 : "ERMapper .ers Labelled" );
1415 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
1416 561 : "frmt_ers.html" );
1417 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONDATATYPES,
1418 561 : "Byte Int16 UInt16 Int32 UInt32 Float32 Float64" );
1419 :
1420 : poDriver->SetMetadataItem( GDAL_DMD_CREATIONOPTIONLIST,
1421 : "<CreationOptionList>"
1422 : " <Option name='PIXELTYPE' type='string' description='By setting this to SIGNEDBYTE, a new Byte file can be forced to be written as signed byte'/>"
1423 : " <Option name='PROJ' type='string' description='ERS Projection Name'/>"
1424 : " <Option name='DATUM' type='string' description='ERS Datum Name' />"
1425 : " <Option name='UNITS' type='string-select' description='ERS Projection Units'>"
1426 : " <Value>METERS</Value>"
1427 : " <Value>FEET</Value>"
1428 : " </Option>"
1429 561 : "</CreationOptionList>" );
1430 :
1431 561 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
1432 :
1433 561 : poDriver->pfnOpen = ERSDataset::Open;
1434 561 : poDriver->pfnIdentify = ERSDataset::Identify;
1435 561 : poDriver->pfnCreate = ERSDataset::Create;
1436 :
1437 561 : GetGDALDriverManager()->RegisterDriver( poDriver );
1438 : }
1439 582 : }
|
