1 : /******************************************************************************
2 : * $Id: irisdataset.cpp 25348 2012-12-26 11:37:51Z rouault $
3 : *
4 : * Project: IRIS Reader
5 : * Purpose: All code for IRIS format Reader
6 : * Author: Roger Veciana, rveciana@gmail.com
7 : * Portions are adapted from code copyright (C) 2005-2012
8 : * Chris Veness under a CC-BY 3.0 licence
9 : *
10 : ******************************************************************************
11 : * Copyright (c) 2012, Roger Veciana <rveciana@gmail.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 : #ifndef DEG2RAD
33 : # define DEG2RAD (M_PI/180.0)
34 : #endif
35 :
36 : #ifndef RAD2DEG
37 : # define RAD2DEG (180.0/M_PI)
38 : #endif
39 :
40 : #include "gdal_pam.h"
41 : #include "ogr_spatialref.h"
42 : #include <sstream>
43 :
44 :
45 : CPL_CVSID("$Id: irisdataset.cpp 25348 2012-12-26 11:37:51Z rouault $");
46 :
47 : CPL_C_START
48 : void GDALRegister_IRIS(void);
49 : CPL_C_END
50 :
51 : #define ARRAY_ELEMENT_COUNT(x) ((sizeof(x))/sizeof(x[0]))
52 :
53 : /************************************************************************/
54 : /* ==================================================================== */
55 : /* IRISDataset */
56 : /* ==================================================================== */
57 : /************************************************************************/
58 :
59 : class IRISRasterBand;
60 :
61 : class IRISDataset : public GDALPamDataset
62 : {
63 : friend class IRISRasterBand;
64 :
65 : VSILFILE *fp;
66 : GByte abyHeader[640];
67 : int bNoDataSet;
68 : double dfNoDataValue;
69 : static const char* const aszProductNames[];
70 : static const char* const aszDataTypeCodes[];
71 : static const char* const aszDataTypes[];
72 : static const char* const aszProjections[];
73 : unsigned short nProductCode;
74 : unsigned short nDataTypeCode;
75 : unsigned char nProjectionCode;
76 : float fNyquistVelocity;
77 : char* pszSRS_WKT;
78 : double adfGeoTransform[6];
79 : int bHasLoadedProjection;
80 : void LoadProjection();
81 : std::pair <double,double> GeodesicCalculation(float fLat, float fLon, float fAngle, float fDist, float fEquatorialRadius, float fPolarRadius, float fFlattening);
82 : public:
83 : IRISDataset();
84 : ~IRISDataset();
85 :
86 : static GDALDataset *Open( GDALOpenInfo * );
87 : static int Identify( GDALOpenInfo * );
88 :
89 : CPLErr GetGeoTransform( double * padfTransform );
90 : const char *GetProjectionRef();
91 :
92 : };
93 :
94 : const char* const IRISDataset::aszProductNames[]= {
95 : "", "PPI", "RHI", "CAPPI", "CROSS", "TOPS", "TRACK", "RAIN1", "RAINN",
96 : "VVP", "VIL", "SHEAR", "WARN", "CATCH", "RTI", "RAW", "MAX", "USER",
97 : "USERV", "OTHER", "STATUS", "SLINE", "WIND", "BEAM", "TEXT", "FCAST",
98 : "NDOP", "IMAGE", "COMP", "TDWR", "GAGE", "DWELL", "SRI", "BASE", "HMAX"};
99 :
100 : const char* const IRISDataset::aszDataTypeCodes[]={
101 : "XHDR", "DBT" ,"dBZ", "VEL", "WIDTH", "ZDR", "ORAIN", "dBZC", "DBT2",
102 : "dBZ2", "VEL2", "WIDTH2", "ZDR2", "RAINRATE2", "KDP", "KDP2", "PHIDP",
103 : "VELC", "SQI", "RHOHV", "RHOHV2", "dBZC2", "VELC2", "SQI2", "PHIDP2",
104 : "LDRH", "LDRH2", "LDRV", "LDRV2", "FLAGS", "FLAGS2", "FLOAT32", "HEIGHT",
105 : "VIL2", "NULL", "SHEAR", "DIVERGE2", "FLIQUID2", "USER", "OTHER", "DEFORM2",
106 : "VVEL2", "HVEL2", "HDIR2", "AXDIL2", "TIME2", "RHOH", "RHOH2", "RHOV",
107 : "RHOV2", "PHIH", "PHIH2", "PHIV", "PHIV2", "USER2", "HCLASS", "HCLASS2",
108 : "ZDRC", "ZDRC2", "TEMPERATURE16", "VIR16", "DBTV8", "DBTV16", "DBZV8",
109 : "DBZV16", "SNR8", "SNR16", "ALBEDO8", "ALBEDO16", "VILD16", "TURB16"};
110 : const char* const IRISDataset::aszDataTypes[]={
111 : "Extended Headers","Total H power (1 byte)","Clutter Corrected H reflectivity (1 byte)",
112 : "Velocity (1 byte)","Width (1 byte)","Differential reflectivity (1 byte)",
113 : "Old Rainfall rate (stored as dBZ)","Fully corrected reflectivity (1 byte)",
114 : "Uncorrected reflectivity (2 byte)","Corrected reflectivity (2 byte)",
115 : "Velocity (2 byte)","Width (2 byte)","Differential reflectivity (2 byte)",
116 : "Rainfall rate (2 byte)","Kdp (specific differential phase)(1 byte)",
117 : "Kdp (specific differential phase)(2 byte)","PHIdp (differential phase)(1 byte)",
118 : "Corrected Velocity (1 byte)","SQI (1 byte)","RhoHV(0) (1 byte)","RhoHV(0) (2 byte)",
119 : "Fully corrected reflectivity (2 byte)","Corrected Velocity (2 byte)","SQI (2 byte)",
120 : "PHIdp (differential phase)(2 byte)","LDR H to V (1 byte)","LDR H to V (2 byte)",
121 : "LDR V to H (1 byte)","LDR V to H (2 byte)","Individual flag bits for each bin","",
122 : "Test of floating format", "Height (1/10 km) (1 byte)", "Linear liquid (.001mm) (2 byte)",
123 : "Data type is not applicable", "Wind Shear (1 byte)", "Divergence (.001 10**-4) (2-byte)",
124 : "Floated liquid (2 byte)", "User type, unspecified data (1 byte)",
125 : "Unspecified data, no color legend", "Deformation (.001 10**-4) (2-byte)",
126 : "Vertical velocity (.01 m/s) (2-byte)", "Horizontal velocity (.01 m/s) (2-byte)",
127 : "Horizontal wind direction (.1 degree) (2-byte)", "Axis of Dillitation (.1 degree) (2-byte)",
128 : "Time of data (seconds) (2-byte)", "Rho H to V (1 byte)", "Rho H to V (2 byte)",
129 : "Rho V to H (1 byte)", "Rho V to H (2 byte)", "Phi H to V (1 byte)", "Phi H to V (2 byte)",
130 : "Phi V to H (1 byte)", "Phi V to H (2 byte)", "User type, unspecified data (2 byte)",
131 : "Hydrometeor class (1 byte)", "Hydrometeor class (2 byte)", "Corrected Differential reflectivity (1 byte)",
132 : "Corrected Differential reflectivity (2 byte)", "Temperature (2 byte)",
133 : "Vertically Integrated Reflectivity (2 byte)", "Total V Power (1 byte)", "Total V Power (2 byte)",
134 : "Clutter Corrected V Reflectivity (1 byte)", "Clutter Corrected V Reflectivity (2 byte)",
135 : "Signal to Noise ratio (1 byte)", "Signal to Noise ratio (2 byte)", "Albedo (1 byte)",
136 : "Albedo (2 byte)", "VIL Density (2 byte)", "Turbulence (2 byte)"};
137 : const char* const IRISDataset::aszProjections[]={
138 : "Azimutal equidistant","Mercator","Polar Stereographic","UTM",
139 : "Prespective from geosync","Equidistant cylindrical","Gnomonic",
140 : "Gauss conformal","Lambert conformal conic"};
141 :
142 : /************************************************************************/
143 : /* ==================================================================== */
144 : /* IRISRasterBand */
145 : /* ==================================================================== */
146 : /************************************************************************/
147 :
148 : class IRISRasterBand : public GDALPamRasterBand
149 : {
150 : friend class IRISDataset;
151 :
152 :
153 : unsigned char* pszRecord;
154 : int bBufferAllocFailed;
155 :
156 : public:
157 : IRISRasterBand( IRISDataset *, int );
158 : ~IRISRasterBand();
159 :
160 : virtual CPLErr IReadBlock( int, int, void * );
161 :
162 : virtual double GetNoDataValue( int * );
163 : virtual CPLErr SetNoDataValue( double );
164 : };
165 :
166 :
167 : /************************************************************************/
168 : /* IRISRasterBand() */
169 : /************************************************************************/
170 :
171 2 : IRISRasterBand::IRISRasterBand( IRISDataset *poDS, int nBand )
172 : {
173 2 : this->poDS = poDS;
174 2 : this->nBand = nBand;
175 2 : eDataType = GDT_Float32;
176 2 : nBlockXSize = poDS->GetRasterXSize();
177 2 : nBlockYSize = 1;
178 2 : pszRecord = NULL;
179 2 : bBufferAllocFailed = FALSE;
180 2 : }
181 :
182 2 : IRISRasterBand::~IRISRasterBand()
183 : {
184 2 : VSIFree(pszRecord);
185 2 : }
186 :
187 : /************************************************************************/
188 : /* IReadBlock() */
189 : /************************************************************************/
190 :
191 263 : CPLErr IRISRasterBand::IReadBlock( int nBlockXOff, int nBlockYOff,
192 : void * pImage )
193 :
194 : {
195 263 : IRISDataset *poGDS = (IRISDataset *) poDS;
196 :
197 : //printf("hola %d %s\n",poGDS->dataTypeCode,poGDS->aszDataTypeCodes[poGDS->dataTypeCode]);
198 : //Every product type has it's own size. TODO: Move it like dataType
199 263 : int nDataLength = 1;
200 263 : if(poGDS->nDataTypeCode == 2){nDataLength=1;}
201 0 : else if(poGDS->nDataTypeCode == 37){nDataLength=2;}
202 0 : else if(poGDS->nDataTypeCode == 33){nDataLength=2;}
203 0 : else if(poGDS->nDataTypeCode == 32){nDataLength=1;}
204 :
205 : int i;
206 : //We allocate space for storing a record:
207 263 : if (pszRecord == NULL)
208 : {
209 2 : if (bBufferAllocFailed)
210 0 : return CE_Failure;
211 :
212 2 : pszRecord = (unsigned char *) VSIMalloc(nBlockXSize*nDataLength);
213 :
214 2 : if (pszRecord == NULL)
215 : {
216 : CPLError(CE_Failure, CPLE_OutOfMemory,
217 0 : "Cannot allocate scanline buffer");
218 0 : bBufferAllocFailed = TRUE;
219 0 : return CE_Failure;
220 : }
221 : }
222 :
223 : //Prepare to read (640 is the header size in bytes) and read (the y axis in the IRIS files in the inverse direction)
224 : //The previous bands are also added as an offset
225 :
226 : VSIFSeekL( poGDS->fp, 640 + (vsi_l_offset)nDataLength*poGDS->GetRasterXSize()*poGDS->GetRasterYSize()*(this->nBand-1) +
227 263 : (vsi_l_offset)nBlockXSize*nDataLength*(poGDS->GetRasterYSize()-1-nBlockYOff), SEEK_SET );
228 :
229 263 : if( (int)VSIFReadL( pszRecord, nBlockXSize*nDataLength, 1, poGDS->fp ) != 1 )
230 0 : return CE_Failure;
231 :
232 : //If datatype is dbZ:
233 : //See point 3.3.5 at page 3.42 of the manual
234 263 : if(poGDS->nDataTypeCode == 2){
235 : float fVal;
236 68384 : for (i=0;i<nBlockXSize;i++){
237 68121 : fVal = (((float) *(pszRecord+i*nDataLength)) -64)/2.0;
238 68121 : if (fVal == 95.5)
239 0 : fVal = -9999;
240 68121 : ((float *) pImage)[i] = fVal;
241 : }
242 : //Fliquid2 (Rain1 & Rainn products)
243 : //See point 3.3.11 at page 3.43 of the manual
244 0 : } else if(poGDS->nDataTypeCode == 37){
245 : unsigned short nVal, nExp, nMantissa;
246 0 : float fVal2=0;
247 0 : for (i=0;i<nBlockXSize;i++){
248 0 : nVal = CPL_LSBUINT16PTR(pszRecord+i*nDataLength);
249 0 : nExp = nVal>>12;
250 0 : nMantissa = nVal - (nExp<<12);
251 0 : if (nVal == 65535)
252 0 : fVal2 = -9999;
253 0 : else if (nExp == 0)
254 0 : fVal2 = (float) nMantissa / 1000.0;
255 : else
256 0 : fVal2 = (float)((nMantissa+4096)<<(nExp-1))/1000.0;
257 0 : ((float *) pImage)[i] = fVal2;
258 : }
259 : //VIL2 (VIL products)
260 : //See point 3.3.41 at page 3.54 of the manual
261 0 : } else if(poGDS->nDataTypeCode == 33){
262 : float fVal;
263 0 : for (i=0;i<nBlockXSize;i++){
264 0 : fVal = (float) CPL_LSBUINT16PTR(pszRecord+i*nDataLength);
265 0 : if (fVal == 65535)
266 0 : ((float *) pImage)[i] = -9999;
267 0 : else if (fVal == 0)
268 0 : ((float *) pImage)[i] = -1;
269 : else
270 0 : ((float *) pImage)[i] = (fVal-1)/1000;
271 : }
272 : //HEIGTH (TOPS products)
273 : //See point 3.3.14 at page 3.46 of the manual
274 0 : } else if(poGDS->nDataTypeCode == 32){
275 : unsigned char nVal;
276 0 : for (i=0;i<nBlockXSize;i++){
277 0 : nVal = *(pszRecord+i*nDataLength) ;
278 0 : if (nVal == 255)
279 0 : ((float *) pImage)[i] = -9999;
280 0 : else if (nVal == 0)
281 0 : ((float *) pImage)[i] = -1;
282 : else
283 0 : ((float *) pImage)[i] = ((float) nVal - 1) / 10;
284 : }
285 : //VEL (Velocity 1-Byte in PPI & others)
286 : //See point 3.3.37 at page 3.53 of the manual
287 0 : } else if(poGDS->nDataTypeCode == 3){
288 : float fVal;
289 0 : for (i=0;i<nBlockXSize;i++){
290 0 : fVal = (float) *(pszRecord+i*nDataLength);
291 0 : if (fVal == 0)
292 0 : fVal = -9997;
293 0 : else if(fVal == 1)
294 0 : fVal = -9998;
295 0 : else if(fVal == 255)
296 0 : fVal = -9999;
297 : else
298 0 : fVal = poGDS->fNyquistVelocity * (fVal - 128)/127;
299 0 : ((float *) pImage)[i] = fVal;
300 : }
301 : }
302 :
303 263 : return CE_None;
304 : }
305 :
306 : /************************************************************************/
307 : /* SetNoDataValue() */
308 : /************************************************************************/
309 :
310 2 : CPLErr IRISRasterBand::SetNoDataValue( double dfNoData )
311 :
312 : {
313 2 : IRISDataset *poGDS = (IRISDataset *) poDS;
314 : // if( poGDS->bNoDataSet && poGDS->dfNoDataValue == dfNoData )
315 : // return CE_None;
316 :
317 :
318 2 : poGDS->bNoDataSet = TRUE;
319 2 : poGDS->dfNoDataValue = dfNoData;
320 :
321 2 : return CE_None;
322 : }
323 :
324 : /************************************************************************/
325 : /* GetNoDataValue() */
326 : /************************************************************************/
327 :
328 0 : double IRISRasterBand::GetNoDataValue( int * pbSuccess )
329 :
330 : {
331 0 : IRISDataset *poGDS = (IRISDataset *) poDS;
332 :
333 :
334 0 : if( poGDS->bNoDataSet )
335 : {
336 0 : if( pbSuccess )
337 0 : *pbSuccess = TRUE;
338 :
339 0 : return poGDS->dfNoDataValue;
340 : }
341 :
342 0 : return GDALPamRasterBand::GetNoDataValue( pbSuccess );
343 : }
344 :
345 :
346 : /************************************************************************/
347 : /* ==================================================================== */
348 : /* IRISDataset */
349 : /* ==================================================================== */
350 : /************************************************************************/
351 :
352 : /************************************************************************/
353 : /* IRISDataset() */
354 : /************************************************************************/
355 :
356 2 : IRISDataset::IRISDataset()
357 :
358 : {
359 2 : bHasLoadedProjection = FALSE;
360 2 : fp = NULL;
361 2 : pszSRS_WKT = NULL;
362 2 : adfGeoTransform[0] = 0.0;
363 2 : adfGeoTransform[1] = 1.0;
364 2 : adfGeoTransform[2] = 0.0;
365 2 : adfGeoTransform[3] = 0.0;
366 2 : adfGeoTransform[4] = 0.0;
367 2 : adfGeoTransform[5] = 1.0;
368 2 : }
369 :
370 : /************************************************************************/
371 : /* ~IRISDataset() */
372 : /************************************************************************/
373 :
374 2 : IRISDataset::~IRISDataset()
375 :
376 : {
377 2 : FlushCache();
378 2 : if( fp != NULL )
379 2 : VSIFCloseL( fp );
380 2 : CPLFree( pszSRS_WKT );
381 2 : }
382 :
383 : /************************************************************************/
384 : /* Calculates the projection and Geotransform */
385 : /************************************************************************/
386 1 : void IRISDataset::LoadProjection()
387 : {
388 1 : bHasLoadedProjection = TRUE;
389 1 : float fEquatorialRadius = float( (CPL_LSBUINT32PTR (abyHeader+220+320+12)))/100; //They give it in cm
390 1 : float fInvFlattening = float( (CPL_LSBUINT32PTR (abyHeader+224+320+12)))/1000000; //Point 3.2.27 pag 3-15
391 : float fFlattening;
392 : float fPolarRadius;
393 :
394 1 : if(fEquatorialRadius == 0){ // if Radius is 0, change to 6371000 Point 3.2.27 pag 3-15 (old IRIS verions)
395 0 : fEquatorialRadius = 6371000;
396 0 : fPolarRadius = fEquatorialRadius;
397 0 : fInvFlattening = 0;
398 0 : fFlattening = 0;
399 : } else {
400 1 : if (fInvFlattening == 0){ //When inverse flattening is infinite, they use 0
401 1 : fFlattening = 0;
402 1 : fPolarRadius = fEquatorialRadius;
403 : } else {
404 0 : fFlattening = 1/fInvFlattening;
405 0 : fPolarRadius = fEquatorialRadius * (1-fFlattening);
406 : }
407 : }
408 :
409 1 : float fCenterLon = 360 * float((CPL_LSBUINT32PTR (abyHeader+112+320+12))) / 4294967295LL;
410 1 : float fCenterLat = 360 * float((CPL_LSBUINT32PTR (abyHeader+108+320+12))) / 4294967295LL;
411 :
412 1 : float fProjRefLon = 360 * float((CPL_LSBUINT32PTR (abyHeader+244+320+12))) / 4294967295LL;
413 1 : float fProjRefLat = 360 * float((CPL_LSBUINT32PTR (abyHeader+240+320+12))) / 4294967295LL;
414 :
415 : float fRadarLocX, fRadarLocY, fScaleX, fScaleY;
416 :
417 1 : fRadarLocX = float (CPL_LSBSINT32PTR (abyHeader + 112 + 12 )) / 1000;
418 1 : fRadarLocY = float (CPL_LSBSINT32PTR (abyHeader + 116 + 12 )) / 1000;
419 :
420 1 : fScaleX = float (CPL_LSBSINT32PTR (abyHeader + 88 + 12 )) / 100;
421 1 : fScaleY = float (CPL_LSBSINT32PTR (abyHeader + 92 + 12 )) / 100;
422 :
423 1 : OGRSpatialReference oSRSOut;
424 :
425 : ////MERCATOR PROJECTION
426 1 : if(EQUAL(aszProjections[nProjectionCode],"Mercator")){
427 1 : OGRCoordinateTransformation *poTransform = NULL;
428 1 : OGRSpatialReference oSRSLatLon;
429 :
430 : oSRSOut.SetGeogCS("unnamed ellipse",
431 : "unknown",
432 : "unnamed",
433 : fEquatorialRadius, fInvFlattening,
434 : "Greenwich", 0.0,
435 1 : "degree", 0.0174532925199433);
436 :
437 1 : oSRSOut.SetMercator(fProjRefLat,fProjRefLon,1,0,0);
438 1 : oSRSOut.exportToWkt(&pszSRS_WKT);
439 :
440 : //The center coordinates are given in LatLon on the defined ellipsoid. Necessary to calculate geotransform.
441 :
442 : oSRSLatLon.SetGeogCS("unnamed ellipse",
443 : "unknown",
444 : "unnamed",
445 : fEquatorialRadius, fInvFlattening,
446 : "Greenwich", 0.0,
447 1 : "degree", 0.0174532925199433);
448 :
449 : poTransform = OGRCreateCoordinateTransformation( &oSRSLatLon,
450 1 : &oSRSOut );
451 1 : std::pair <double,double> oPositionX2 = GeodesicCalculation(fCenterLat, fCenterLon, 90, fScaleX, fEquatorialRadius, fPolarRadius, fFlattening);
452 1 : std::pair <double,double> oPositionY2 = GeodesicCalculation(fCenterLat, fCenterLon, 0, fScaleY, fEquatorialRadius, fPolarRadius, fFlattening);
453 :
454 : double dfLon2, dfLat2;
455 1 : dfLon2 = oPositionX2.first;
456 1 : dfLat2 = oPositionY2.second;
457 : double dfX, dfY, dfX2, dfY2;
458 1 : dfX = fCenterLon ;
459 1 : dfY = fCenterLat ;
460 1 : dfX2 = dfLon2;
461 1 : dfY2 = dfLat2;
462 :
463 1 : if( poTransform == NULL || !poTransform->Transform( 1, &dfX, &dfY ) )
464 0 : CPLError( CE_Failure, CPLE_None, "Transformation Failed\n" );
465 :
466 1 : if( poTransform == NULL || !poTransform->Transform( 1, &dfX2, &dfY2 ) )
467 0 : CPLError( CE_Failure, CPLE_None, "Transformation Failed\n" );
468 :
469 1 : adfGeoTransform[0] = dfX - (fRadarLocX * (dfX2 - dfX));
470 1 : adfGeoTransform[1] = dfX2 - dfX;
471 1 : adfGeoTransform[2] = 0.0;
472 1 : adfGeoTransform[3] = dfY + (fRadarLocY * (dfY2 - dfY));
473 1 : adfGeoTransform[4] = 0.0;
474 1 : adfGeoTransform[5] = -1*(dfY2 - dfY);
475 :
476 1 : delete poTransform;
477 :
478 0 : }else if(EQUAL(aszProjections[nProjectionCode],"Azimutal equidistant")){
479 :
480 : oSRSOut.SetGeogCS("unnamed ellipse",
481 : "unknown",
482 : "unnamed",
483 : fEquatorialRadius, fInvFlattening,
484 : "Greenwich", 0.0,
485 0 : "degree", 0.0174532925199433);
486 0 : oSRSOut.SetAE(fProjRefLat,fProjRefLon,0,0);
487 0 : oSRSOut.exportToWkt(&pszSRS_WKT) ;
488 0 : adfGeoTransform[0] = -1*(fRadarLocX*fScaleX);
489 0 : adfGeoTransform[1] = fScaleX;
490 0 : adfGeoTransform[2] = 0.0;
491 0 : adfGeoTransform[3] = fRadarLocY*fScaleY;
492 0 : adfGeoTransform[4] = 0.0;
493 0 : adfGeoTransform[5] = -1*fScaleY;
494 : //When the projection is different from Mercator or Azimutal equidistant, we set a standard geotransform
495 : } else {
496 0 : adfGeoTransform[0] = -1*(fRadarLocX*fScaleX);
497 0 : adfGeoTransform[1] = fScaleX;
498 0 : adfGeoTransform[2] = 0.0;
499 0 : adfGeoTransform[3] = fRadarLocY*fScaleY;
500 0 : adfGeoTransform[4] = 0.0;
501 0 : adfGeoTransform[5] = -1*fScaleY;
502 1 : }
503 :
504 1 : }
505 :
506 : /******************************************************************************/
507 : /* The geotransform in Mercator projection must be calculated transforming */
508 : /* distance to degrees over the ellipsoid, using Vincenty's formula. */
509 : /* The following method is ported from a version for Javascript by Chris */
510 : /* Veness distributed under a CC-BY 3.0 licence, whose conditions is that the */
511 : /* following copyright notice is retained as well as the link to : */
512 : /* http://www.movable-type.co.uk/scripts/latlong-vincenty-direct.html */
513 : /******************************************************************************/
514 :
515 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
516 : /* Vincenty Direct Solution of Geodesics on the Ellipsoid (c) Chris Veness 2005-2012 */
517 : /* */
518 : /* from: Vincenty direct formula - T Vincenty, "Direct and Inverse Solutions of Geodesics on the */
519 : /* Ellipsoid with application of nested equations", Survey Review, vol XXII no 176, 1975 */
520 : /* http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf */
521 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
522 :
523 2 : std::pair <double,double> IRISDataset::GeodesicCalculation(float fLat, float fLon, float fAngle, float fDist, float fEquatorialRadius, float fPolarRadius, float fFlattening)
524 : {
525 2 : std::pair <double,double> oOutput;
526 2 : double dfAlpha1 = DEG2RAD * fAngle;
527 2 : double dfSinAlpha1 = sin(dfAlpha1);
528 2 : double dfCosAlpha1 = cos(dfAlpha1);
529 :
530 2 : double dfTanU1 = (1-fFlattening) * tan(fLat*DEG2RAD);
531 2 : double dfCosU1 = 1 / sqrt((1 + dfTanU1*dfTanU1));
532 2 : double dfSinU1 = dfTanU1*dfCosU1;
533 :
534 2 : double dfSigma1 = atan2(dfTanU1, dfCosAlpha1);
535 2 : double dfSinAlpha = dfCosU1 * dfSinAlpha1;
536 2 : double dfCosSqAlpha = 1 - dfSinAlpha*dfSinAlpha;
537 2 : double dfUSq = dfCosSqAlpha * (fEquatorialRadius*fEquatorialRadius - fPolarRadius*fPolarRadius) / (fPolarRadius*fPolarRadius);
538 2 : double dfA = 1 + dfUSq/16384*(4096+dfUSq*(-768+dfUSq*(320-175*dfUSq)));
539 2 : double dfB = dfUSq/1024 * (256+dfUSq*(-128+dfUSq*(74-47*dfUSq)));
540 :
541 2 : double dfSigma = fDist / (fPolarRadius*dfA);
542 2 : double dfSigmaP = 2*M_PI;
543 :
544 : double dfSinSigma;
545 : double dfCosSigma;
546 : double dfCos2SigmaM;
547 : double dfDeltaSigma;
548 :
549 6 : while (fabs(dfSigma-dfSigmaP) > 1e-12) {
550 2 : dfCos2SigmaM = cos(2*dfSigma1 + dfSigma);
551 2 : dfSinSigma = sin(dfSigma);
552 2 : dfCosSigma = cos(dfSigma);
553 : dfDeltaSigma = dfB*dfSinSigma*(dfCos2SigmaM+dfB/4*(dfCosSigma*(-1+2*dfCos2SigmaM*dfCos2SigmaM)-
554 2 : dfB/6*dfCos2SigmaM*(-3+4*dfSinSigma*dfSinSigma)*(-3+4*dfCos2SigmaM*dfCos2SigmaM)));
555 2 : dfSigmaP = dfSigma;
556 2 : dfSigma = fDist / (fPolarRadius*dfA) + dfDeltaSigma;
557 : }
558 :
559 2 : double dfTmp = dfSinU1*dfSinSigma - dfCosU1*dfCosSigma*dfCosAlpha1;
560 : double dfLat2 = atan2(dfSinU1*dfCosSigma + dfCosU1*dfSinSigma*dfCosAlpha1,
561 2 : (1-fFlattening)*sqrt(dfSinAlpha*dfSinAlpha + dfTmp*dfTmp));
562 2 : double dfLambda = atan2(dfSinSigma*dfSinAlpha1, dfCosU1*dfCosSigma - dfSinU1*dfSinSigma*dfCosAlpha1);
563 2 : double dfC = fFlattening/16*dfCosSqAlpha*(4+fFlattening*(4-3*dfCosSqAlpha));
564 : double dfL = dfLambda - (1-dfC) * fFlattening * dfSinAlpha *
565 2 : (dfSigma + dfC*dfSinSigma*(dfCos2SigmaM+dfC*dfCosSigma*(-1+2*dfCos2SigmaM*dfCos2SigmaM)));
566 2 : double dfLon2 = fLon*DEG2RAD+dfL;
567 2 : if (dfLon2 > M_PI)
568 0 : dfLon2 = dfLon2 - 2*M_PI;
569 2 : if (dfLon2 < -1*M_PI)
570 0 : dfLon2 = dfLon2 + 2*M_PI;
571 2 : oOutput.first = dfLon2*RAD2DEG;
572 2 : oOutput.second = dfLat2*RAD2DEG;
573 :
574 2 : return oOutput;
575 : }
576 :
577 : /************************************************************************/
578 : /* GetGeoTransform() */
579 : /************************************************************************/
580 :
581 1 : CPLErr IRISDataset::GetGeoTransform( double * padfTransform )
582 :
583 : {
584 1 : if (!bHasLoadedProjection)
585 0 : LoadProjection();
586 1 : memcpy( padfTransform, adfGeoTransform, sizeof(double)*6 );
587 1 : return CE_None;
588 : }
589 :
590 : /************************************************************************/
591 : /* GetProjectionRef() */
592 : /************************************************************************/
593 :
594 1 : const char *IRISDataset::GetProjectionRef(){
595 1 : if (!bHasLoadedProjection)
596 1 : LoadProjection();
597 1 : return pszSRS_WKT;
598 : }
599 :
600 : /************************************************************************/
601 : /* Identify() */
602 : /************************************************************************/
603 :
604 11416 : int IRISDataset::Identify( GDALOpenInfo * poOpenInfo )
605 :
606 : {
607 :
608 : /* -------------------------------------------------------------------- */
609 : /* Confirm that the file is an IRIS file */
610 : /* -------------------------------------------------------------------- */
611 : //Si no el posem, peta al fer el translate, quan s'obre Identify des de GDALIdentifyDriver
612 11416 : if( poOpenInfo->nHeaderBytes < 640 )
613 11296 : return FALSE;
614 :
615 :
616 120 : short nId1 = CPL_LSBSINT16PTR(poOpenInfo->pabyHeader);
617 120 : short nId2 = CPL_LSBSINT16PTR(poOpenInfo->pabyHeader+12);
618 120 : unsigned short nType = CPL_LSBUINT16PTR (poOpenInfo->pabyHeader+24);
619 :
620 : /*Check if the two headers are 27 (product hdr) & 26 (product configuration), and the product type is in the range 1 -> 34*/
621 120 : if( !(nId1 == 27 && nId2 == 26 && nType > 0 && nType < 35) )
622 118 : return FALSE;
623 :
624 2 : return TRUE;
625 : }
626 :
627 : /************************************************************************/
628 : /* FillString() */
629 : /************************************************************************/
630 :
631 14 : static void FillString(char* szBuffer, size_t nBufferSize, void* pSrcBuffer)
632 : {
633 190 : for(size_t i = 0; i < nBufferSize - 1; i++)
634 176 : szBuffer[i] = ((char*)pSrcBuffer)[i];
635 14 : szBuffer[nBufferSize-1] = '\0';
636 14 : }
637 :
638 : /************************************************************************/
639 : /* Open() */
640 : /************************************************************************/
641 :
642 1314 : GDALDataset *IRISDataset::Open( GDALOpenInfo * poOpenInfo )
643 :
644 : {
645 1314 : if (!Identify(poOpenInfo))
646 1312 : return NULL;
647 : /* -------------------------------------------------------------------- */
648 : /* Confirm the requested access is supported. */
649 : /* -------------------------------------------------------------------- */
650 2 : if( poOpenInfo->eAccess == GA_Update )
651 : {
652 : CPLError( CE_Failure, CPLE_NotSupported,
653 : "The IRIS driver does not support update access to existing"
654 0 : " datasets.\n" );
655 0 : return NULL;
656 : }
657 :
658 : /* -------------------------------------------------------------------- */
659 : /* Create a corresponding GDALDataset. */
660 : /* -------------------------------------------------------------------- */
661 : IRISDataset *poDS;
662 :
663 2 : poDS = new IRISDataset();
664 :
665 2 : poDS->fp = VSIFOpenL( poOpenInfo->pszFilename, "rb" );
666 2 : if (poDS->fp == NULL)
667 : {
668 0 : delete poDS;
669 0 : return NULL;
670 : }
671 :
672 : /* -------------------------------------------------------------------- */
673 : /* Read the header. */
674 : /* -------------------------------------------------------------------- */
675 2 : VSIFReadL( poDS->abyHeader, 1, 640, poDS->fp );
676 2 : int nXSize = CPL_LSBSINT32PTR(poDS->abyHeader+100+12);
677 2 : int nYSize = CPL_LSBSINT32PTR(poDS->abyHeader+104+12);
678 2 : int nNumBands = CPL_LSBSINT32PTR(poDS->abyHeader+108+12);
679 :
680 2 : poDS->nRasterXSize = nXSize;
681 :
682 2 : poDS->nRasterYSize = nYSize;
683 2 : if (poDS->nRasterXSize <= 0 || poDS->nRasterYSize <= 0 )
684 : {
685 : CPLError( CE_Failure, CPLE_AppDefined,
686 : "Invalid dimensions : %d x %d",
687 0 : poDS->nRasterXSize, poDS->nRasterYSize);
688 0 : delete poDS;
689 0 : return NULL;
690 : }
691 :
692 2 : if( !GDALCheckBandCount(nNumBands, TRUE) )
693 : {
694 0 : delete poDS;
695 0 : return NULL;
696 : }
697 :
698 : /* -------------------------------------------------------------------- */
699 : /* Create band information objects. */
700 : /* -------------------------------------------------------------------- */
701 8 : for (int iBandNum = 1; iBandNum <= nNumBands; iBandNum++) {
702 2 : poDS->SetBand( iBandNum, new IRISRasterBand( poDS, iBandNum ));
703 2 : poDS->GetRasterBand(iBandNum)->SetNoDataValue(-9999);
704 : }
705 :
706 : /* -------------------------------------------------------------------- */
707 : /* Setting the Metadata */
708 : /* -------------------------------------------------------------------- */
709 : //See point 3.2.26 at page 3.12 of the manual
710 2 : poDS->nProductCode = CPL_LSBUINT16PTR (poDS->abyHeader+12+12);
711 2 : poDS->SetMetadataItem( "PRODUCT_ID", CPLString().Printf("%d", poDS->nProductCode ));
712 2 : if( poDS->nProductCode >= ARRAY_ELEMENT_COUNT(poDS->aszProductNames) )
713 : {
714 0 : delete poDS;
715 0 : return NULL;
716 : }
717 :
718 2 : poDS->SetMetadataItem( "PRODUCT",poDS->aszProductNames[poDS->nProductCode]);
719 :
720 2 : poDS->nDataTypeCode = CPL_LSBUINT16PTR (poDS->abyHeader+130+12);
721 2 : if( poDS->nDataTypeCode >= ARRAY_ELEMENT_COUNT(poDS->aszDataTypeCodes) )
722 : {
723 0 : delete poDS;
724 0 : return NULL;
725 : }
726 2 : poDS->SetMetadataItem( "DATA_TYPE_CODE",poDS->aszDataTypeCodes[poDS->nDataTypeCode]);
727 :
728 2 : if( poDS->nDataTypeCode >= ARRAY_ELEMENT_COUNT(poDS->aszDataTypes) )
729 : {
730 0 : delete poDS;
731 0 : return NULL;
732 : }
733 2 : poDS->SetMetadataItem( "DATA_TYPE",poDS->aszDataTypes[poDS->nDataTypeCode]);
734 :
735 2 : unsigned short nDataTypeInputCode = CPL_LSBUINT16PTR (poDS->abyHeader+144+12);
736 2 : if( nDataTypeInputCode >= ARRAY_ELEMENT_COUNT(poDS->aszDataTypeCodes) )
737 : {
738 0 : delete poDS;
739 0 : return NULL;
740 : }
741 2 : poDS->SetMetadataItem( "DATA_TYPE_INPUT_CODE",poDS->aszDataTypeCodes[nDataTypeInputCode]);
742 :
743 2 : unsigned short nDataTypeInput = CPL_LSBUINT16PTR (poDS->abyHeader+144+12);
744 2 : if( nDataTypeInput >= ARRAY_ELEMENT_COUNT(poDS->aszDataTypes) )
745 : {
746 0 : delete poDS;
747 0 : return NULL;
748 : }
749 2 : poDS->SetMetadataItem( "DATA_TYPE_INPUT",poDS->aszDataTypes[nDataTypeInput]);
750 :
751 2 : poDS->nProjectionCode = * (unsigned char *) (poDS->abyHeader+146+12);
752 2 : if( poDS->nProjectionCode >= ARRAY_ELEMENT_COUNT(poDS->aszProjections) )
753 : {
754 0 : delete poDS;
755 0 : return NULL;
756 : }
757 :
758 : ////TIMES
759 2 : int nSeconds = CPL_LSBSINT32PTR(poDS->abyHeader+20+12);
760 :
761 2 : int nHour = (nSeconds - (nSeconds%3600)) /3600;
762 2 : int nMinute = ((nSeconds - nHour * 3600) - (nSeconds - nHour * 3600)%60)/ 60;
763 2 : int nSecond = nSeconds - nHour * 3600 - nMinute * 60;
764 :
765 2 : short nYear = CPL_LSBSINT16PTR(poDS->abyHeader+26+12);
766 2 : short nMonth = CPL_LSBSINT16PTR(poDS->abyHeader+28+12);
767 2 : short nDay = CPL_LSBSINT16PTR(poDS->abyHeader+30+12);
768 :
769 2 : poDS->SetMetadataItem( "TIME_PRODUCT_GENERATED", CPLString().Printf("%d-%02d-%02d %02d:%02d:%02d", nYear, nMonth, nDay, nHour, nMinute, nSecond ) );
770 :
771 :
772 2 : nSeconds = CPL_LSBSINT32PTR(poDS->abyHeader+32+12);
773 :
774 2 : nHour = (nSeconds - (nSeconds%3600)) /3600;
775 2 : nMinute = ((nSeconds - nHour * 3600) - (nSeconds - nHour * 3600)%60)/ 60;
776 2 : nSecond = nSeconds - nHour * 3600 - nMinute * 60;
777 :
778 2 : nYear = CPL_LSBSINT16PTR(poDS->abyHeader+26+12);
779 2 : nMonth = CPL_LSBSINT16PTR(poDS->abyHeader+28+12);
780 2 : nDay = CPL_LSBSINT16PTR(poDS->abyHeader+30+12);
781 :
782 4 : poDS->SetMetadataItem( "TIME_INPUT_INGEST_SWEEP", CPLString().Printf("%d-%02d-%02d %02d:%02d:%02d", nYear, nMonth, nDay, nHour, nMinute, nSecond ) );
783 :
784 : ///Site and task information
785 :
786 2 : char szSiteName[17] = ""; //Must have one extra char for string end!
787 2 : char szVersionName[9] = "";
788 :
789 2 : FillString(szSiteName, sizeof(szSiteName), poDS->abyHeader+320+12);
790 2 : FillString(szVersionName, sizeof(szVersionName), poDS->abyHeader+16+320+12);
791 2 : poDS->SetMetadataItem( "PRODUCT_SITE_NAME",szSiteName);
792 2 : poDS->SetMetadataItem( "PRODUCT_SITE_IRIS_VERSION",szVersionName);
793 :
794 2 : FillString(szSiteName, sizeof(szSiteName), poDS->abyHeader+90+320+12);
795 2 : FillString(szVersionName, sizeof(szVersionName), poDS->abyHeader+24+320+12);
796 2 : poDS->SetMetadataItem( "INGEST_SITE_NAME",szSiteName);
797 2 : poDS->SetMetadataItem( "INGEST_SITE_IRIS_VERSION",szVersionName);
798 :
799 2 : FillString(szSiteName, sizeof(szSiteName), poDS->abyHeader+74+320+12);
800 2 : poDS->SetMetadataItem( "INGEST_HARDWARE_NAME",szSiteName);
801 :
802 2 : char szConfigFile[13] = "";
803 2 : FillString(szConfigFile, sizeof(szConfigFile), poDS->abyHeader+62+12);
804 2 : poDS->SetMetadataItem( "PRODUCT_CONFIGURATION_NAME",szConfigFile);
805 :
806 2 : char szTaskName[13] = "";
807 2 : FillString(szTaskName, sizeof(szTaskName), poDS->abyHeader+74+12);
808 2 : poDS->SetMetadataItem( "TASK_NAME",szTaskName);
809 :
810 2 : short nRadarHeight = CPL_LSBSINT16PTR(poDS->abyHeader+284+320+12);
811 4 : poDS->SetMetadataItem( "RADAR_HEIGHT",CPLString().Printf("%d m",nRadarHeight));
812 2 : short nGroundHeight = CPL_LSBSINT16PTR(poDS->abyHeader+118+320+12);
813 4 : poDS->SetMetadataItem( "GROUND_HEIGHT",CPLString().Printf("%d m",nRadarHeight-nGroundHeight)); //Ground height over the sea level
814 :
815 2 : unsigned short nFlags = CPL_LSBUINT16PTR (poDS->abyHeader+86+12);
816 : //Get eleventh bit
817 2 : nFlags=nFlags<<4;
818 2 : nFlags=nFlags>>15;
819 2 : if (nFlags == 1){
820 1 : poDS->SetMetadataItem( "COMPOSITED_PRODUCT","YES");
821 1 : unsigned int compositedMask = CPL_LSBUINT32PTR (poDS->abyHeader+232+320+12);
822 1 : poDS->SetMetadataItem( "COMPOSITED_PRODUCT_MASK",CPLString().Printf("0x%08x",compositedMask));
823 : } else{
824 1 : poDS->SetMetadataItem( "COMPOSITED_PRODUCT","NO");
825 : }
826 :
827 : //Wave values
828 2 : poDS->SetMetadataItem( "PRF",CPLString().Printf("%d Hz",CPL_LSBSINT32PTR(poDS->abyHeader+120+320+12)));
829 4 : poDS->SetMetadataItem( "WAVELENGTH",CPLString().Printf("%4.2f cm",(float) CPL_LSBSINT32PTR(poDS->abyHeader+148+320+12)/100));
830 2 : unsigned short nPolarizationType = CPL_LSBUINT16PTR (poDS->abyHeader+172+320+12);
831 2 : float fNyquist = (CPL_LSBSINT32PTR(poDS->abyHeader+120+320+12))*((float) CPL_LSBSINT32PTR(poDS->abyHeader+148+320+12)/10000)/4; //See section 3.3.37 & 3.2.54
832 2 : if (nPolarizationType == 1)
833 0 : fNyquist = fNyquist * 2;
834 2 : else if(nPolarizationType == 2)
835 0 : fNyquist = fNyquist * 3;
836 2 : else if(nPolarizationType == 3)
837 0 : fNyquist = fNyquist * 4;
838 2 : poDS->fNyquistVelocity = fNyquist;
839 2 : poDS->SetMetadataItem( "NYQUIST_VELOCITY",CPLString().Printf("%.2f m/s",fNyquist));
840 :
841 : ///Product dependent metadata (stored in 80 bytes fromm 162 bytes at the product header) See point 3.2.30 at page 3.19 of the manual
842 : //See point 3.2.25 at page 3.12 of the manual
843 2 : if (EQUAL(poDS->aszProductNames[poDS->nProductCode],"PPI")){
844 : //Degrees = 360 * (Binary Angle)*2^N
845 : //float fElevation = 360 * float((CPL_LSBUINT16PTR (poDS->abyHeader+164+12))) / 65536;
846 1 : float fElevation = 360 * float((CPL_LSBSINT16PTR (poDS->abyHeader+164+12))) / 65536;
847 :
848 1 : poDS->SetMetadataItem( "PPI_ELEVATION_ANGLE",CPLString().Printf("%f",fElevation));
849 1 : if (EQUAL(poDS->aszDataTypeCodes[poDS->nDataTypeCode],"dBZ"))
850 1 : poDS->SetMetadataItem( "DATA_TYPE_UNITS","dBZ");
851 : else
852 0 : poDS->SetMetadataItem( "DATA_TYPE_UNITS","m/s");
853 : //See point 3.2.2 at page 3.2 of the manual
854 1 : } else if (EQUAL(poDS->aszProductNames[poDS->nProductCode],"CAPPI")){
855 1 : float fElevation = ((float) CPL_LSBSINT32PTR(poDS->abyHeader+4+164+12))/100;
856 1 : poDS->SetMetadataItem( "CAPPI_HEIGHT",CPLString().Printf("%.1f m",fElevation));
857 1 : float fAzimuthSmoothingForShear = 360 * float((CPL_LSBUINT16PTR (poDS->abyHeader+10+164+12))) / 65536;
858 2 : poDS->SetMetadataItem( "AZIMUTH_SMOOTHING_FOR_SHEAR" ,CPLString().Printf("%.1f", fAzimuthSmoothingForShear));
859 1 : unsigned int nMaxAgeVVPCorrection = CPL_LSBUINT32PTR (poDS->abyHeader+24+164+12);
860 2 : poDS->SetMetadataItem( "MAX_AGE_FOR_SHEAR_VVP_CORRECTION" ,CPLString().Printf("%d s", nMaxAgeVVPCorrection));
861 1 : if (EQUAL(poDS->aszDataTypeCodes[poDS->nDataTypeCode],"dBZ"))
862 1 : poDS->SetMetadataItem( "DATA_TYPE_UNITS","dBZ");
863 : else
864 0 : poDS->SetMetadataItem( "DATA_TYPE_UNITS","m/s");
865 : //See point 3.2.32 at page 3.19 of the manual
866 0 : } else if (EQUAL(poDS->aszProductNames[poDS->nProductCode],"RAIN1") || EQUAL(poDS->aszProductNames[poDS->nProductCode],"RAINN")){
867 0 : short nNumProducts = CPL_LSBSINT16PTR(poDS->abyHeader+170+320+12);
868 0 : poDS->SetMetadataItem( "NUM_FILES_USED",CPLString().Printf("%d",nNumProducts));
869 :
870 0 : float fMinZAcum= (float)((CPL_LSBUINT32PTR (poDS->abyHeader+164+12))-32768)/1000;
871 0 : poDS->SetMetadataItem( "MINIMUM_Z_TO_ACUMULATE",CPLString().Printf("%f",fMinZAcum));
872 :
873 0 : unsigned short nSecondsOfAccumulation = CPL_LSBUINT16PTR (poDS->abyHeader+6+164+12);
874 0 : poDS->SetMetadataItem( "SECONDS_OF_ACCUMULATION",CPLString().Printf("%d s",nSecondsOfAccumulation));
875 :
876 0 : unsigned int nSpanInputFiles = CPL_LSBUINT32PTR (poDS->abyHeader+24+164+12);
877 0 : poDS->SetMetadataItem( "SPAN_OF_INPUT_FILES",CPLString().Printf("%d s",nSpanInputFiles));
878 0 : poDS->SetMetadataItem( "DATA_TYPE_UNITS","mm");
879 :
880 0 : char szInputProductName[13] = "";
881 0 : for(int k=0; k<12;k++)
882 0 : szInputProductName[k] = * (char *) (poDS->abyHeader+k+12+164+12);
883 0 : poDS->SetMetadataItem( "INPUT_PRODUCT_NAME",CPLString().Printf("%s",szInputProductName));
884 :
885 0 : if (EQUAL(poDS->aszProductNames[poDS->nProductCode],"RAINN"))
886 0 : poDS->SetMetadataItem( "NUM_HOURS_ACCUMULATE",CPLString().Printf("%d",CPL_LSBUINT16PTR (poDS->abyHeader+10+164+12)));
887 :
888 : //See point 3.2.73 at page 3.36 of the manual
889 0 : } else if (EQUAL(poDS->aszProductNames[poDS->nProductCode],"VIL")){
890 0 : float fBottomHeigthInterval = (float) CPL_LSBSINT32PTR(poDS->abyHeader+4+164+12) / 100;
891 0 : poDS->SetMetadataItem( "BOTTOM_OF_HEIGTH_INTERVAL",CPLString().Printf("%.1f m",fBottomHeigthInterval));
892 0 : float fTopHeigthInterval = (float) CPL_LSBSINT32PTR(poDS->abyHeader+8+164+12) / 100;
893 0 : poDS->SetMetadataItem( "TOP_OF_HEIGTH_INTERVAL",CPLString().Printf("%.1f m",fTopHeigthInterval));
894 0 : poDS->SetMetadataItem( "VIL_DENSITY_NOT_AVAILABLE_VALUE","-1");
895 0 : poDS->SetMetadataItem( "DATA_TYPE_UNITS","mm");
896 : //See point 3.2.68 at page 3.36 of the manual
897 0 : } else if (EQUAL(poDS->aszProductNames[poDS->nProductCode],"TOPS")){
898 0 : float fZThreshold = (float) CPL_LSBSINT16PTR(poDS->abyHeader+4+164+12) / 16;
899 0 : poDS->SetMetadataItem( "Z_THRESHOLD",CPLString().Printf("%.1f dBZ",fZThreshold));
900 0 : poDS->SetMetadataItem( "ECHO_TOPS_NOT_AVAILABLE_VALUE","-1");
901 0 : poDS->SetMetadataItem( "DATA_TYPE_UNITS","km");
902 : //See point 3.2.20 at page 3.10 of the manual
903 0 : } else if (EQUAL(poDS->aszProductNames[poDS->nProductCode],"MAX")){
904 0 : float fBottomInterval = (float) CPL_LSBSINT32PTR(poDS->abyHeader+4+164+12) / 100;
905 0 : poDS->SetMetadataItem( "BOTTOM_OF_INTERVAL",CPLString().Printf("%.1f m",fBottomInterval));
906 0 : float fTopInterval = (float) CPL_LSBSINT32PTR(poDS->abyHeader+8+164+12) / 100;
907 0 : poDS->SetMetadataItem( "TOP_OF_INTERVAL",CPLString().Printf("%.1f m",fTopInterval));
908 0 : int nNumPixelsSidePanels = CPL_LSBSINT32PTR(poDS->abyHeader+12+164+12);
909 0 : poDS->SetMetadataItem( "NUM_PIXELS_SIDE_PANELS",CPLString().Printf("%d",nNumPixelsSidePanels));
910 0 : short nHorizontalSmootherSidePanels = CPL_LSBSINT16PTR(poDS->abyHeader+16+164+12);
911 0 : poDS->SetMetadataItem( "HORIZONTAL_SMOOTHER_SIDE_PANELS",CPLString().Printf("%d",nHorizontalSmootherSidePanels));
912 0 : short nVerticalSmootherSidePanels = CPL_LSBSINT16PTR(poDS->abyHeader+18+164+12);
913 0 : poDS->SetMetadataItem( "VERTICAL_SMOOTHER_SIDE_PANELS",CPLString().Printf("%d",nVerticalSmootherSidePanels));
914 : }
915 :
916 : /* -------------------------------------------------------------------- */
917 : /* Initialize any PAM information. */
918 : /* -------------------------------------------------------------------- */
919 2 : poDS->SetDescription( poOpenInfo->pszFilename );
920 2 : poDS->TryLoadXML();
921 :
922 : /* -------------------------------------------------------------------- */
923 : /* Check for overviews. */
924 : /* -------------------------------------------------------------------- */
925 2 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
926 :
927 2 : return( poDS );
928 : }
929 :
930 : /************************************************************************/
931 : /* GDALRegister_IRIS() */
932 : /************************************************************************/
933 :
934 610 : void GDALRegister_IRIS()
935 :
936 : {
937 : GDALDriver *poDriver;
938 :
939 610 : if( GDALGetDriverByName( "IRIS" ) == NULL )
940 : {
941 588 : poDriver = new GDALDriver();
942 :
943 588 : poDriver->SetDescription( "IRIS" );
944 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
945 588 : "IRIS data (.PPI, .CAPPi etc)" );
946 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
947 588 : "frmt_various.html#IRIS" );
948 588 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "ppi" );
949 588 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
950 :
951 588 : poDriver->pfnOpen = IRISDataset::Open;
952 588 : poDriver->pfnIdentify = IRISDataset::Identify;
953 :
954 588 : GetGDALDriverManager()->RegisterDriver( poDriver );
955 : }
956 610 : }
|
