1 : /******************************************************************************
2 : * $Id: grddataset.cpp 18354 2009-12-20 10:42:15Z rouault $
3 : *
4 : * Project: GRD Reader
5 : * Purpose: GDAL driver for Northwood Grid Format
6 : * Author: Perry Casson
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2006, Waypoint Information Technology
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 "gdal_pam.h"
31 : #include "northwood.h"
32 :
33 : #ifdef OGR_ENABLED
34 : #ifdef MSVC
35 : #include "..\..\ogr\ogrsf_frmts\mitab\mitab.h"
36 : #else
37 : #include "../../ogr/ogrsf_frmts/mitab/mitab.h"
38 : #endif
39 : #endif
40 :
41 :
42 : CPL_C_START void GDALRegister_NWT_GRD( void );
43 : CPL_C_END
44 : /************************************************************************/
45 : /* ==================================================================== */
46 : /* NWT_GRDDataset */
47 : /* ==================================================================== */
48 : /************************************************************************/
49 : class NWT_GRDRasterBand;
50 :
51 : class NWT_GRDDataset:public GDALPamDataset
52 : {
53 : friend class NWT_GRDRasterBand;
54 :
55 : FILE *fp;
56 : GByte abyHeader[1024];
57 : NWT_GRID *pGrd;
58 : NWT_RGB ColorMap[4096];
59 : char *pszProjection;
60 :
61 : public:
62 : NWT_GRDDataset();
63 : ~NWT_GRDDataset();
64 :
65 : static GDALDataset *Open( GDALOpenInfo * );
66 :
67 : CPLErr GetGeoTransform( double *padfTransform );
68 : const char *GetProjectionRef();
69 : };
70 :
71 : /************************************************************************/
72 : /* ==================================================================== */
73 : /* NWT_GRDRasterBand */
74 : /* ==================================================================== */
75 : /************************************************************************/
76 :
77 : class NWT_GRDRasterBand:public GDALPamRasterBand
78 24 : {
79 : friend class NWT_GRDDataset;
80 :
81 : int bHaveOffsetScale;
82 : double dfOffset;
83 : double dfScale;
84 :
85 : public:
86 :
87 : NWT_GRDRasterBand( NWT_GRDDataset *, int );
88 :
89 : virtual CPLErr IReadBlock( int, int, void * );
90 : virtual double GetNoDataValue( int *pbSuccess );
91 :
92 : virtual double GetOffset( int *pbSuccess = NULL );
93 : virtual CPLErr SetOffset( double dfNewValue );
94 : virtual double GetScale( int *pbSuccess = NULL );
95 : virtual CPLErr SetScale( double dfNewValue );
96 :
97 : virtual GDALColorInterp GetColorInterpretation();
98 : };
99 :
100 :
101 : /************************************************************************/
102 : /* NWT_GRDRasterBand() */
103 : /************************************************************************/
104 12 : NWT_GRDRasterBand::NWT_GRDRasterBand( NWT_GRDDataset * poDS, int nBand )
105 : {
106 12 : this->poDS = poDS;
107 12 : this->nBand = nBand;
108 :
109 12 : if( nBand == 4 )
110 : {
111 3 : bHaveOffsetScale = TRUE;
112 3 : dfOffset = poDS->pGrd->fZMin;
113 :
114 3 : if( poDS->pGrd->cFormat == 0x01 )
115 : {
116 0 : eDataType = GDT_Float32;
117 0 : dfScale =( poDS->pGrd->fZMax - poDS->pGrd->fZMin ) / 4294967294.0;
118 : }
119 : else
120 : {
121 3 : eDataType = GDT_Float32;
122 3 : dfScale =( poDS->pGrd->fZMax - poDS->pGrd->fZMin ) / 65534.0;
123 : }
124 : }
125 : else
126 : {
127 9 : bHaveOffsetScale = FALSE;
128 9 : dfOffset = 0;
129 9 : dfScale = 1.0;
130 9 : eDataType = GDT_Byte;
131 : }
132 12 : nBlockXSize = poDS->GetRasterXSize();
133 12 : nBlockYSize = 1;
134 12 : }
135 :
136 0 : double NWT_GRDRasterBand::GetNoDataValue( int *pbSuccess )
137 : {
138 0 : if( pbSuccess != NULL )
139 0 : *pbSuccess = TRUE;
140 :
141 0 : return 0; //Northwood grid 0 is always null
142 : }
143 :
144 0 : GDALColorInterp NWT_GRDRasterBand::GetColorInterpretation()
145 : {
146 : //return GCI_RGB;
147 0 : if( nBand == 4 )
148 0 : return GCI_Undefined;
149 0 : else if( nBand == 1 )
150 0 : return GCI_RedBand;
151 0 : else if( nBand == 2 )
152 0 : return GCI_GreenBand;
153 0 : else if( nBand == 3 )
154 0 : return GCI_BlueBand;
155 : else
156 0 : return GCI_Undefined;
157 : }
158 :
159 : /************************************************************************/
160 : /* IReadBlock() */
161 : /************************************************************************/
162 141 : CPLErr NWT_GRDRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff, void *pImage )
163 : {
164 141 : NWT_GRDDataset *poGDS = (NWT_GRDDataset *) poDS;
165 : char *pszRecord;
166 141 : int nRecordSize = nBlockXSize * 2;
167 : int i;
168 : unsigned short raw1;
169 :
170 141 : VSIFSeek( poGDS->fp, 1024 + nRecordSize * nBlockYOff, SEEK_SET );
171 :
172 141 : pszRecord = (char *) CPLMalloc( nRecordSize );
173 141 : VSIFRead( pszRecord, 1, nRecordSize, poGDS->fp );
174 :
175 141 : if( nBand == 4 ) //Z values
176 : {
177 0 : for( i = 0; i < nBlockXSize; i++ )
178 : {
179 0 : memcpy( (void *) &raw1, (void *)(pszRecord + 2 * i), 2 );
180 : CPL_LSBPTR16(&raw1);
181 0 : if( raw1 == 0 )
182 : {
183 0 : ((float *)pImage)[i] = -1.e37; // null value
184 : }
185 : else
186 : {
187 0 : ((float *)pImage)[i] = dfOffset + ((raw1 - 1) * dfScale);
188 : }
189 : }
190 : }
191 141 : else if( nBand == 1 ) // red values
192 : {
193 3008 : for( i = 0; i < nBlockXSize; i++ )
194 : {
195 2961 : memcpy( (void *) &raw1, (void *)(pszRecord + 2 * i), 2 );
196 : CPL_LSBPTR16(&raw1);
197 2961 : ((char *)pImage)[i] = poGDS->ColorMap[raw1 / 16].r;
198 : }
199 : }
200 94 : else if( nBand == 2 ) // green
201 : {
202 3008 : for( i = 0; i < nBlockXSize; i++ )
203 : {
204 2961 : memcpy( (void *) &raw1, (void *)(pszRecord + 2 * i), 2 );
205 : CPL_LSBPTR16(&raw1);
206 2961 : ((char *) pImage)[i] = poGDS->ColorMap[raw1 / 16].g;
207 : }
208 : }
209 47 : else if( nBand == 3 ) // blue
210 : {
211 3008 : for( i = 0; i < nBlockXSize; i++ )
212 : {
213 2961 : memcpy( (void *) &raw1, (void *) (pszRecord + 2 * i), 2 );
214 : CPL_LSBPTR16(&raw1);
215 2961 : ((char *) pImage)[i] = poGDS->ColorMap[raw1 / 16].b;
216 : }
217 : }
218 : else
219 : {
220 : CPLError( CE_Failure, CPLE_IllegalArg,
221 : "No band number %d",
222 0 : nBand );
223 0 : if( pszRecord != NULL )
224 0 : CPLFree( pszRecord );
225 0 : return CE_Failure;
226 : }
227 141 : if( pszRecord != NULL )
228 : {
229 141 : CPLFree( pszRecord );
230 : }
231 141 : return CE_None;
232 : }
233 :
234 : /************************************************************************/
235 : /* GetOffset() */
236 : /************************************************************************/
237 0 : double NWT_GRDRasterBand::GetOffset( int *pbSuccess )
238 : {
239 0 : if( pbSuccess )
240 0 : *pbSuccess = bHaveOffsetScale;
241 0 : return dfOffset;
242 : }
243 :
244 : /************************************************************************/
245 : /* SetOffset() */
246 : /************************************************************************/
247 0 : CPLErr NWT_GRDRasterBand::SetOffset( double dfNewValue )
248 : {
249 : //poGDS->bMetadataChanged = TRUE;
250 :
251 0 : bHaveOffsetScale = TRUE;
252 0 : dfOffset = dfNewValue;
253 0 : return CE_None;
254 : }
255 :
256 : /************************************************************************/
257 : /* GetScale() */
258 : /************************************************************************/
259 0 : double NWT_GRDRasterBand::GetScale( int *pbSuccess )
260 : {
261 0 : if( pbSuccess )
262 0 : *pbSuccess = bHaveOffsetScale;
263 0 : return dfScale;
264 : }
265 :
266 : /************************************************************************/
267 : /* SetScale() */
268 : /************************************************************************/
269 0 : CPLErr NWT_GRDRasterBand::SetScale( double dfNewValue )
270 : {
271 0 : bHaveOffsetScale = TRUE;
272 0 : dfScale = dfNewValue;
273 0 : return CE_None;
274 : }
275 :
276 : /************************************************************************/
277 : /* ==================================================================== */
278 : /* NWT_GRDDataset */
279 : /* ==================================================================== */
280 : /************************************************************************/
281 3 : NWT_GRDDataset::NWT_GRDDataset()
282 : {
283 3 : pszProjection = NULL;
284 : //poCT = NULL;
285 3 : }
286 :
287 :
288 : /************************************************************************/
289 : /* ~NWT_GRDDataset() */
290 : /************************************************************************/
291 :
292 6 : NWT_GRDDataset::~NWT_GRDDataset()
293 : {
294 3 : FlushCache();
295 3 : pGrd->fp = NULL; // this prevents nwtCloseGrid from closing the fp
296 3 : nwtCloseGrid( pGrd );
297 :
298 3 : if( fp != NULL )
299 3 : VSIFClose( fp );
300 :
301 3 : if( pszProjection != NULL )
302 : {
303 0 : CPLFree( pszProjection );
304 : }
305 : /*if( poCT != NULL )
306 : delete poCT;*/
307 6 : }
308 :
309 : /************************************************************************/
310 : /* GetGeoTransform() */
311 : /************************************************************************/
312 0 : CPLErr NWT_GRDDataset::GetGeoTransform( double *padfTransform )
313 : {
314 0 : padfTransform[0] = pGrd->dfMinX - ( pGrd->dfStepSize * 0.5 );
315 0 : padfTransform[3] = pGrd->dfMaxY + ( pGrd->dfStepSize * 0.5 );
316 0 : padfTransform[1] = pGrd->dfStepSize;
317 0 : padfTransform[2] = 0.0;
318 :
319 0 : padfTransform[4] = 0.0;
320 0 : padfTransform[5] = -1 * pGrd->dfStepSize;
321 :
322 0 : return CE_None;
323 : }
324 :
325 : /************************************************************************/
326 : /* GetProjectionRef() */
327 : /************************************************************************/
328 0 : const char *NWT_GRDDataset::GetProjectionRef()
329 : {
330 : #ifdef OGR_ENABLED
331 0 : if (pszProjection == NULL)
332 : {
333 : OGRSpatialReference *poSpatialRef;
334 0 : poSpatialRef = MITABCoordSys2SpatialRef( pGrd->cMICoordSys );
335 0 : if (poSpatialRef)
336 : {
337 0 : poSpatialRef->exportToWkt( &pszProjection );
338 0 : poSpatialRef->Release();
339 : }
340 : }
341 : #endif
342 0 : return( (const char *)pszProjection );
343 : }
344 :
345 : /************************************************************************/
346 : /* Open() */
347 : /************************************************************************/
348 8447 : GDALDataset *NWT_GRDDataset::Open( GDALOpenInfo * poOpenInfo )
349 : {
350 : /* -------------------------------------------------------------------- */
351 : /* Look for the header */
352 : /* -------------------------------------------------------------------- */
353 8447 : if( poOpenInfo->fp == NULL || poOpenInfo->nHeaderBytes < 50 )
354 8315 : return NULL;
355 :
356 148 : if( poOpenInfo->pabyHeader[0] != 'H' ||
357 4 : poOpenInfo->pabyHeader[1] != 'G' ||
358 4 : poOpenInfo->pabyHeader[2] != 'P' ||
359 8 : poOpenInfo->pabyHeader[3] != 'C' || poOpenInfo->pabyHeader[4] != '1' )
360 129 : return NULL;
361 :
362 : /* -------------------------------------------------------------------- */
363 : /* Create a corresponding GDALDataset. */
364 : /* -------------------------------------------------------------------- */
365 : NWT_GRDDataset *poDS;
366 :
367 3 : poDS = new NWT_GRDDataset();
368 :
369 3 : poDS->fp = poOpenInfo->fp;
370 3 : poOpenInfo->fp = NULL;
371 :
372 : /* -------------------------------------------------------------------- */
373 : /* Read the header. */
374 : /* -------------------------------------------------------------------- */
375 3 : VSIFSeek( poDS->fp, 0, SEEK_SET );
376 3 : VSIFRead( poDS->abyHeader, 1, 1024, poDS->fp );
377 3 : poDS->pGrd = (NWT_GRID *) malloc(sizeof(NWT_GRID));
378 6 : if (!nwt_ParseHeader( poDS->pGrd, (char *) poDS->abyHeader ) ||
379 : !GDALCheckDatasetDimensions(poDS->pGrd->nXSide, poDS->pGrd->nYSide) )
380 : {
381 0 : delete poDS;
382 0 : return NULL;
383 : }
384 :
385 3 : poDS->nRasterXSize = poDS->pGrd->nXSide;
386 3 : poDS->nRasterYSize = poDS->pGrd->nYSide;
387 :
388 : // create a colorTable
389 : // if( poDS->pGrd->iNumColorInflections > 0 )
390 : // poDS->CreateColorTable();
391 3 : nwt_LoadColors( poDS->ColorMap, 4096, poDS->pGrd );
392 : /* -------------------------------------------------------------------- */
393 : /* Create band information objects. */
394 : /* -------------------------------------------------------------------- */
395 3 : poDS->SetBand( 1, new NWT_GRDRasterBand( poDS, 1 ) ); //r
396 6 : poDS->SetBand( 2, new NWT_GRDRasterBand( poDS, 2 ) ); //g
397 6 : poDS->SetBand( 3, new NWT_GRDRasterBand( poDS, 3 ) ); //b
398 6 : poDS->SetBand( 4, new NWT_GRDRasterBand( poDS, 4 ) ); //z
399 :
400 : /* -------------------------------------------------------------------- */
401 : /* Initialize any PAM information. */
402 : /* -------------------------------------------------------------------- */
403 3 : poDS->SetDescription( poOpenInfo->pszFilename );
404 3 : poDS->TryLoadXML();
405 :
406 3 : return (poDS);
407 : }
408 :
409 :
410 : /************************************************************************/
411 : /* GDALRegister_GRD() */
412 : /************************************************************************/
413 338 : void GDALRegister_NWT_GRD()
414 : {
415 : GDALDriver *poDriver;
416 :
417 338 : if( GDALGetDriverByName( "NWT_GRD" ) == NULL )
418 : {
419 336 : poDriver = new GDALDriver();
420 :
421 336 : poDriver->SetDescription( "NWT_GRD" );
422 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
423 336 : "Northwood Numeric Grid Format .grd/.tab" );
424 336 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC, "frmt_various.html#grd");
425 336 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "grd" );
426 :
427 336 : poDriver->pfnOpen = NWT_GRDDataset::Open;
428 :
429 336 : GetGDALDriverManager()->RegisterDriver( poDriver );
430 : }
431 338 : }
|
