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