1 : /******************************************************************************
2 : * $Id: ndfdataset.cpp 21749 2011-02-19 12:25:29Z rouault $
3 : *
4 : * Project: NDF Driver
5 : * Purpose: Implementation of NLAPS Data Format read support.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2005, Frank Warmerdam
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 :
34 : CPL_CVSID("$Id: ndfdataset.cpp 21749 2011-02-19 12:25:29Z rouault $");
35 :
36 : /************************************************************************/
37 : /* ==================================================================== */
38 : /* NDFDataset */
39 : /* ==================================================================== */
40 : /************************************************************************/
41 :
42 : class NDFDataset : public RawDataset
43 : {
44 : double adfGeoTransform[6];
45 :
46 : char *pszProjection;
47 : char **papszExtraFiles;
48 :
49 : char **papszHeader;
50 : const char *Get( const char *pszKey, const char *pszDefault);
51 :
52 : public:
53 : NDFDataset();
54 : ~NDFDataset();
55 :
56 : virtual CPLErr GetGeoTransform( double * padfTransform );
57 : virtual const char *GetProjectionRef(void);
58 : virtual char **GetFileList(void);
59 :
60 : static GDALDataset *Open( GDALOpenInfo * );
61 : };
62 :
63 : /************************************************************************/
64 : /* NDFDataset() */
65 : /************************************************************************/
66 :
67 1 : NDFDataset::NDFDataset()
68 : {
69 1 : pszProjection = CPLStrdup("");
70 :
71 1 : papszHeader = NULL;
72 1 : papszExtraFiles = NULL;
73 :
74 1 : adfGeoTransform[0] = 0.0;
75 1 : adfGeoTransform[1] = 1.0;
76 1 : adfGeoTransform[2] = 0.0;
77 1 : adfGeoTransform[3] = 0.0;
78 1 : adfGeoTransform[4] = 0.0;
79 1 : adfGeoTransform[5] = 1.0;
80 1 : }
81 :
82 : /************************************************************************/
83 : /* ~NDFDataset() */
84 : /************************************************************************/
85 :
86 1 : NDFDataset::~NDFDataset()
87 :
88 : {
89 1 : FlushCache();
90 1 : CPLFree( pszProjection );
91 1 : CSLDestroy( papszHeader );
92 1 : CSLDestroy( papszExtraFiles );
93 :
94 2 : for( int i = 0; i < GetRasterCount(); i++ )
95 : {
96 1 : VSIFCloseL( ((RawRasterBand *) GetRasterBand(i+1))->GetFPL() );
97 : }
98 1 : }
99 :
100 : /************************************************************************/
101 : /* GetProjectionRef() */
102 : /************************************************************************/
103 :
104 1 : const char *NDFDataset::GetProjectionRef()
105 :
106 : {
107 1 : return pszProjection;
108 : }
109 :
110 : /************************************************************************/
111 : /* GetGeoTransform() */
112 : /************************************************************************/
113 :
114 1 : CPLErr NDFDataset::GetGeoTransform( double * padfTransform )
115 :
116 : {
117 1 : memcpy( padfTransform, adfGeoTransform, sizeof(double) * 6 );
118 1 : return CE_None;
119 : }
120 :
121 : /************************************************************************/
122 : /* Get() */
123 : /* */
124 : /* Fetch a value from the header by keyword. */
125 : /************************************************************************/
126 :
127 13 : const char *NDFDataset::Get( const char *pszKey, const char *pszDefault )
128 :
129 : {
130 13 : const char *pszResult = CSLFetchNameValue( papszHeader, pszKey );
131 :
132 13 : if( pszResult == NULL )
133 0 : return pszDefault;
134 : else
135 13 : return pszResult;
136 : }
137 :
138 : /************************************************************************/
139 : /* GetFileList() */
140 : /************************************************************************/
141 :
142 0 : char **NDFDataset::GetFileList()
143 :
144 : {
145 0 : char **papszFileList = NULL;
146 :
147 : // Main data file, etc.
148 0 : papszFileList = GDALPamDataset::GetFileList();
149 :
150 : // Header file.
151 : papszFileList = CSLInsertStrings( papszFileList, -1,
152 0 : papszExtraFiles );
153 :
154 0 : return papszFileList;
155 : }
156 :
157 : /************************************************************************/
158 : /* Open() */
159 : /************************************************************************/
160 :
161 11956 : GDALDataset *NDFDataset::Open( GDALOpenInfo * poOpenInfo )
162 :
163 : {
164 : /* -------------------------------------------------------------------- */
165 : /* The user must select the header file (ie. .H1). */
166 : /* -------------------------------------------------------------------- */
167 11956 : if( poOpenInfo->nHeaderBytes < 50 )
168 11364 : return NULL;
169 :
170 592 : if( !EQUALN((const char *)poOpenInfo->pabyHeader,"NDF_REVISION=2",14)
171 : && !EQUALN((const char *)poOpenInfo->pabyHeader,"NDF_REVISION=0",14) )
172 591 : return NULL;
173 :
174 : /* -------------------------------------------------------------------- */
175 : /* Read and process the header into a local name/value */
176 : /* stringlist. We just take off the trailing semicolon. The */
177 : /* keyword is already seperated from the value by an equal */
178 : /* sign. */
179 : /* -------------------------------------------------------------------- */
180 :
181 1 : VSILFILE* fp = VSIFOpenL(poOpenInfo->pszFilename, "rb");
182 1 : if (fp == NULL)
183 0 : return NULL;
184 :
185 : const char *pszLine;
186 1 : const int nHeaderMax = 1000;
187 1 : int nHeaderLines = 0;
188 1 : char **papszHeader = (char **) CPLMalloc(sizeof(char *) * (nHeaderMax+1));
189 :
190 54 : while( nHeaderLines < nHeaderMax
191 : && (pszLine = CPLReadLineL( fp )) != NULL
192 : && !EQUAL(pszLine,"END_OF_HDR;") )
193 : {
194 : char *pszFixed;
195 :
196 52 : if( strstr(pszLine,"=") == NULL )
197 0 : break;
198 :
199 52 : pszFixed = CPLStrdup( pszLine );
200 52 : if( pszFixed[strlen(pszFixed)-1] == ';' )
201 52 : pszFixed[strlen(pszFixed)-1] = '\0';
202 :
203 52 : papszHeader[nHeaderLines++] = pszFixed;
204 52 : papszHeader[nHeaderLines] = NULL;
205 : }
206 1 : VSIFCloseL(fp);
207 1 : fp = NULL;
208 :
209 1 : if( CSLFetchNameValue( papszHeader, "PIXELS_PER_LINE" ) == NULL
210 : || CSLFetchNameValue( papszHeader, "LINES_PER_DATA_FILE" ) == NULL
211 : || CSLFetchNameValue( papszHeader, "BITS_PER_PIXEL" ) == NULL
212 : || CSLFetchNameValue( papszHeader, "PIXEL_FORMAT" ) == NULL )
213 : {
214 : CPLError( CE_Failure, CPLE_AppDefined,
215 0 : "Dataset appears to be NDF but is missing a required field.");
216 0 : CSLDestroy( papszHeader );
217 0 : return NULL;
218 : }
219 :
220 1 : if( !EQUAL(CSLFetchNameValue( papszHeader, "PIXEL_FORMAT"),
221 : "BYTE" )
222 : || !EQUAL(CSLFetchNameValue( papszHeader, "BITS_PER_PIXEL"),"8") )
223 : {
224 : CPLError( CE_Failure, CPLE_AppDefined,
225 0 : "Currently NDF driver supports only 8bit BYTE format." );
226 0 : CSLDestroy( papszHeader );
227 0 : return NULL;
228 : }
229 :
230 : /* -------------------------------------------------------------------- */
231 : /* Confirm the requested access is supported. */
232 : /* -------------------------------------------------------------------- */
233 1 : if( poOpenInfo->eAccess == GA_Update )
234 : {
235 0 : CSLDestroy( papszHeader );
236 : CPLError( CE_Failure, CPLE_NotSupported,
237 : "The NDF driver does not support update access to existing"
238 0 : " datasets.\n" );
239 0 : return NULL;
240 : }
241 :
242 : /* -------------------------------------------------------------------- */
243 : /* Create a corresponding GDALDataset. */
244 : /* -------------------------------------------------------------------- */
245 : NDFDataset *poDS;
246 :
247 1 : poDS = new NDFDataset();
248 1 : poDS->papszHeader = papszHeader;
249 :
250 1 : poDS->nRasterXSize = atoi(poDS->Get("PIXELS_PER_LINE",""));
251 1 : poDS->nRasterYSize = atoi(poDS->Get("LINES_PER_DATA_FILE",""));
252 :
253 : /* -------------------------------------------------------------------- */
254 : /* Create a raw raster band for each file. */
255 : /* -------------------------------------------------------------------- */
256 : int iBand;
257 : const char* pszBand = CSLFetchNameValue(papszHeader,
258 1 : "NUMBER_OF_BANDS_IN_VOLUME");
259 1 : if (pszBand == NULL)
260 : {
261 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot find band count");
262 0 : delete poDS;
263 0 : return NULL;
264 : }
265 1 : int nBands = atoi(pszBand);
266 :
267 1 : if (!GDALCheckDatasetDimensions(poDS->nRasterXSize, poDS->nRasterYSize) ||
268 : !GDALCheckBandCount(nBands, FALSE))
269 : {
270 0 : delete poDS;
271 0 : return NULL;
272 : }
273 :
274 1 : for( iBand = 0; iBand < nBands; iBand++ )
275 : {
276 : char szKey[100];
277 1 : CPLString osFilename;
278 :
279 1 : sprintf( szKey, "BAND%d_FILENAME", iBand+1 );
280 1 : osFilename = poDS->Get(szKey,"");
281 :
282 : // NDF1 file do not include the band filenames.
283 1 : if( osFilename.size() == 0 )
284 : {
285 : char szBandExtension[15];
286 0 : sprintf( szBandExtension, "I%d", iBand+1 );
287 : osFilename = CPLResetExtension( poOpenInfo->pszFilename,
288 0 : szBandExtension );
289 : }
290 : else
291 : {
292 1 : CPLString osBasePath = CPLGetPath(poOpenInfo->pszFilename);
293 1 : osFilename = CPLFormFilename( osBasePath, osFilename, NULL);
294 : }
295 :
296 1 : VSILFILE *fpRaw = VSIFOpenL( osFilename, "rb" );
297 1 : if( fpRaw == NULL )
298 : {
299 : CPLError( CE_Failure, CPLE_AppDefined,
300 : "Failed to open band file: %s",
301 0 : osFilename.c_str() );
302 0 : delete poDS;
303 0 : return NULL;
304 : }
305 : poDS->papszExtraFiles =
306 : CSLAddString( poDS->papszExtraFiles,
307 1 : osFilename );
308 :
309 : RawRasterBand *poBand =
310 : new RawRasterBand( poDS, iBand+1, fpRaw, 0, 1, poDS->nRasterXSize,
311 1 : GDT_Byte, TRUE, TRUE );
312 :
313 1 : sprintf( szKey, "BAND%d_NAME", iBand+1 );
314 2 : poBand->SetDescription( poDS->Get(szKey, "") );
315 :
316 1 : sprintf( szKey, "BAND%d_WAVELENGTHS", iBand+1 );
317 1 : poBand->SetMetadataItem( "WAVELENGTHS", poDS->Get(szKey,"") );
318 :
319 1 : sprintf( szKey, "BAND%d_RADIOMETRIC_GAINS/BIAS", iBand+1 );
320 : poBand->SetMetadataItem( "RADIOMETRIC_GAINS_BIAS",
321 1 : poDS->Get(szKey,"") );
322 :
323 1 : poDS->SetBand( iBand+1, poBand );
324 : }
325 :
326 : /* -------------------------------------------------------------------- */
327 : /* Fetch and parse USGS projection parameters. */
328 : /* -------------------------------------------------------------------- */
329 1 : double adfUSGSParms[15] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
330 : char **papszParmTokens =
331 : CSLTokenizeStringComplex( poDS->Get( "USGS_PROJECTION_NUMBER", "" ),
332 1 : ",", FALSE, TRUE );
333 :
334 1 : if( CSLCount( papszParmTokens ) >= 15 )
335 : {
336 : int i;
337 0 : for( i = 0; i < 15; i++ )
338 0 : adfUSGSParms[i] = atof(papszParmTokens[i]);
339 : }
340 1 : CSLDestroy(papszParmTokens);
341 1 : papszParmTokens = NULL;
342 :
343 : /* -------------------------------------------------------------------- */
344 : /* Minimal georef support ... should add full USGS style */
345 : /* support at some point. */
346 : /* -------------------------------------------------------------------- */
347 1 : OGRSpatialReference oSRS;
348 1 : int nUSGSProjection = atoi(poDS->Get( "USGS_PROJECTION_NUMBER", "" ));
349 1 : int nZone = atoi(poDS->Get("USGS_MAP_ZONE","0"));
350 :
351 1 : oSRS.importFromUSGS( nUSGSProjection, nZone, adfUSGSParms, 12 );
352 :
353 1 : CPLString osDatum = poDS->Get( "HORIZONTAL_DATUM", "" );
354 1 : if( EQUAL(osDatum,"WGS84") || EQUAL(osDatum,"NAD83")
355 : || EQUAL(osDatum,"NAD27") )
356 : {
357 1 : oSRS.SetWellKnownGeogCS( osDatum );
358 : }
359 0 : else if( EQUALN(osDatum,"NAD27",5) )
360 : {
361 0 : oSRS.SetWellKnownGeogCS( "NAD27" );
362 : }
363 : else
364 : {
365 : CPLError( CE_Warning, CPLE_AppDefined,
366 : "Unrecognised datum name in NLAPS/NDF file:%s, assuming WGS84.",
367 0 : osDatum.c_str() );
368 0 : oSRS.SetWellKnownGeogCS( "WGS84" );
369 : }
370 :
371 1 : if( oSRS.GetRoot() != NULL )
372 : {
373 1 : CPLFree( poDS->pszProjection );
374 1 : poDS->pszProjection = NULL;
375 1 : oSRS.exportToWkt( &(poDS->pszProjection) );
376 : }
377 :
378 : /* -------------------------------------------------------------------- */
379 : /* Get geotransform. */
380 : /* -------------------------------------------------------------------- */
381 : char **papszUL = CSLTokenizeString2(
382 1 : poDS->Get("UPPER_LEFT_CORNER",""), ",", 0 );
383 : char **papszUR = CSLTokenizeString2(
384 1 : poDS->Get("UPPER_RIGHT_CORNER",""), ",", 0 );
385 : char **papszLL = CSLTokenizeString2(
386 1 : poDS->Get("LOWER_LEFT_CORNER",""), ",", 0 );
387 :
388 1 : if( CSLCount(papszUL) == 4
389 : && CSLCount(papszUR) == 4
390 : && CSLCount(papszLL) == 4 )
391 : {
392 1 : poDS->adfGeoTransform[0] = atof(papszUL[2]);
393 : poDS->adfGeoTransform[1] =
394 1 : (atof(papszUR[2]) - atof(papszUL[2])) / (poDS->nRasterXSize-1);
395 : poDS->adfGeoTransform[2] =
396 1 : (atof(papszUR[3]) - atof(papszUL[3])) / (poDS->nRasterXSize-1);
397 :
398 1 : poDS->adfGeoTransform[3] = atof(papszUL[3]);
399 : poDS->adfGeoTransform[4] =
400 1 : (atof(papszLL[2]) - atof(papszUL[2])) / (poDS->nRasterYSize-1);
401 : poDS->adfGeoTransform[5] =
402 1 : (atof(papszLL[3]) - atof(papszUL[3])) / (poDS->nRasterYSize-1);
403 :
404 : // Move origin up-left half a pixel.
405 1 : poDS->adfGeoTransform[0] -= poDS->adfGeoTransform[1] * 0.5;
406 1 : poDS->adfGeoTransform[0] -= poDS->adfGeoTransform[4] * 0.5;
407 1 : poDS->adfGeoTransform[3] -= poDS->adfGeoTransform[2] * 0.5;
408 1 : poDS->adfGeoTransform[3] -= poDS->adfGeoTransform[5] * 0.5;
409 : }
410 :
411 1 : CSLDestroy( papszUL );
412 1 : CSLDestroy( papszLL );
413 1 : CSLDestroy( papszUR );
414 :
415 : /* -------------------------------------------------------------------- */
416 : /* Initialize any PAM information. */
417 : /* -------------------------------------------------------------------- */
418 1 : poDS->SetDescription( poOpenInfo->pszFilename );
419 1 : poDS->TryLoadXML();
420 :
421 : /* -------------------------------------------------------------------- */
422 : /* Check for overviews. */
423 : /* -------------------------------------------------------------------- */
424 1 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
425 :
426 1 : return( poDS );
427 : }
428 :
429 : /************************************************************************/
430 : /* GDALRegister_NDF() */
431 : /************************************************************************/
432 :
433 582 : void GDALRegister_NDF()
434 :
435 : {
436 : GDALDriver *poDriver;
437 :
438 582 : if( GDALGetDriverByName( "NDF" ) == NULL )
439 : {
440 561 : poDriver = new GDALDriver();
441 :
442 561 : poDriver->SetDescription( "NDF" );
443 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
444 561 : "NLAPS Data Format" );
445 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
446 561 : "frmt_various.html#NDF" );
447 561 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
448 :
449 561 : poDriver->pfnOpen = NDFDataset::Open;
450 :
451 561 : GetGDALDriverManager()->RegisterDriver( poDriver );
452 : }
453 582 : }
454 :
|
