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