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