Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: IRIS Reader
4 : * Purpose: All code for IRIS format Reader
5 : * Author: Roger Veciana, rveciana@gmail.com
6 : * Portions are adapted from code copyright (C) 2005-2012
7 : * Chris Veness under a CC-BY 3.0 licence
8 : *
9 : ******************************************************************************
10 : * Copyright (c) 2012, Roger Veciana <rveciana@gmail.com>
11 : * Copyright (c) 2012-2013, Even Rouault <even dot rouault at spatialys.com>
12 : *
13 : * SPDX-License-Identifier: MIT
14 : ****************************************************************************/
15 :
16 : #include "cpl_port.h"
17 : #include "gdal_frmts.h"
18 : #include "gdal_pam.h"
19 : #include "ogr_spatialref.h"
20 :
21 : #include <algorithm>
22 : #include <sstream>
23 :
24 : static double DEG2RAD = M_PI / 180.0;
25 : static double RAD2DEG = 180.0 / M_PI;
26 :
27 : /************************************************************************/
28 : /* ==================================================================== */
29 : /* IRISDataset */
30 : /* ==================================================================== */
31 : /************************************************************************/
32 :
33 : class IRISRasterBand;
34 :
35 : class IRISDataset final : public GDALPamDataset
36 : {
37 : friend class IRISRasterBand;
38 :
39 : VSILFILE *fp;
40 : GByte abyHeader[640];
41 : bool bNoDataSet;
42 : double dfNoDataValue;
43 : static const char *const aszProductNames[];
44 : static const char *const aszDataTypeCodes[];
45 : static const char *const aszDataTypes[];
46 : static const char *const aszProjections[];
47 : unsigned short nProductCode;
48 : unsigned short nDataTypeCode;
49 : unsigned char nProjectionCode;
50 : float fNyquistVelocity;
51 : mutable OGRSpatialReference m_oSRS{};
52 : mutable double adfGeoTransform[6];
53 : mutable bool bHasLoadedProjection;
54 : void LoadProjection() const;
55 : static bool GeodesicCalculation(double fLat, double fLon, double fAngle,
56 : double fDist, double fEquatorialRadius,
57 : double fPolarRadius, double fFlattening,
58 : std::pair<double, double> &oOutPair);
59 :
60 : public:
61 : IRISDataset();
62 : virtual ~IRISDataset();
63 :
64 : static GDALDataset *Open(GDALOpenInfo *);
65 : static int Identify(GDALOpenInfo *);
66 :
67 : CPLErr GetGeoTransform(double *padfTransform) override;
68 : const OGRSpatialReference *GetSpatialRef() const override;
69 : };
70 :
71 : const char *const IRISDataset::aszProductNames[] = {
72 : "", "PPI", "RHI", "CAPPI", "CROSS", "TOPS", "TRACK",
73 : "RAIN1", "RAINN", "VVP", "VIL", "SHEAR", "WARN", "CATCH",
74 : "RTI", "RAW", "MAX", "USER", "USERV", "OTHER", "STATUS",
75 : "SLINE", "WIND", "BEAM", "TEXT", "FCAST", "NDOP", "IMAGE",
76 : "COMP", "TDWR", "GAGE", "DWELL", "SRI", "BASE", "HMAX"};
77 :
78 : const char *const IRISDataset::aszDataTypeCodes[] = {
79 : "XHDR", "DBT", "dBZ", "VEL", "WIDTH", "ZDR",
80 : "ORAIN", "dBZC", "DBT2", "dBZ2", "VEL2", "WIDTH2",
81 : "ZDR2", "RAINRATE2", "KDP", "KDP2", "PHIDP", "VELC",
82 : "SQI", "RHOHV", "RHOHV2", "dBZC2", "VELC2", "SQI2",
83 : "PHIDP2", "LDRH", "LDRH2", "LDRV", "LDRV2", "FLAGS",
84 : "FLAGS2", "FLOAT32", "HEIGHT", "VIL2", "NULL", "SHEAR",
85 : "DIVERGE2", "FLIQUID2", "USER", "OTHER", "DEFORM2", "VVEL2",
86 : "HVEL2", "HDIR2", "AXDIL2", "TIME2", "RHOH", "RHOH2",
87 : "RHOV", "RHOV2", "PHIH", "PHIH2", "PHIV", "PHIV2",
88 : "USER2", "HCLASS", "HCLASS2", "ZDRC", "ZDRC2", "TEMPERATURE16",
89 : "VIR16", "DBTV8", "DBTV16", "DBZV8", "DBZV16", "SNR8",
90 : "SNR16", "ALBEDO8", "ALBEDO16", "VILD16", "TURB16"};
91 :
92 : const char *const IRISDataset::aszDataTypes[] = {
93 : "Extended Headers",
94 : "Total H power (1 byte)",
95 : "Clutter Corrected H reflectivity (1 byte)",
96 : "Velocity (1 byte)",
97 : "Width (1 byte)",
98 : "Differential reflectivity (1 byte)",
99 : "Old Rainfall rate (stored as dBZ)",
100 : "Fully corrected reflectivity (1 byte)",
101 : "Uncorrected reflectivity (2 byte)",
102 : "Corrected reflectivity (2 byte)",
103 : "Velocity (2 byte)",
104 : "Width (2 byte)",
105 : "Differential reflectivity (2 byte)",
106 : "Rainfall rate (2 byte)",
107 : "Kdp (specific differential phase)(1 byte)",
108 : "Kdp (specific differential phase)(2 byte)",
109 : "PHIdp (differential phase)(1 byte)",
110 : "Corrected Velocity (1 byte)",
111 : "SQI (1 byte)",
112 : "RhoHV(0) (1 byte)",
113 : "RhoHV(0) (2 byte)",
114 : "Fully corrected reflectivity (2 byte)",
115 : "Corrected Velocity (2 byte)",
116 : "SQI (2 byte)",
117 : "PHIdp (differential phase)(2 byte)",
118 : "LDR H to V (1 byte)",
119 : "LDR H to V (2 byte)",
120 : "LDR V to H (1 byte)",
121 : "LDR V to H (2 byte)",
122 : "Individual flag bits for each bin",
123 : "",
124 : "Test of floating format",
125 : "Height (1/10 km) (1 byte)",
126 : "Linear liquid (.001mm) (2 byte)",
127 : "Data type is not applicable",
128 : "Wind Shear (1 byte)",
129 : "Divergence (.001 10**-4) (2-byte)",
130 : "Floated liquid (2 byte)",
131 : "User type, unspecified data (1 byte)",
132 : "Unspecified data, no color legend",
133 : "Deformation (.001 10**-4) (2-byte)",
134 : "Vertical velocity (.01 m/s) (2-byte)",
135 : "Horizontal velocity (.01 m/s) (2-byte)",
136 : "Horizontal wind direction (.1 degree) (2-byte)",
137 : "Axis of Dillitation (.1 degree) (2-byte)",
138 : "Time of data (seconds) (2-byte)",
139 : "Rho H to V (1 byte)",
140 : "Rho H to V (2 byte)",
141 : "Rho V to H (1 byte)",
142 : "Rho V to H (2 byte)",
143 : "Phi H to V (1 byte)",
144 : "Phi H to V (2 byte)",
145 : "Phi V to H (1 byte)",
146 : "Phi V to H (2 byte)",
147 : "User type, unspecified data (2 byte)",
148 : "Hydrometeor class (1 byte)",
149 : "Hydrometeor class (2 byte)",
150 : "Corrected Differential reflectivity (1 byte)",
151 : "Corrected Differential reflectivity (2 byte)",
152 : "Temperature (2 byte)",
153 : "Vertically Integrated Reflectivity (2 byte)",
154 : "Total V Power (1 byte)",
155 : "Total V Power (2 byte)",
156 : "Clutter Corrected V Reflectivity (1 byte)",
157 : "Clutter Corrected V Reflectivity (2 byte)",
158 : "Signal to Noise ratio (1 byte)",
159 : "Signal to Noise ratio (2 byte)",
160 : "Albedo (1 byte)",
161 : "Albedo (2 byte)",
162 : "VIL Density (2 byte)",
163 : "Turbulence (2 byte)"};
164 :
165 : const char *const IRISDataset::aszProjections[] = {
166 : "Azimutal equidistant", "Mercator", "Polar Stereographic", "UTM",
167 : // FIXME: is it a typo here or in IRIS itself: Perspective or Prespective ?
168 : "Perspective from geosync", "Equidistant cylindrical", "Gnomonic",
169 : "Gauss conformal", "Lambert conformal conic"};
170 :
171 : /************************************************************************/
172 : /* ==================================================================== */
173 : /* IRISRasterBand */
174 : /* ==================================================================== */
175 : /************************************************************************/
176 :
177 : class IRISRasterBand final : public GDALPamRasterBand
178 : {
179 : friend class IRISDataset;
180 :
181 : unsigned char *pszRecord;
182 : bool bBufferAllocFailed;
183 :
184 : public:
185 : IRISRasterBand(IRISDataset *, int);
186 : virtual ~IRISRasterBand();
187 :
188 : virtual CPLErr IReadBlock(int, int, void *) override;
189 :
190 : virtual double GetNoDataValue(int *) override;
191 : virtual CPLErr SetNoDataValue(double) override;
192 : };
193 :
194 : /************************************************************************/
195 : /* IRISRasterBand() */
196 : /************************************************************************/
197 :
198 3 : IRISRasterBand::IRISRasterBand(IRISDataset *poDSIn, int nBandIn)
199 3 : : pszRecord(nullptr), bBufferAllocFailed(false)
200 : {
201 3 : poDS = poDSIn;
202 3 : nBand = nBandIn;
203 3 : eDataType = GDT_Float32;
204 3 : nBlockXSize = poDS->GetRasterXSize();
205 3 : nBlockYSize = 1;
206 3 : }
207 :
208 6 : IRISRasterBand::~IRISRasterBand()
209 : {
210 3 : VSIFree(pszRecord);
211 6 : }
212 :
213 : /************************************************************************/
214 : /* IReadBlock() */
215 : /************************************************************************/
216 :
217 263 : CPLErr IRISRasterBand::IReadBlock(int /* nBlockXOff */, int nBlockYOff,
218 : void *pImage)
219 :
220 : {
221 263 : IRISDataset *poGDS = static_cast<IRISDataset *>(poDS);
222 :
223 : // Every product type has its own size. TODO: Move it like dataType.
224 263 : int nDataLength = 1;
225 263 : if (poGDS->nDataTypeCode == 2)
226 263 : nDataLength = 1;
227 0 : else if (poGDS->nDataTypeCode == 8)
228 0 : nDataLength = 2;
229 0 : else if (poGDS->nDataTypeCode == 9)
230 0 : nDataLength = 2;
231 0 : else if (poGDS->nDataTypeCode == 37)
232 0 : nDataLength = 2;
233 0 : else if (poGDS->nDataTypeCode == 33)
234 0 : nDataLength = 2;
235 0 : else if (poGDS->nDataTypeCode == 32)
236 0 : nDataLength = 1;
237 :
238 : // We allocate space for storing a record:
239 263 : if (pszRecord == nullptr)
240 : {
241 2 : if (bBufferAllocFailed)
242 0 : return CE_Failure;
243 :
244 2 : pszRecord = static_cast<unsigned char *>(
245 2 : VSI_MALLOC_VERBOSE(static_cast<size_t>(nBlockXSize) * nDataLength));
246 :
247 2 : if (pszRecord == nullptr)
248 : {
249 0 : bBufferAllocFailed = true;
250 0 : return CE_Failure;
251 : }
252 : }
253 :
254 : // Prepare to read (640 is the header size in bytes) and read (the
255 : // y axis in the IRIS files in the inverse direction). The
256 : // previous bands are also added as an offset
257 :
258 263 : VSIFSeekL(poGDS->fp,
259 : 640 +
260 526 : static_cast<vsi_l_offset>(nDataLength) *
261 263 : poGDS->GetRasterXSize() * poGDS->GetRasterYSize() *
262 526 : (this->nBand - 1) +
263 526 : static_cast<vsi_l_offset>(nBlockXSize) * nDataLength *
264 263 : (poGDS->GetRasterYSize() - 1 - nBlockYOff),
265 : SEEK_SET);
266 :
267 263 : if (static_cast<int>(
268 263 : VSIFReadL(pszRecord, static_cast<size_t>(nBlockXSize) * nDataLength,
269 263 : 1, poGDS->fp)) != 1)
270 0 : return CE_Failure;
271 :
272 : // If datatype is dbZ or dBT:
273 : // See point 3.3.3 at page 3.33 of the manual.
274 263 : if (poGDS->nDataTypeCode == 2 || poGDS->nDataTypeCode == 1)
275 : {
276 68384 : for (int i = 0; i < nBlockXSize; i++)
277 : {
278 68121 : float fVal = ((*(pszRecord + i * nDataLength)) - 64.0f) / 2.0f;
279 68121 : if (fVal == 95.5f)
280 0 : fVal = -9999.0f;
281 68121 : ((float *)pImage)[i] = fVal;
282 263 : }
283 : // If datatype is dbZ2 or dBT2:
284 : // See point 3.3.4 at page 3.33 of the manual.
285 : }
286 0 : else if (poGDS->nDataTypeCode == 8 || poGDS->nDataTypeCode == 9)
287 : {
288 0 : for (int i = 0; i < nBlockXSize; i++)
289 : {
290 0 : float fVal =
291 0 : (CPL_LSBUINT16PTR(pszRecord + i * nDataLength) - 32768.0f) /
292 : 100.0f;
293 0 : if (fVal == 327.67f)
294 0 : fVal = -9999.0f;
295 0 : ((float *)pImage)[i] = fVal;
296 0 : }
297 : // Fliquid2 (Rain1 & Rainn products)
298 : // See point 3.3.11 at page 3.43 of the manual.
299 : }
300 0 : else if (poGDS->nDataTypeCode == 37)
301 : {
302 0 : for (int i = 0; i < nBlockXSize; i++)
303 : {
304 0 : const unsigned short nVal =
305 0 : CPL_LSBUINT16PTR(pszRecord + i * nDataLength);
306 0 : const unsigned short nExp = nVal >> 12;
307 0 : const unsigned short nMantissa = nVal - (nExp << 12);
308 0 : float fVal2 = 0.0f;
309 0 : if (nVal == 65535)
310 0 : fVal2 = -9999.0f;
311 0 : else if (nExp == 0)
312 0 : fVal2 = nMantissa / 1000.0f;
313 : else
314 0 : fVal2 = ((nMantissa + 4096) << (nExp - 1)) / 1000.0f;
315 0 : ((float *)pImage)[i] = fVal2;
316 : }
317 : // VIL2 (VIL products)
318 : // See point 3.3.41 at page 3.54 of the manual.
319 : }
320 0 : else if (poGDS->nDataTypeCode == 33)
321 : {
322 0 : for (int i = 0; i < nBlockXSize; i++)
323 : {
324 0 : float fVal = static_cast<float>(
325 0 : CPL_LSBUINT16PTR(pszRecord + i * nDataLength));
326 0 : if (fVal == 65535.0f)
327 0 : ((float *)pImage)[i] = -9999.0f;
328 0 : else if (fVal == 0.0f)
329 0 : ((float *)pImage)[i] = -1.0f;
330 : else
331 0 : ((float *)pImage)[i] = (fVal - 1) / 1000.0f;
332 : }
333 : // HEIGHT (TOPS products)
334 : // See point 3.3.14 at page 3.46 of the manual.
335 : }
336 0 : else if (poGDS->nDataTypeCode == 32)
337 : {
338 0 : for (int i = 0; i < nBlockXSize; i++)
339 : {
340 0 : unsigned char nVal = *(pszRecord + i * nDataLength);
341 0 : if (nVal == 255)
342 0 : ((float *)pImage)[i] = -9999.0f;
343 0 : else if (nVal == 0)
344 0 : ((float *)pImage)[i] = -1.0f;
345 : else
346 0 : ((float *)pImage)[i] = (nVal - 1.0f) / 10.0f;
347 : }
348 : // VEL (Velocity 1-Byte in PPI & others)
349 : // See point 3.3.37 at page 3.53 of the manual.
350 : }
351 0 : else if (poGDS->nDataTypeCode == 3)
352 : {
353 0 : for (int i = 0; i < nBlockXSize; i++)
354 : {
355 0 : float fVal = static_cast<float>(*(pszRecord + i * nDataLength));
356 0 : if (fVal == 0.0f)
357 0 : fVal = -9997.0f;
358 0 : else if (fVal == 1.0f)
359 0 : fVal = -9998.0f;
360 0 : else if (fVal == 255.0f)
361 0 : fVal = -9999.0f;
362 : else
363 0 : fVal = poGDS->fNyquistVelocity * (fVal - 128.0f) / 127.0f;
364 0 : ((float *)pImage)[i] = fVal;
365 : }
366 : // SHEAR (1-Byte Shear)
367 : // See point 3.3.23 at page 3.39 of the manual.
368 : }
369 0 : else if (poGDS->nDataTypeCode == 35)
370 : {
371 0 : for (int i = 0; i < nBlockXSize; i++)
372 : {
373 0 : float fVal = static_cast<float>(*(pszRecord + i * nDataLength));
374 0 : if (fVal == 0.0f)
375 0 : fVal = -9998.0f;
376 0 : else if (fVal == 255.0f)
377 0 : fVal = -9999.0f;
378 : else
379 0 : fVal = (fVal - 128.0f) * 0.2f;
380 0 : ((float *)pImage)[i] = fVal;
381 : }
382 : }
383 :
384 263 : return CE_None;
385 : }
386 :
387 : /************************************************************************/
388 : /* SetNoDataValue() */
389 : /************************************************************************/
390 :
391 3 : CPLErr IRISRasterBand::SetNoDataValue(double dfNoData)
392 :
393 : {
394 3 : IRISDataset *poGDS = static_cast<IRISDataset *>(poDS);
395 :
396 3 : poGDS->bNoDataSet = true;
397 3 : poGDS->dfNoDataValue = dfNoData;
398 :
399 3 : return CE_None;
400 : }
401 :
402 : /************************************************************************/
403 : /* GetNoDataValue() */
404 : /************************************************************************/
405 :
406 0 : double IRISRasterBand::GetNoDataValue(int *pbSuccess)
407 :
408 : {
409 0 : IRISDataset *poGDS = static_cast<IRISDataset *>(poDS);
410 :
411 0 : if (poGDS->bNoDataSet)
412 : {
413 0 : if (pbSuccess)
414 0 : *pbSuccess = TRUE;
415 :
416 0 : return poGDS->dfNoDataValue;
417 : }
418 :
419 0 : return GDALPamRasterBand::GetNoDataValue(pbSuccess);
420 : }
421 :
422 : /************************************************************************/
423 : /* ==================================================================== */
424 : /* IRISDataset */
425 : /* ==================================================================== */
426 : /************************************************************************/
427 :
428 : /************************************************************************/
429 : /* IRISDataset() */
430 : /************************************************************************/
431 :
432 3 : IRISDataset::IRISDataset()
433 : : fp(nullptr), bNoDataSet(false), dfNoDataValue(0.0), nProductCode(0),
434 : nDataTypeCode(0), nProjectionCode(0), fNyquistVelocity(0.0),
435 3 : bHasLoadedProjection(false)
436 : {
437 3 : m_oSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
438 3 : std::fill_n(abyHeader, CPL_ARRAYSIZE(abyHeader), static_cast<GByte>(0));
439 3 : adfGeoTransform[0] = 0.0;
440 3 : adfGeoTransform[1] = 1.0;
441 3 : adfGeoTransform[2] = 0.0;
442 3 : adfGeoTransform[3] = 0.0;
443 3 : adfGeoTransform[4] = 0.0;
444 3 : adfGeoTransform[5] = 1.0;
445 3 : }
446 :
447 : /************************************************************************/
448 : /* ~IRISDataset() */
449 : /************************************************************************/
450 :
451 6 : IRISDataset::~IRISDataset()
452 :
453 : {
454 3 : FlushCache(true);
455 3 : if (fp != nullptr)
456 3 : VSIFCloseL(fp);
457 6 : }
458 :
459 : /************************************************************************/
460 : /* Calculates the projection and Geotransform */
461 : /************************************************************************/
462 1 : void IRISDataset::LoadProjection() const
463 : {
464 1 : bHasLoadedProjection = true;
465 : // They give the radius in cm.
466 1 : double dfEquatorialRadius =
467 1 : CPL_LSBUINT32PTR(abyHeader + 220 + 320 + 12) / 100.0;
468 : // Point 3.2.27 pag 3-15.
469 1 : double dfInvFlattening =
470 1 : CPL_LSBUINT32PTR(abyHeader + 224 + 320 + 12) / 1000000.0;
471 1 : double dfFlattening = 0.0;
472 1 : double dfPolarRadius = 0.0;
473 :
474 1 : if (dfEquatorialRadius == 0.0)
475 : {
476 : // If Radius is 0, change to 6371000 Point 3.2.27 pag 3-15 (old IRIS
477 : // versions).
478 0 : dfEquatorialRadius = 6371000.0;
479 0 : dfPolarRadius = dfEquatorialRadius;
480 0 : dfInvFlattening = 0.0;
481 0 : dfFlattening = 0.0;
482 : }
483 : else
484 : {
485 1 : if (dfInvFlattening == 0.0)
486 : {
487 : // When inverse flattening is infinite, they use 0.
488 1 : dfFlattening = 0.0;
489 1 : dfPolarRadius = dfEquatorialRadius;
490 : }
491 : else
492 : {
493 0 : dfFlattening = 1.0 / dfInvFlattening;
494 0 : dfPolarRadius = dfEquatorialRadius * (1.0 - dfFlattening);
495 : }
496 : }
497 :
498 1 : constexpr GUInt32 knUINT32_MAX = 0xFFFFFFFFU;
499 1 : const double dfCenterLon =
500 1 : CPL_LSBUINT32PTR(abyHeader + 112 + 320 + 12) * 360.0 / knUINT32_MAX;
501 1 : const double dfCenterLat =
502 1 : CPL_LSBUINT32PTR(abyHeader + 108 + 320 + 12) * 360.0 / knUINT32_MAX;
503 :
504 1 : const double dfProjRefLon =
505 1 : CPL_LSBUINT32PTR(abyHeader + 244 + 320 + 12) * 360.0 / knUINT32_MAX;
506 1 : const double dfProjRefLat =
507 1 : CPL_LSBUINT32PTR(abyHeader + 240 + 320 + 12) * 360.0 / knUINT32_MAX;
508 :
509 1 : const double dfRadarLocX = CPL_LSBSINT32PTR(abyHeader + 112 + 12) / 1000.0;
510 1 : const double dfRadarLocY = CPL_LSBSINT32PTR(abyHeader + 116 + 12) / 1000.0;
511 :
512 1 : const double dfScaleX = CPL_LSBSINT32PTR(abyHeader + 88 + 12) / 100.0;
513 1 : const double dfScaleY = CPL_LSBSINT32PTR(abyHeader + 92 + 12) / 100.0;
514 1 : if (dfScaleX <= 0.0 || dfScaleY <= 0.0 || dfScaleX >= dfPolarRadius ||
515 : dfScaleY >= dfPolarRadius)
516 0 : return;
517 :
518 : // Mercator projection.
519 1 : if (EQUAL(aszProjections[nProjectionCode], "Mercator"))
520 : {
521 1 : std::pair<double, double> oPositionX2;
522 1 : if (!GeodesicCalculation(dfCenterLat, dfCenterLon, 90.0, dfScaleX,
523 : dfEquatorialRadius, dfPolarRadius,
524 : dfFlattening, oPositionX2))
525 0 : return;
526 1 : std::pair<double, double> oPositionY2;
527 1 : if (!GeodesicCalculation(dfCenterLat, dfCenterLon, 0.0, dfScaleY,
528 : dfEquatorialRadius, dfPolarRadius,
529 : dfFlattening, oPositionY2))
530 0 : return;
531 :
532 1 : m_oSRS.SetGeogCS("unnamed ellipse", "unknown", "unnamed",
533 : dfEquatorialRadius, dfInvFlattening, "Greenwich", 0.0,
534 : "degree", 0.0174532925199433);
535 :
536 1 : m_oSRS.SetMercator(dfProjRefLat, dfProjRefLon, 1.0, 0.0, 0.0);
537 1 : m_oSRS.SetLinearUnits("Metre", 1.0);
538 :
539 : // The center coordinates are given in LatLon on the defined
540 : // ellipsoid. Necessary to calculate geotransform.
541 :
542 2 : OGRSpatialReference oSRSLatLon;
543 1 : oSRSLatLon.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
544 1 : oSRSLatLon.SetGeogCS("unnamed ellipse", "unknown", "unnamed",
545 : dfEquatorialRadius, dfInvFlattening, "Greenwich",
546 : 0.0, "degree", 0.0174532925199433);
547 :
548 : OGRCoordinateTransformation *poTransform =
549 1 : OGRCreateCoordinateTransformation(&oSRSLatLon, &m_oSRS);
550 :
551 1 : const double dfLon2 = oPositionX2.first;
552 1 : const double dfLat2 = oPositionY2.second;
553 :
554 1 : double dfX = dfCenterLon;
555 1 : double dfY = dfCenterLat;
556 1 : if (poTransform == nullptr || !poTransform->Transform(1, &dfX, &dfY))
557 0 : CPLError(CE_Failure, CPLE_None, "Transformation Failed");
558 :
559 1 : double dfX2 = dfLon2;
560 1 : double dfY2 = dfLat2;
561 1 : if (poTransform == nullptr || !poTransform->Transform(1, &dfX2, &dfY2))
562 0 : CPLError(CE_Failure, CPLE_None, "Transformation Failed");
563 :
564 1 : adfGeoTransform[0] = dfX - (dfRadarLocX * (dfX2 - dfX));
565 1 : adfGeoTransform[1] = dfX2 - dfX;
566 1 : adfGeoTransform[2] = 0.0;
567 1 : adfGeoTransform[3] = dfY + (dfRadarLocY * (dfY2 - dfY));
568 1 : adfGeoTransform[4] = 0.0;
569 1 : adfGeoTransform[5] = -1 * (dfY2 - dfY);
570 :
571 1 : delete poTransform;
572 : }
573 0 : else if (EQUAL(aszProjections[nProjectionCode], "Azimutal equidistant"))
574 : {
575 0 : m_oSRS.SetGeogCS("unnamed ellipse", "unknown", "unnamed",
576 : dfEquatorialRadius, dfInvFlattening, "Greenwich", 0.0,
577 : "degree", 0.0174532925199433);
578 0 : m_oSRS.SetAE(dfProjRefLat, dfProjRefLon, 0.0, 0.0);
579 :
580 0 : adfGeoTransform[0] = -1 * (dfRadarLocX * dfScaleX);
581 0 : adfGeoTransform[1] = dfScaleX;
582 0 : adfGeoTransform[2] = 0.0;
583 0 : adfGeoTransform[3] = dfRadarLocY * dfScaleY;
584 0 : adfGeoTransform[4] = 0.0;
585 0 : adfGeoTransform[5] = -1 * dfScaleY;
586 : // When the projection is different from Mercator or Azimutal
587 : // equidistant, we set a standard geotransform.
588 : }
589 : else
590 : {
591 0 : adfGeoTransform[0] = -1 * (dfRadarLocX * dfScaleX);
592 0 : adfGeoTransform[1] = dfScaleX;
593 0 : adfGeoTransform[2] = 0.0;
594 0 : adfGeoTransform[3] = dfRadarLocY * dfScaleY;
595 0 : adfGeoTransform[4] = 0.0;
596 0 : adfGeoTransform[5] = -1 * dfScaleY;
597 : }
598 : }
599 :
600 : /******************************************************************************/
601 : /* The geotransform in Mercator projection must be calculated transforming */
602 : /* distance to degrees over the ellipsoid, using Vincenty's formula. */
603 : /* The following method is ported from a version for Javascript by Chris */
604 : /* Veness distributed under a CC-BY 3.0 licence, whose conditions is that the */
605 : /* following copyright notice is retained as well as the link to : */
606 : /* http://www.movable-type.co.uk/scripts/latlong-vincenty-direct.html */
607 : /******************************************************************************/
608 :
609 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
610 : * - - - - - - - - */
611 : /* Vincenty Direct Solution of Geodesics on the Ellipsoid (c) Chris Veness
612 : * 2005-2012 */
613 : /* */
614 : /* from: Vincenty direct formula - T Vincenty, "Direct and Inverse Solutions of
615 : * Geodesics on the */
616 : /* Ellipsoid with application of nested equations", Survey Review, vol
617 : * XXII no 176, 1975 */
618 : /* http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf */
619 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
620 : * - - - - - - - - */
621 :
622 2 : bool IRISDataset::GeodesicCalculation(double fLat, double fLon, double fAngle,
623 : double fDist, double fEquatorialRadius,
624 : double fPolarRadius, double fFlattening,
625 : std::pair<double, double> &oOutPair)
626 : {
627 2 : const double dfAlpha1 = DEG2RAD * fAngle;
628 2 : const double dfSinAlpha1 = sin(dfAlpha1);
629 2 : const double dfCosAlpha1 = cos(dfAlpha1);
630 :
631 2 : const double dfTanU1 = (1 - fFlattening) * tan(fLat * DEG2RAD);
632 2 : const double dfCosU1 = 1 / sqrt((1 + dfTanU1 * dfTanU1));
633 2 : const double dfSinU1 = dfTanU1 * dfCosU1;
634 :
635 2 : const double dfSigma1 = atan2(dfTanU1, dfCosAlpha1);
636 2 : const double dfSinAlpha = dfCosU1 * dfSinAlpha1;
637 2 : const double dfCosSqAlpha = 1 - dfSinAlpha * dfSinAlpha;
638 2 : const double dfUSq =
639 2 : dfCosSqAlpha *
640 2 : (fEquatorialRadius * fEquatorialRadius - fPolarRadius * fPolarRadius) /
641 2 : (fPolarRadius * fPolarRadius);
642 2 : const double dfA =
643 : 1 +
644 2 : dfUSq / 16384 * (4096 + dfUSq * (-768 + dfUSq * (320 - 175 * dfUSq)));
645 2 : const double dfB =
646 2 : dfUSq / 1024 * (256 + dfUSq * (-128 + dfUSq * (74 - 47 * dfUSq)));
647 :
648 2 : double dfSigma = fDist / (fPolarRadius * dfA);
649 2 : double dfSigmaP = 2 * M_PI;
650 :
651 2 : double dfSinSigma = 0.0;
652 2 : double dfCosSigma = 0.0;
653 2 : double dfCos2SigmaM = 0.0;
654 :
655 2 : int nIter = 0;
656 4 : while (fabs(dfSigma - dfSigmaP) > 1e-12)
657 : {
658 2 : dfCos2SigmaM = cos(2 * dfSigma1 + dfSigma);
659 2 : dfSinSigma = sin(dfSigma);
660 2 : dfCosSigma = cos(dfSigma);
661 2 : const double dfDeltaSigma =
662 2 : dfB * dfSinSigma *
663 2 : (dfCos2SigmaM +
664 2 : dfB / 4 *
665 2 : (dfCosSigma * (-1 + 2 * dfCos2SigmaM * dfCos2SigmaM) -
666 2 : dfB / 6 * dfCos2SigmaM * (-3 + 4 * dfSinSigma * dfSinSigma) *
667 2 : (-3 + 4 * dfCos2SigmaM * dfCos2SigmaM)));
668 2 : dfSigmaP = dfSigma;
669 2 : dfSigma = fDist / (fPolarRadius * dfA) + dfDeltaSigma;
670 2 : nIter++;
671 2 : if (nIter == 100)
672 0 : return false;
673 : }
674 :
675 2 : const double dfTmp =
676 2 : dfSinU1 * dfSinSigma - dfCosU1 * dfCosSigma * dfCosAlpha1;
677 2 : const double dfLat2 = atan2(
678 2 : dfSinU1 * dfCosSigma + dfCosU1 * dfSinSigma * dfCosAlpha1,
679 2 : (1 - fFlattening) * sqrt(dfSinAlpha * dfSinAlpha + dfTmp * dfTmp));
680 : const double dfLambda =
681 2 : atan2(dfSinSigma * dfSinAlpha1,
682 2 : dfCosU1 * dfCosSigma - dfSinU1 * dfSinSigma * dfCosAlpha1);
683 2 : const double dfC = fFlattening / 16 * dfCosSqAlpha *
684 2 : (4 + fFlattening * (4 - 3 * dfCosSqAlpha));
685 2 : const double dfL =
686 : dfLambda -
687 2 : (1 - dfC) * fFlattening * dfSinAlpha *
688 2 : (dfSigma +
689 2 : dfC * dfSinSigma *
690 2 : (dfCos2SigmaM +
691 2 : dfC * dfCosSigma * (-1 + 2 * dfCos2SigmaM * dfCos2SigmaM)));
692 2 : double dfLon2 = fLon * DEG2RAD + dfL;
693 :
694 2 : if (dfLon2 > M_PI)
695 0 : dfLon2 = dfLon2 - 2 * M_PI;
696 2 : if (dfLon2 < -1 * M_PI)
697 0 : dfLon2 = dfLon2 + 2 * M_PI;
698 :
699 2 : oOutPair = std::pair<double, double>(dfLon2 * RAD2DEG, dfLat2 * RAD2DEG);
700 :
701 2 : return true;
702 : }
703 :
704 : /************************************************************************/
705 : /* GetGeoTransform() */
706 : /************************************************************************/
707 :
708 1 : CPLErr IRISDataset::GetGeoTransform(double *padfTransform)
709 :
710 : {
711 1 : if (!bHasLoadedProjection)
712 0 : LoadProjection();
713 1 : memcpy(padfTransform, adfGeoTransform, sizeof(double) * 6);
714 1 : return CE_None;
715 : }
716 :
717 : /************************************************************************/
718 : /* GetSpatialRef() */
719 : /************************************************************************/
720 :
721 1 : const OGRSpatialReference *IRISDataset::GetSpatialRef() const
722 : {
723 1 : if (!bHasLoadedProjection)
724 1 : LoadProjection();
725 1 : return m_oSRS.IsEmpty() ? nullptr : &m_oSRS;
726 : }
727 :
728 : /************************************************************************/
729 : /* Identify() */
730 : /************************************************************************/
731 :
732 50638 : int IRISDataset::Identify(GDALOpenInfo *poOpenInfo)
733 :
734 : {
735 : /* -------------------------------------------------------------------- */
736 : /* Confirm that the file is an IRIS file */
737 : /* -------------------------------------------------------------------- */
738 50638 : if (poOpenInfo->nHeaderBytes < 640)
739 48796 : return FALSE;
740 :
741 1842 : const short nId1 = CPL_LSBSINT16PTR(poOpenInfo->pabyHeader);
742 1842 : const short nId2 = CPL_LSBSINT16PTR(poOpenInfo->pabyHeader + 12);
743 1842 : unsigned short nProductCode =
744 1842 : CPL_LSBUINT16PTR(poOpenInfo->pabyHeader + 12 + 12);
745 1842 : const short nYear = CPL_LSBSINT16PTR(poOpenInfo->pabyHeader + 26 + 12);
746 1842 : const short nMonth = CPL_LSBSINT16PTR(poOpenInfo->pabyHeader + 28 + 12);
747 1842 : const short nDay = CPL_LSBSINT16PTR(poOpenInfo->pabyHeader + 30 + 12);
748 :
749 : // Check if the two headers are 27 (product hdr) & 26 (product
750 : // configuration), and other metadata
751 1848 : if (!(nId1 == 27 && nId2 == 26 && nProductCode > 0 &&
752 6 : nProductCode < CPL_ARRAYSIZE(aszProductNames) && nYear >= 1900 &&
753 6 : nYear < 2100 && nMonth >= 1 && nMonth <= 12 && nDay >= 1 &&
754 : nDay <= 31))
755 1836 : return FALSE;
756 :
757 6 : return TRUE;
758 : }
759 :
760 : /************************************************************************/
761 : /* FillString() */
762 : /************************************************************************/
763 :
764 21 : static void FillString(char *szBuffer, size_t nBufferSize, void *pSrcBuffer)
765 : {
766 285 : for (size_t i = 0; i < nBufferSize - 1; i++)
767 264 : szBuffer[i] = static_cast<char *>(pSrcBuffer)[i];
768 21 : szBuffer[nBufferSize - 1] = '\0';
769 21 : }
770 :
771 : /************************************************************************/
772 : /* Open() */
773 : /************************************************************************/
774 :
775 3 : GDALDataset *IRISDataset::Open(GDALOpenInfo *poOpenInfo)
776 :
777 : {
778 3 : if (!Identify(poOpenInfo) || poOpenInfo->fpL == nullptr)
779 0 : return nullptr;
780 : /* -------------------------------------------------------------------- */
781 : /* Confirm the requested access is supported. */
782 : /* -------------------------------------------------------------------- */
783 3 : if (poOpenInfo->eAccess == GA_Update)
784 : {
785 0 : CPLError(CE_Failure, CPLE_NotSupported,
786 : "The IRIS driver does not support update access to existing"
787 : " datasets.");
788 0 : return nullptr;
789 : }
790 :
791 : /* -------------------------------------------------------------------- */
792 : /* Create a corresponding GDALDataset. */
793 : /* -------------------------------------------------------------------- */
794 3 : IRISDataset *poDS = new IRISDataset();
795 3 : poDS->fp = poOpenInfo->fpL;
796 3 : poOpenInfo->fpL = nullptr;
797 :
798 : /* -------------------------------------------------------------------- */
799 : /* Read the header. */
800 : /* -------------------------------------------------------------------- */
801 3 : VSIFReadL(poDS->abyHeader, 1, 640, poDS->fp);
802 3 : const int nXSize = CPL_LSBSINT32PTR(poDS->abyHeader + 100 + 12);
803 3 : const int nYSize = CPL_LSBSINT32PTR(poDS->abyHeader + 104 + 12);
804 3 : const int nNumBands = CPL_LSBSINT32PTR(poDS->abyHeader + 108 + 12);
805 :
806 3 : poDS->nRasterXSize = nXSize;
807 :
808 3 : poDS->nRasterYSize = nYSize;
809 3 : if (poDS->nRasterXSize <= 0 || poDS->nRasterYSize <= 0)
810 : {
811 0 : CPLError(CE_Failure, CPLE_AppDefined, "Invalid dimensions : %d x %d",
812 : poDS->nRasterXSize, poDS->nRasterYSize);
813 0 : delete poDS;
814 0 : return nullptr;
815 : }
816 :
817 3 : if (!GDALCheckBandCount(nNumBands, TRUE))
818 : {
819 0 : delete poDS;
820 0 : return nullptr;
821 : }
822 :
823 : /* -------------------------------------------------------------------- */
824 : /* Setting the Metadata */
825 : /* -------------------------------------------------------------------- */
826 : // See point 3.2.26 at page 3.12 of the manual.
827 3 : poDS->nProductCode = CPL_LSBUINT16PTR(poDS->abyHeader + 12 + 12);
828 3 : poDS->SetMetadataItem("PRODUCT_ID",
829 6 : CPLString().Printf("%d", poDS->nProductCode));
830 3 : if (poDS->nProductCode >= CPL_ARRAYSIZE(aszProductNames))
831 : {
832 0 : delete poDS;
833 0 : return nullptr;
834 : }
835 :
836 3 : poDS->SetMetadataItem("PRODUCT", aszProductNames[poDS->nProductCode]);
837 :
838 3 : poDS->nDataTypeCode = CPL_LSBUINT16PTR(poDS->abyHeader + 130 + 12);
839 3 : if (poDS->nDataTypeCode >= CPL_ARRAYSIZE(aszDataTypeCodes))
840 : {
841 0 : delete poDS;
842 0 : return nullptr;
843 : }
844 3 : poDS->SetMetadataItem("DATA_TYPE_CODE",
845 3 : aszDataTypeCodes[poDS->nDataTypeCode]);
846 :
847 3 : if (poDS->nDataTypeCode >= CPL_ARRAYSIZE(aszDataTypes))
848 : {
849 0 : delete poDS;
850 0 : return nullptr;
851 : }
852 3 : poDS->SetMetadataItem("DATA_TYPE", aszDataTypes[poDS->nDataTypeCode]);
853 :
854 3 : const unsigned short nDataTypeInputCode =
855 3 : CPL_LSBUINT16PTR(poDS->abyHeader + 144 + 12);
856 3 : if (nDataTypeInputCode >= CPL_ARRAYSIZE(aszDataTypeCodes))
857 : {
858 0 : delete poDS;
859 0 : return nullptr;
860 : }
861 3 : poDS->SetMetadataItem("DATA_TYPE_INPUT_CODE",
862 3 : aszDataTypeCodes[nDataTypeInputCode]);
863 :
864 3 : const unsigned short nDataTypeInput =
865 3 : CPL_LSBUINT16PTR(poDS->abyHeader + 144 + 12);
866 3 : if (nDataTypeInput >= CPL_ARRAYSIZE(aszDataTypes))
867 : {
868 0 : delete poDS;
869 0 : return nullptr;
870 : }
871 3 : poDS->SetMetadataItem("DATA_TYPE_INPUT", aszDataTypes[nDataTypeInput]);
872 :
873 3 : poDS->nProjectionCode =
874 3 : *static_cast<unsigned char *>(poDS->abyHeader + 146 + 12);
875 3 : if (poDS->nProjectionCode >= CPL_ARRAYSIZE(aszProjections))
876 : {
877 0 : delete poDS;
878 0 : return nullptr;
879 : }
880 :
881 : // Times.
882 : {
883 3 : const int nSeconds = CPL_LSBSINT32PTR(poDS->abyHeader + 20 + 12);
884 :
885 3 : const int nHour = (nSeconds - (nSeconds % 3600)) / 3600;
886 3 : const int nMinute =
887 3 : ((nSeconds - nHour * 3600) - (nSeconds - nHour * 3600) % 60) / 60;
888 3 : const int nSecond = nSeconds - nHour * 3600 - nMinute * 60;
889 :
890 3 : const short nYear = CPL_LSBSINT16PTR(poDS->abyHeader + 26 + 12);
891 3 : const short nMonth = CPL_LSBSINT16PTR(poDS->abyHeader + 28 + 12);
892 3 : const short nDay = CPL_LSBSINT16PTR(poDS->abyHeader + 30 + 12);
893 :
894 3 : poDS->SetMetadataItem("TIME_PRODUCT_GENERATED",
895 3 : CPLString().Printf("%d-%02d-%02d %02d:%02d:%02d",
896 : nYear, nMonth, nDay, nHour,
897 3 : nMinute, nSecond));
898 : }
899 :
900 : {
901 3 : const int nSeconds = CPL_LSBSINT32PTR(poDS->abyHeader + 32 + 12);
902 :
903 3 : const int nHour = (nSeconds - (nSeconds % 3600)) / 3600;
904 3 : const int nMinute =
905 3 : ((nSeconds - nHour * 3600) - (nSeconds - nHour * 3600) % 60) / 60;
906 3 : const int nSecond = nSeconds - nHour * 3600 - nMinute * 60;
907 :
908 3 : const short nYear = CPL_LSBSINT16PTR(poDS->abyHeader + 26 + 12);
909 3 : const short nMonth = CPL_LSBSINT16PTR(poDS->abyHeader + 28 + 12);
910 3 : const short nDay = CPL_LSBSINT16PTR(poDS->abyHeader + 30 + 12);
911 :
912 3 : poDS->SetMetadataItem("TIME_INPUT_INGEST_SWEEP",
913 3 : CPLString().Printf("%d-%02d-%02d %02d:%02d:%02d",
914 : nYear, nMonth, nDay, nHour,
915 3 : nMinute, nSecond));
916 : }
917 :
918 : // Site and task information.
919 :
920 3 : char szSiteName[17] = {}; // Must have one extra char for string end.
921 3 : char szVersionName[9] = {};
922 :
923 3 : FillString(szSiteName, sizeof(szSiteName), poDS->abyHeader + 320 + 12);
924 3 : FillString(szVersionName, sizeof(szVersionName),
925 3 : poDS->abyHeader + 16 + 320 + 12);
926 3 : poDS->SetMetadataItem("PRODUCT_SITE_NAME", szSiteName);
927 3 : poDS->SetMetadataItem("PRODUCT_SITE_IRIS_VERSION", szVersionName);
928 :
929 3 : FillString(szSiteName, sizeof(szSiteName), poDS->abyHeader + 90 + 320 + 12);
930 3 : FillString(szVersionName, sizeof(szVersionName),
931 3 : poDS->abyHeader + 24 + 320 + 12);
932 3 : poDS->SetMetadataItem("INGEST_SITE_NAME", szSiteName);
933 3 : poDS->SetMetadataItem("INGEST_SITE_IRIS_VERSION", szVersionName);
934 :
935 3 : FillString(szSiteName, sizeof(szSiteName), poDS->abyHeader + 74 + 320 + 12);
936 3 : poDS->SetMetadataItem("INGEST_HARDWARE_NAME", szSiteName);
937 :
938 3 : char szConfigFile[13] = {};
939 3 : FillString(szConfigFile, sizeof(szConfigFile), poDS->abyHeader + 62 + 12);
940 3 : poDS->SetMetadataItem("PRODUCT_CONFIGURATION_NAME", szConfigFile);
941 :
942 3 : char szTaskName[13] = {};
943 3 : FillString(szTaskName, sizeof(szTaskName), poDS->abyHeader + 74 + 12);
944 3 : poDS->SetMetadataItem("TASK_NAME", szTaskName);
945 :
946 3 : const short nRadarHeight =
947 3 : CPL_LSBSINT16PTR(poDS->abyHeader + 284 + 320 + 12);
948 3 : poDS->SetMetadataItem("RADAR_HEIGHT",
949 6 : CPLString().Printf("%d m", nRadarHeight));
950 : // Ground height over the sea level.
951 3 : const short nGroundHeight =
952 3 : CPL_LSBSINT16PTR(poDS->abyHeader + 118 + 320 + 12);
953 3 : poDS->SetMetadataItem(
954 : "GROUND_HEIGHT",
955 6 : CPLString().Printf("%d m", nRadarHeight - nGroundHeight));
956 :
957 3 : unsigned short nFlags = CPL_LSBUINT16PTR(poDS->abyHeader + 86 + 12);
958 : // Get eleventh bit.
959 3 : nFlags = nFlags << 4;
960 3 : nFlags = nFlags >> 15;
961 3 : if (nFlags == 1)
962 : {
963 1 : poDS->SetMetadataItem("COMPOSITED_PRODUCT", "YES");
964 1 : const unsigned int compositedMask =
965 1 : CPL_LSBUINT32PTR(poDS->abyHeader + 232 + 320 + 12);
966 1 : poDS->SetMetadataItem("COMPOSITED_PRODUCT_MASK",
967 2 : CPLString().Printf("0x%08x", compositedMask));
968 : }
969 : else
970 : {
971 2 : poDS->SetMetadataItem("COMPOSITED_PRODUCT", "NO");
972 : }
973 :
974 : // Wave values.
975 3 : poDS->SetMetadataItem(
976 3 : "PRF", CPLString().Printf("%d Hz", CPL_LSBSINT32PTR(poDS->abyHeader +
977 3 : 120 + 320 + 12)));
978 3 : poDS->SetMetadataItem(
979 : "WAVELENGTH",
980 3 : CPLString().Printf("%4.2f cm",
981 3 : CPL_LSBSINT32PTR(poDS->abyHeader + 148 + 320 + 12) /
982 3 : 100.0f));
983 3 : const unsigned short nPolarizationType =
984 3 : CPL_LSBUINT16PTR(poDS->abyHeader + 172 + 320 + 12);
985 :
986 : // See section 3.3.37 & 3.2.54.
987 3 : float fNyquist = (CPL_LSBSINT32PTR(poDS->abyHeader + 120 + 320 + 12)) *
988 3 : (static_cast<float>(
989 3 : CPL_LSBSINT32PTR(poDS->abyHeader + 148 + 320 + 12)) /
990 : 10000.0f) /
991 : 4.0f;
992 3 : if (nPolarizationType == 1)
993 0 : fNyquist = fNyquist * 2.0f;
994 3 : else if (nPolarizationType == 2)
995 0 : fNyquist = fNyquist * 3.0f;
996 3 : else if (nPolarizationType == 3)
997 0 : fNyquist = fNyquist * 4.0f;
998 3 : poDS->fNyquistVelocity = fNyquist;
999 3 : poDS->SetMetadataItem("NYQUIST_VELOCITY",
1000 6 : CPLString().Printf("%.2f m/s", fNyquist));
1001 :
1002 : // Product dependent metadata (stored in 80 bytes from 162 bytes
1003 : // at the product header) See point 3.2.30 at page 3.19 of the
1004 : // manual.
1005 : // See point 3.2.25 at page 3.12 of the manual.
1006 3 : if (EQUAL(aszProductNames[poDS->nProductCode], "PPI"))
1007 : {
1008 : // Degrees = 360 * (Binary Angle)*2^N
1009 2 : const float fElevation =
1010 2 : CPL_LSBSINT16PTR(poDS->abyHeader + 164 + 12) * 360.0f / 65536.0f;
1011 :
1012 2 : poDS->SetMetadataItem("PPI_ELEVATION_ANGLE",
1013 4 : CPLString().Printf("%f", fElevation));
1014 2 : if (EQUAL(aszDataTypeCodes[poDS->nDataTypeCode], "dBZ"))
1015 2 : poDS->SetMetadataItem("DATA_TYPE_UNITS", "dBZ");
1016 : else
1017 0 : poDS->SetMetadataItem("DATA_TYPE_UNITS", "m/s");
1018 : // See point 3.2.2 at page 3.2 of the manual.
1019 : }
1020 1 : else if (EQUAL(aszProductNames[poDS->nProductCode], "CAPPI"))
1021 : {
1022 1 : const float fElevation =
1023 1 : CPL_LSBSINT32PTR(poDS->abyHeader + 4 + 164 + 12) / 100.0f;
1024 1 : poDS->SetMetadataItem("CAPPI_BOTTOM_HEIGHT",
1025 2 : CPLString().Printf("%.1f m", fElevation));
1026 1 : const float fAzimuthSmoothingForShear =
1027 1 : CPL_LSBUINT16PTR(poDS->abyHeader + 10 + 164 + 12) * 360.0f /
1028 : 65536.0f;
1029 1 : poDS->SetMetadataItem(
1030 : "AZIMUTH_SMOOTHING_FOR_SHEAR",
1031 2 : CPLString().Printf("%.1f", fAzimuthSmoothingForShear));
1032 1 : const unsigned int nMaxAgeVVPCorrection =
1033 1 : CPL_LSBUINT32PTR(poDS->abyHeader + 24 + 164 + 12);
1034 1 : poDS->SetMetadataItem("MAX_AGE_FOR_SHEAR_VVP_CORRECTION",
1035 2 : CPLString().Printf("%d s", nMaxAgeVVPCorrection));
1036 1 : if (EQUAL(aszDataTypeCodes[poDS->nDataTypeCode], "dBZ"))
1037 1 : poDS->SetMetadataItem("DATA_TYPE_UNITS", "dBZ");
1038 : else
1039 0 : poDS->SetMetadataItem("DATA_TYPE_UNITS", "m/s");
1040 : // See point 3.2.32 at page 3.19 of the manual.
1041 : }
1042 0 : else if (EQUAL(aszProductNames[poDS->nProductCode], "RAIN1") ||
1043 0 : EQUAL(aszProductNames[poDS->nProductCode], "RAINN"))
1044 : {
1045 0 : const short nNumProducts =
1046 0 : CPL_LSBSINT16PTR(poDS->abyHeader + 170 + 320 + 12);
1047 0 : poDS->SetMetadataItem("NUM_FILES_USED",
1048 0 : CPLString().Printf("%d", nNumProducts));
1049 :
1050 0 : const float fMinZAcum =
1051 0 : (CPL_LSBUINT32PTR(poDS->abyHeader + 164 + 12) - 32768.0f) /
1052 : 10000.0f;
1053 0 : poDS->SetMetadataItem("MINIMUM_Z_TO_ACCUMULATE",
1054 0 : CPLString().Printf("%f", fMinZAcum));
1055 :
1056 0 : const unsigned short nSecondsOfAccumulation =
1057 0 : CPL_LSBUINT16PTR(poDS->abyHeader + 6 + 164 + 12);
1058 0 : poDS->SetMetadataItem(
1059 : "SECONDS_OF_ACCUMULATION",
1060 0 : CPLString().Printf("%d s", nSecondsOfAccumulation));
1061 :
1062 0 : const unsigned int nSpanInputFiles =
1063 0 : CPL_LSBUINT32PTR(poDS->abyHeader + 24 + 164 + 12);
1064 0 : poDS->SetMetadataItem("SPAN_OF_INPUT_FILES",
1065 0 : CPLString().Printf("%d s", nSpanInputFiles));
1066 0 : poDS->SetMetadataItem("DATA_TYPE_UNITS", "mm");
1067 :
1068 0 : char szInputProductName[13] = "";
1069 0 : for (int k = 0; k < 12; k++)
1070 0 : szInputProductName[k] =
1071 0 : *reinterpret_cast<char *>(poDS->abyHeader + k + 12 + 164 + 12);
1072 :
1073 0 : poDS->SetMetadataItem("INPUT_PRODUCT_NAME",
1074 0 : CPLString().Printf("%s", szInputProductName));
1075 :
1076 0 : if (EQUAL(aszProductNames[poDS->nProductCode], "RAINN"))
1077 0 : poDS->SetMetadataItem(
1078 : "NUM_HOURS_ACCUMULATE",
1079 0 : CPLString().Printf(
1080 0 : "%d", CPL_LSBUINT16PTR(poDS->abyHeader + 10 + 164 + 12)));
1081 :
1082 : // See point 3.2.73 at page 3.36 of the manual.
1083 : }
1084 0 : else if (EQUAL(aszProductNames[poDS->nProductCode], "VIL"))
1085 : {
1086 0 : const float fBottomHeightInterval =
1087 0 : CPL_LSBSINT32PTR(poDS->abyHeader + 4 + 164 + 12) / 100.0f;
1088 : // TYPO in metadata key: FIXME ?
1089 0 : poDS->SetMetadataItem(
1090 : "BOTTOM_OF_HEIGTH_INTERVAL",
1091 0 : CPLString().Printf("%.1f m", fBottomHeightInterval));
1092 0 : const float fTopHeightInterval =
1093 0 : CPL_LSBSINT32PTR(poDS->abyHeader + 8 + 164 + 12) / 100.0f;
1094 : // TYPO in metadata key: FIXME ?
1095 0 : poDS->SetMetadataItem("TOP_OF_HEIGTH_INTERVAL",
1096 0 : CPLString().Printf("%.1f m", fTopHeightInterval));
1097 0 : poDS->SetMetadataItem("VIL_DENSITY_NOT_AVAILABLE_VALUE", "-1");
1098 0 : poDS->SetMetadataItem("DATA_TYPE_UNITS", "mm");
1099 : // See point 3.2.68 at page 3.36 of the manual
1100 : }
1101 0 : else if (EQUAL(aszProductNames[poDS->nProductCode], "TOPS"))
1102 : {
1103 0 : const float fZThreshold =
1104 0 : CPL_LSBSINT16PTR(poDS->abyHeader + 4 + 164 + 12) / 16.0f;
1105 0 : poDS->SetMetadataItem("Z_THRESHOLD",
1106 0 : CPLString().Printf("%.1f dBZ", fZThreshold));
1107 0 : poDS->SetMetadataItem("ECHO_TOPS_NOT_AVAILABLE_VALUE", "-1");
1108 0 : poDS->SetMetadataItem("DATA_TYPE_UNITS", "km");
1109 : // See point 3.2.20 at page 3.10 of the manual.
1110 : }
1111 0 : else if (EQUAL(aszProductNames[poDS->nProductCode], "MAX"))
1112 : {
1113 0 : const float fBottomInterval =
1114 0 : CPL_LSBSINT32PTR(poDS->abyHeader + 4 + 164 + 12) / 100.0f;
1115 0 : poDS->SetMetadataItem("BOTTOM_OF_INTERVAL",
1116 0 : CPLString().Printf("%.1f m", fBottomInterval));
1117 0 : const float fTopInterval =
1118 0 : CPL_LSBSINT32PTR(poDS->abyHeader + 8 + 164 + 12) / 100.0f;
1119 0 : poDS->SetMetadataItem("TOP_OF_INTERVAL",
1120 0 : CPLString().Printf("%.1f m", fTopInterval));
1121 0 : const int nNumPixelsSidePanels =
1122 0 : CPL_LSBSINT32PTR(poDS->abyHeader + 12 + 164 + 12);
1123 0 : poDS->SetMetadataItem("NUM_PIXELS_SIDE_PANELS",
1124 0 : CPLString().Printf("%d", nNumPixelsSidePanels));
1125 0 : const short nHorizontalSmootherSidePanels =
1126 0 : CPL_LSBSINT16PTR(poDS->abyHeader + 16 + 164 + 12);
1127 0 : poDS->SetMetadataItem(
1128 : "HORIZONTAL_SMOOTHER_SIDE_PANELS",
1129 0 : CPLString().Printf("%d", nHorizontalSmootherSidePanels));
1130 0 : const short nVerticalSmootherSidePanels =
1131 0 : CPL_LSBSINT16PTR(poDS->abyHeader + 18 + 164 + 12);
1132 0 : poDS->SetMetadataItem(
1133 : "VERTICAL_SMOOTHER_SIDE_PANELS",
1134 0 : CPLString().Printf("%d", nVerticalSmootherSidePanels));
1135 : }
1136 :
1137 : /* -------------------------------------------------------------------- */
1138 : /* Create band information objects. */
1139 : /* -------------------------------------------------------------------- */
1140 : // coverity[tainted_data]
1141 6 : for (int iBandNum = 1; iBandNum <= nNumBands; iBandNum++)
1142 : {
1143 3 : poDS->SetBand(iBandNum, new IRISRasterBand(poDS, iBandNum));
1144 :
1145 3 : poDS->GetRasterBand(iBandNum)->SetNoDataValue(-9999);
1146 : // Calculating the band height to include it in the band metadata. Only
1147 : // for the CAPPI product.
1148 3 : if (EQUAL(aszProductNames[poDS->nProductCode], "CAPPI"))
1149 : {
1150 1 : const float fScaleZ =
1151 1 : CPL_LSBSINT32PTR(poDS->abyHeader + 96 + 12) / 100.0f;
1152 1 : const float fOffset =
1153 1 : CPL_LSBSINT32PTR(poDS->abyHeader + 4 + 164 + 12) / 100.0f;
1154 :
1155 2 : poDS->GetRasterBand(iBandNum)->SetMetadataItem(
1156 1 : "height", CPLString().Printf(
1157 1 : "%.0f m", fOffset + fScaleZ * (iBandNum - 1)));
1158 : }
1159 : }
1160 : /* -------------------------------------------------------------------- */
1161 : /* Initialize any PAM information. */
1162 : /* -------------------------------------------------------------------- */
1163 3 : poDS->SetDescription(poOpenInfo->pszFilename);
1164 3 : poDS->TryLoadXML();
1165 :
1166 : /* -------------------------------------------------------------------- */
1167 : /* Check for overviews. */
1168 : /* -------------------------------------------------------------------- */
1169 3 : poDS->oOvManager.Initialize(poDS, poOpenInfo->pszFilename);
1170 :
1171 3 : return poDS;
1172 : }
1173 :
1174 : /************************************************************************/
1175 : /* GDALRegister_IRIS() */
1176 : /************************************************************************/
1177 :
1178 1595 : void GDALRegister_IRIS()
1179 :
1180 : {
1181 1595 : if (GDALGetDriverByName("IRIS") != nullptr)
1182 302 : return;
1183 :
1184 1293 : GDALDriver *poDriver = new GDALDriver();
1185 :
1186 1293 : poDriver->SetDescription("IRIS");
1187 1293 : poDriver->SetMetadataItem(GDAL_DCAP_RASTER, "YES");
1188 1293 : poDriver->SetMetadataItem(GDAL_DMD_LONGNAME,
1189 1293 : "IRIS data (.PPI, .CAPPi etc)");
1190 1293 : poDriver->SetMetadataItem(GDAL_DMD_HELPTOPIC, "drivers/raster/iris.html");
1191 1293 : poDriver->SetMetadataItem(GDAL_DMD_EXTENSION, "ppi");
1192 1293 : poDriver->SetMetadataItem(GDAL_DCAP_VIRTUALIO, "YES");
1193 :
1194 1293 : poDriver->pfnOpen = IRISDataset::Open;
1195 1293 : poDriver->pfnIdentify = IRISDataset::Identify;
1196 :
1197 1293 : GetGDALDriverManager()->RegisterDriver(poDriver);
1198 : }
|
