1 : /******************************************************************************
2 : * $Id: sdtsrasterreader.cpp 17213 2009-06-07 11:55:50Z rouault $
3 : *
4 : * Project: SDTS Translator
5 : * Purpose: Implementation of SDTSRasterReader class.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 1999, 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 "sdts_al.h"
31 :
32 : CPL_CVSID("$Id: sdtsrasterreader.cpp 17213 2009-06-07 11:55:50Z rouault $");
33 :
34 : /************************************************************************/
35 : /* SDTSRasterReader() */
36 : /************************************************************************/
37 :
38 2 : SDTSRasterReader::SDTSRasterReader()
39 :
40 : {
41 2 : nXSize = 0;
42 2 : nYSize = 0;
43 2 : nXBlockSize = 0;
44 2 : nYBlockSize = 0;
45 2 : nXStart = 0;
46 2 : nYStart = 0;
47 :
48 2 : strcpy( szINTR, "CE" );
49 2 : }
50 :
51 : /************************************************************************/
52 : /* ~SDTSRasterReader() */
53 : /************************************************************************/
54 :
55 2 : SDTSRasterReader::~SDTSRasterReader()
56 : {
57 2 : }
58 :
59 : /************************************************************************/
60 : /* Close() */
61 : /************************************************************************/
62 :
63 0 : void SDTSRasterReader::Close()
64 :
65 : {
66 0 : oDDFModule.Close();
67 0 : }
68 :
69 : /************************************************************************/
70 : /* Open() */
71 : /* */
72 : /* Open the requested cell file, and collect required */
73 : /* information. */
74 : /************************************************************************/
75 :
76 2 : int SDTSRasterReader::Open( SDTS_CATD * poCATD, SDTS_IREF * poIREF,
77 : const char * pszModule )
78 :
79 : {
80 2 : strncpy( szModule, pszModule, sizeof(szModule) );
81 2 : szModule[sizeof(szModule) - 1] = '\0';
82 :
83 : /* ==================================================================== */
84 : /* Search the LDEF module for the requested cell module. */
85 : /* ==================================================================== */
86 2 : DDFModule oLDEF;
87 : DDFRecord *poRecord;
88 :
89 : /* -------------------------------------------------------------------- */
90 : /* Open the LDEF module, and report failure if it is missing. */
91 : /* -------------------------------------------------------------------- */
92 2 : if( poCATD->GetModuleFilePath("LDEF") == NULL )
93 : {
94 : CPLError( CE_Failure, CPLE_AppDefined,
95 : "Can't find LDEF entry in CATD module ... "
96 0 : "can't treat as raster.\n" );
97 0 : return FALSE;
98 : }
99 :
100 2 : if( !oLDEF.Open( poCATD->GetModuleFilePath("LDEF") ) )
101 0 : return FALSE;
102 :
103 : /* -------------------------------------------------------------------- */
104 : /* Read each record, till we find what we want. */
105 : /* -------------------------------------------------------------------- */
106 4 : while( (poRecord = oLDEF.ReadRecord() ) != NULL )
107 : {
108 2 : if( EQUAL(poRecord->GetStringSubfield("LDEF",0,"CMNM",0), pszModule) )
109 2 : break;
110 : }
111 :
112 2 : if( poRecord == NULL )
113 : {
114 : CPLError( CE_Failure, CPLE_AppDefined,
115 : "Can't find module `%s' in LDEF file.\n",
116 0 : pszModule );
117 0 : return FALSE;
118 : }
119 :
120 : /* -------------------------------------------------------------------- */
121 : /* Extract raster dimensions, and origin offset (0/1). */
122 : /* -------------------------------------------------------------------- */
123 2 : nXSize = poRecord->GetIntSubfield( "LDEF", 0, "NCOL", 0 );
124 2 : nYSize = poRecord->GetIntSubfield( "LDEF", 0, "NROW", 0 );
125 :
126 2 : nXStart = poRecord->GetIntSubfield( "LDEF", 0, "SOCI", 0 );
127 2 : nYStart = poRecord->GetIntSubfield( "LDEF", 0, "SORI", 0 );
128 :
129 : /* -------------------------------------------------------------------- */
130 : /* Get the point in the pixel that the origin defines. We only */
131 : /* support top left and center. */
132 : /* -------------------------------------------------------------------- */
133 2 : strcpy( szINTR, poRecord->GetStringSubfield( "LDEF", 0, "INTR", 0 ) );
134 2 : if( EQUAL(szINTR,"") )
135 0 : strcpy( szINTR, "CE" );
136 :
137 2 : if( !EQUAL(szINTR,"CE") && !EQUAL(szINTR,"TL") )
138 : {
139 : CPLError( CE_Warning, CPLE_AppDefined,
140 : "Unsupported INTR value of `%s', assume CE.\n"
141 : "Positions may be off by one pixel.\n",
142 0 : szINTR );
143 0 : strcpy( szINTR, "CE" );
144 : }
145 :
146 : /* -------------------------------------------------------------------- */
147 : /* Record the LDEF record number we used so we can find the */
148 : /* corresponding RSDF record. */
149 : /* -------------------------------------------------------------------- */
150 : int nLDEF_RCID;
151 :
152 2 : nLDEF_RCID = poRecord->GetIntSubfield( "LDEF", 0, "RCID", 0 );
153 :
154 2 : oLDEF.Close();
155 :
156 : /* ==================================================================== */
157 : /* Search the RSDF module for the requested cell module. */
158 : /* ==================================================================== */
159 2 : DDFModule oRSDF;
160 :
161 : /* -------------------------------------------------------------------- */
162 : /* Open the RSDF module, and report failure if it is missing. */
163 : /* -------------------------------------------------------------------- */
164 2 : if( poCATD->GetModuleFilePath("RSDF") == NULL )
165 : {
166 : CPLError( CE_Failure, CPLE_AppDefined,
167 : "Can't find RSDF entry in CATD module ... "
168 0 : "can't treat as raster.\n" );
169 0 : return FALSE;
170 : }
171 :
172 2 : if( !oRSDF.Open( poCATD->GetModuleFilePath("RSDF") ) )
173 0 : return FALSE;
174 :
175 : /* -------------------------------------------------------------------- */
176 : /* Read each record, till we find what we want. */
177 : /* -------------------------------------------------------------------- */
178 4 : while( (poRecord = oRSDF.ReadRecord() ) != NULL )
179 : {
180 2 : if( poRecord->GetIntSubfield("LYID",0,"RCID",0) == nLDEF_RCID )
181 2 : break;
182 : }
183 :
184 2 : if( poRecord == NULL )
185 : {
186 : CPLError( CE_Failure, CPLE_AppDefined,
187 : "Can't find LDEF:%d record in RSDF file.\n",
188 0 : nLDEF_RCID );
189 0 : return FALSE;
190 : }
191 :
192 : /* -------------------------------------------------------------------- */
193 : /* Establish the raster pixel/line to georef transformation. */
194 : /* -------------------------------------------------------------------- */
195 : double dfZ;
196 :
197 2 : if( poRecord->FindField( "SADR" ) == NULL )
198 : {
199 : CPLError( CE_Failure, CPLE_AppDefined,
200 0 : "Can't find SADR field in RSDF record.\n" );
201 0 : return FALSE;
202 : }
203 :
204 : poIREF->GetSADR( poRecord->FindField( "SADR" ), 1,
205 2 : adfTransform + 0, adfTransform + 3, &dfZ );
206 :
207 2 : adfTransform[1] = poIREF->dfXRes;
208 2 : adfTransform[2] = 0.0;
209 2 : adfTransform[4] = 0.0;
210 2 : adfTransform[5] = -1 * poIREF->dfYRes;
211 :
212 : /* -------------------------------------------------------------------- */
213 : /* If the origin is the center of the pixel, then shift it back */
214 : /* half a pixel to the top left of the top left. */
215 : /* -------------------------------------------------------------------- */
216 2 : if( EQUAL(szINTR,"CE") )
217 : {
218 2 : adfTransform[0] -= adfTransform[1] * 0.5;
219 2 : adfTransform[3] -= adfTransform[5] * 0.5;
220 : }
221 :
222 : /* -------------------------------------------------------------------- */
223 : /* Verify some other assumptions. */
224 : /* -------------------------------------------------------------------- */
225 : const char *pszString;
226 :
227 2 : pszString = poRecord->GetStringSubfield( "RSDF", 0, "OBRP", 0);
228 2 : if( !EQUAL(pszString,"G2") )
229 : {
230 : CPLError( CE_Failure, CPLE_AppDefined,
231 : "OBRP value of `%s' not expected 2D raster code (G2).\n",
232 0 : pszString );
233 0 : return FALSE;
234 : }
235 :
236 2 : pszString = poRecord->GetStringSubfield( "RSDF", 0, "SCOR", 0);
237 2 : if( !EQUAL(pszString,"TL") )
238 : {
239 : CPLError( CE_Warning, CPLE_AppDefined,
240 : "SCOR (origin) is `%s' instead of expected top left.\n"
241 : "Georef coordinates will likely be incorrect.\n",
242 0 : pszString );
243 : }
244 :
245 2 : oRSDF.Close();
246 :
247 : /* -------------------------------------------------------------------- */
248 : /* For now we will assume that the block size is one scanline. */
249 : /* We will blow a gasket later while reading the cell file if */
250 : /* this isn't the case. */
251 : /* */
252 : /* This isn't a very flexible raster implementation! */
253 : /* -------------------------------------------------------------------- */
254 2 : nXBlockSize = nXSize;
255 2 : nYBlockSize = 1;
256 :
257 : /* ==================================================================== */
258 : /* Fetch the data type used for the raster, and the units from */
259 : /* the data dictionary/schema record (DDSH). */
260 : /* ==================================================================== */
261 2 : DDFModule oDDSH;
262 :
263 : /* -------------------------------------------------------------------- */
264 : /* Open the DDSH module, and report failure if it is missing. */
265 : /* -------------------------------------------------------------------- */
266 2 : if( poCATD->GetModuleFilePath("DDSH") == NULL )
267 : {
268 : CPLError( CE_Failure, CPLE_AppDefined,
269 : "Can't find DDSH entry in CATD module ... "
270 0 : "can't treat as raster.\n" );
271 0 : return FALSE;
272 : }
273 :
274 2 : if( !oDDSH.Open( poCATD->GetModuleFilePath("DDSH") ) )
275 0 : return FALSE;
276 :
277 : /* -------------------------------------------------------------------- */
278 : /* Read each record, till we find what we want. */
279 : /* -------------------------------------------------------------------- */
280 4 : while( (poRecord = oDDSH.ReadRecord() ) != NULL )
281 : {
282 2 : if( EQUAL(poRecord->GetStringSubfield("DDSH",0,"NAME",0),pszModule) )
283 2 : break;
284 : }
285 :
286 2 : if( poRecord == NULL )
287 : {
288 : CPLError( CE_Failure, CPLE_AppDefined,
289 : "Can't find DDSH record for %s.\n",
290 0 : pszModule );
291 0 : return FALSE;
292 : }
293 :
294 : /* -------------------------------------------------------------------- */
295 : /* Get some values we are interested in. */
296 : /* -------------------------------------------------------------------- */
297 2 : if( poRecord->GetStringSubfield("DDSH",0,"FMT",0) != NULL )
298 2 : strcpy( szFMT, poRecord->GetStringSubfield("DDSH",0,"FMT",0) );
299 : else
300 0 : strcpy( szFMT, "BUI16" );
301 :
302 2 : if( poRecord->GetStringSubfield("DDSH",0,"UNIT",0) != NULL )
303 2 : strcpy( szUNITS, poRecord->GetStringSubfield("DDSH",0,"UNIT",0) );
304 : else
305 0 : strcpy( szUNITS, "METERS" );
306 :
307 2 : if( poRecord->GetStringSubfield("DDSH",0,"ATLB",0) != NULL )
308 2 : strcpy( szLabel, poRecord->GetStringSubfield("DDSH",0,"ATLB",0) );
309 : else
310 0 : strcpy( szLabel, "" );
311 :
312 : /* -------------------------------------------------------------------- */
313 : /* Open the cell file. */
314 : /* -------------------------------------------------------------------- */
315 2 : return( oDDFModule.Open( poCATD->GetModuleFilePath(pszModule) ) );
316 : }
317 :
318 : /************************************************************************/
319 : /* GetBlock() */
320 : /* */
321 : /* Read a requested block of raster data from the file. */
322 : /* */
323 : /* Currently we will always use sequential access. In the */
324 : /* future we should modify the iso8211 library to support */
325 : /* seeking, and modify this to seek directly to the right */
326 : /* record once it's location is known. */
327 : /************************************************************************/
328 :
329 : /**
330 : Read a block of raster data from the file.
331 :
332 : @param nXOffset X block offset into the file. Normally zero for scanline
333 : organized raster files.
334 :
335 : @param nYOffset Y block offset into the file. Normally the scanline offset
336 : from top of raster for scanline organized raster files.
337 :
338 : @param pData pointer to GInt16 (signed short) buffer of data into which to
339 : read the raster.
340 :
341 : @return TRUE on success and FALSE on error.
342 :
343 : */
344 :
345 25 : int SDTSRasterReader::GetBlock( int nXOffset, int nYOffset, void * pData )
346 :
347 : {
348 : DDFRecord *poRecord;
349 : int nBytesPerValue;
350 : int iTry;
351 :
352 : CPLAssert( nXOffset == 0 );
353 :
354 : /* -------------------------------------------------------------------- */
355 : /* Analyse the datatype. */
356 : /* -------------------------------------------------------------------- */
357 : CPLAssert( EQUAL(szFMT,"BI16") || EQUAL(szFMT,"BFP32") );
358 :
359 25 : if( EQUAL(szFMT,"BI16") )
360 25 : nBytesPerValue = 2;
361 : else
362 0 : nBytesPerValue = 4;
363 :
364 50 : for(iTry=0;iTry<2;iTry++)
365 : {
366 : /* -------------------------------------------------------------------- */
367 : /* Read through till we find the desired record. */
368 : /* -------------------------------------------------------------------- */
369 25 : CPLErrorReset();
370 25 : while( (poRecord = oDDFModule.ReadRecord()) != NULL )
371 : {
372 25 : if( poRecord->GetIntSubfield( "CELL", 0, "ROWI", 0 )
373 : == nYOffset + nYStart )
374 : {
375 25 : break;
376 : }
377 : }
378 :
379 25 : if( CPLGetLastErrorType() == CE_Failure )
380 0 : return FALSE;
381 :
382 : /* -------------------------------------------------------------------- */
383 : /* If we didn't get what we needed just start over. */
384 : /* -------------------------------------------------------------------- */
385 25 : if( poRecord == NULL )
386 : {
387 0 : if (iTry == 0)
388 0 : oDDFModule.Rewind();
389 : else
390 : {
391 : CPLError( CE_Failure, CPLE_AppDefined,
392 : "Cannot read scanline %d. Raster access failed.\n",
393 0 : nYOffset );
394 0 : return FALSE;
395 : }
396 : }
397 : else
398 : {
399 25 : break;
400 : }
401 : }
402 :
403 : /* -------------------------------------------------------------------- */
404 : /* Validate the records size. Does it represent exactly one */
405 : /* scanline? */
406 : /* -------------------------------------------------------------------- */
407 : DDFField *poCVLS;
408 :
409 25 : poCVLS = poRecord->FindField( "CVLS" );
410 25 : if( poCVLS == NULL )
411 0 : return FALSE;
412 :
413 25 : if( poCVLS->GetRepeatCount() != nXSize )
414 : {
415 : CPLError( CE_Failure, CPLE_AppDefined,
416 : "Cell record is %d long, but we expected %d, the number\n"
417 : "of pixels in a scanline. Raster access failed.\n",
418 0 : poCVLS->GetRepeatCount(), nXSize );
419 0 : return FALSE;
420 : }
421 :
422 : /* -------------------------------------------------------------------- */
423 : /* Does the CVLS field consist of exactly 1 B(16) field? */
424 : /* -------------------------------------------------------------------- */
425 25 : if( poCVLS->GetDataSize() < nBytesPerValue * nXSize
426 : || poCVLS->GetDataSize() > nBytesPerValue * nXSize + 1 )
427 : {
428 : CPLError( CE_Failure, CPLE_AppDefined,
429 : "Cell record is not of expected format. Raster access "
430 0 : "failed.\n" );
431 :
432 0 : return FALSE;
433 : }
434 :
435 : /* -------------------------------------------------------------------- */
436 : /* Copy the data to the application buffer, and byte swap if */
437 : /* required. */
438 : /* -------------------------------------------------------------------- */
439 25 : memcpy( pData, poCVLS->GetData(), nXSize * nBytesPerValue );
440 :
441 : #ifdef CPL_LSB
442 25 : if( nBytesPerValue == 2 )
443 : {
444 8500 : for( int i = 0; i < nXSize; i++ )
445 : {
446 8475 : ((GInt16 *) pData)[i] = CPL_MSBWORD16(((GInt16 *) pData)[i]);
447 : }
448 : }
449 : else
450 : {
451 0 : for( int i = 0; i < nXSize; i++ )
452 : {
453 0 : CPL_MSBPTR32( ((GByte *)pData) + i*4 );
454 : }
455 : }
456 : #endif
457 :
458 25 : return TRUE;
459 : }
460 :
461 : /************************************************************************/
462 : /* GetTransform() */
463 : /************************************************************************/
464 :
465 : /**
466 : Fetch the transformation between pixel/line coordinates and georeferenced
467 : coordinates.
468 :
469 : @param padfTransformOut pointer to an array of six doubles which will be
470 : filled with the georeferencing transform.
471 :
472 : @return TRUE is returned, indicating success.
473 :
474 : The padfTransformOut array consists of six values. The pixel/line coordinate
475 : (Xp,Yp) can be related to a georeferenced coordinate (Xg,Yg) or (Easting,
476 : Northing).
477 :
478 : <pre>
479 : Xg = padfTransformOut[0] + Xp * padfTransform[1] + Yp * padfTransform[2]
480 : Yg = padfTransformOut[3] + Xp * padfTransform[4] + Yp * padfTransform[5]
481 : </pre>
482 :
483 : In other words, for a north up image the top left corner of the top left
484 : pixel is at georeferenced coordinate (padfTransform[0],padfTransform[3])
485 : the pixel width is padfTransform[1], the pixel height is padfTransform[5]
486 : and padfTransform[2] and padfTransform[4] will be zero.
487 :
488 : */
489 :
490 1 : int SDTSRasterReader::GetTransform( double * padfTransformOut )
491 :
492 : {
493 1 : memcpy( padfTransformOut, adfTransform, sizeof(double)*6 );
494 :
495 1 : return TRUE;
496 : }
497 :
498 : /************************************************************************/
499 : /* GetRasterType() */
500 : /************************************************************************/
501 :
502 : /**
503 : * Fetch the pixel data type.
504 : *
505 : * Returns one of SDTS_RT_INT16 (1) or SDTS_RT_FLOAT32 (6) indicating the
506 : * type of buffer that should be passed to GetBlock().
507 : */
508 :
509 2 : int SDTSRasterReader::GetRasterType()
510 :
511 : {
512 2 : if( EQUAL(szFMT,"BFP32") )
513 0 : return 6;
514 : else
515 2 : return 1;
516 : }
517 :
518 : /************************************************************************/
519 : /* GetMinMax() */
520 : /************************************************************************/
521 :
522 : /**
523 : * Fetch the minimum and maximum raster values that occur in the file.
524 : *
525 : * Note this operation current results in a scan of the entire file.
526 : *
527 : * @param pdfMin variable in which the minimum value encountered is returned.
528 : * @param pdfMax variable in which the maximum value encountered is returned.
529 : * @param dfNoData a value to ignore when computing min/max, defaults to
530 : * -32766.
531 : *
532 : * @return TRUE on success, or FALSE if an error occurs.
533 : */
534 :
535 0 : int SDTSRasterReader::GetMinMax( double * pdfMin, double * pdfMax,
536 : double dfNoData )
537 :
538 : {
539 : void *pBuffer;
540 0 : int bFirst = TRUE;
541 0 : int b32Bit = GetRasterType() == SDTS_RT_FLOAT32;
542 :
543 : CPLAssert( GetBlockXSize() == GetXSize() && GetBlockYSize() == 1 );
544 :
545 0 : pBuffer = CPLMalloc(sizeof(float) * GetXSize());
546 :
547 0 : for( int iLine = 0; iLine < GetYSize(); iLine++ )
548 : {
549 0 : if( !GetBlock( 0, iLine, pBuffer ) )
550 : {
551 0 : CPLFree( pBuffer );
552 0 : return FALSE;
553 : }
554 :
555 0 : for( int iPixel = 0; iPixel < GetXSize(); iPixel++ )
556 : {
557 : double dfValue;
558 :
559 0 : if( b32Bit )
560 0 : dfValue = ((float *) pBuffer)[iPixel];
561 : else
562 0 : dfValue = ((short *) pBuffer)[iPixel];
563 :
564 0 : if( dfValue != dfNoData )
565 : {
566 0 : if( bFirst )
567 : {
568 0 : *pdfMin = *pdfMax = dfValue;
569 0 : bFirst = FALSE;
570 : }
571 : else
572 : {
573 0 : *pdfMin = MIN(*pdfMin,dfValue);
574 0 : *pdfMax = MAX(*pdfMax,dfValue);
575 : }
576 : }
577 : }
578 : }
579 :
580 0 : CPLFree( pBuffer );
581 :
582 0 : return !bFirst;
583 : }
584 :
585 :
|
