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