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