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