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