1 : /******************************************************************************
2 : * $Id: irisdataset.cpp 25302 2012-12-13 18:48:59Z 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 25302 2012-12-13 18:48:59Z 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 263 : VSIFSeekL( poGDS->fp, 640 + nDataLength*poGDS->GetRasterXSize()*poGDS->GetRasterYSize()*(this->nBand-1) + nBlockXSize*nDataLength*(poGDS->GetRasterYSize()-1-nBlockYOff), SEEK_SET );
227 :
228 263 : VSIFReadL( pszRecord, 1, nBlockXSize*nDataLength, poGDS->fp );
229 :
230 :
231 : //If datatype is dbZ:
232 : //See point 3.3.5 at page 3.42 of the manual
233 263 : if(poGDS->nDataTypeCode == 2){
234 : float fVal;
235 68384 : for (i=0;i<nBlockXSize;i++){
236 68121 : fVal = (((float) *(pszRecord+i*nDataLength)) -64)/2.0;
237 68121 : if (fVal == 95.5)
238 0 : fVal = -9999;
239 68121 : ((float *) pImage)[i] = fVal;
240 : }
241 : //Fliquid2 (Rain1 & Rainn products)
242 : //See point 3.3.11 at page 3.43 of the manual
243 0 : } else if(poGDS->nDataTypeCode == 37){
244 : unsigned short nVal, nExp, nMantissa;
245 0 : float fVal2=0;
246 0 : for (i=0;i<nBlockXSize;i++){
247 0 : nVal = CPL_LSBUINT16PTR(pszRecord+i*nDataLength);
248 0 : nExp = nVal>>12;
249 0 : nMantissa = nVal - (nExp<<12);
250 0 : if (nVal == 65535)
251 0 : fVal2 = -9999;
252 0 : else if (nExp == 0)
253 0 : fVal2 = (float) nMantissa / 1000.0;
254 : else
255 0 : fVal2 = (float)((nMantissa+4096)<<(nExp-1))/1000.0;
256 0 : ((float *) pImage)[i] = fVal2;
257 : }
258 : //VIL2 (VIL products)
259 : //See point 3.3.41 at page 3.54 of the manual
260 0 : } else if(poGDS->nDataTypeCode == 33){
261 : float fVal;
262 0 : for (i=0;i<nBlockXSize;i++){
263 0 : fVal = (float) CPL_LSBUINT16PTR(pszRecord+i*nDataLength);
264 0 : if (fVal == 65535)
265 0 : ((float *) pImage)[i] = -9999;
266 0 : else if (fVal == 0)
267 0 : ((float *) pImage)[i] = -1;
268 : else
269 0 : ((float *) pImage)[i] = (fVal-1)/1000;
270 : }
271 : //HEIGTH (TOPS products)
272 : //See point 3.3.14 at page 3.46 of the manual
273 0 : } else if(poGDS->nDataTypeCode == 32){
274 : unsigned char nVal;
275 0 : for (i=0;i<nBlockXSize;i++){
276 0 : nVal = *(pszRecord+i*nDataLength) ;
277 0 : if (nVal == 255)
278 0 : ((float *) pImage)[i] = -9999;
279 0 : else if (nVal == 0)
280 0 : ((float *) pImage)[i] = -1;
281 : else
282 0 : ((float *) pImage)[i] = ((float) nVal - 1) / 10;
283 : }
284 : //VEL (Velocity 1-Byte in PPI & others)
285 : //See point 3.3.37 at page 3.53 of the manual
286 0 : } else if(poGDS->nDataTypeCode == 3){
287 : float fVal;
288 0 : for (i=0;i<nBlockXSize;i++){
289 0 : fVal = (float) *(pszRecord+i*nDataLength);
290 0 : if (fVal == 0)
291 0 : fVal = -9997;
292 0 : else if(fVal == 1)
293 0 : fVal = -9998;
294 0 : else if(fVal == 255)
295 0 : fVal = -9999;
296 : else
297 0 : fVal = poGDS->fNyquistVelocity * (fVal - 128)/127;
298 0 : ((float *) pImage)[i] = fVal;
299 : }
300 : }
301 :
302 263 : return CE_None;
303 : }
304 :
305 : /************************************************************************/
306 : /* SetNoDataValue() */
307 : /************************************************************************/
308 :
309 2 : CPLErr IRISRasterBand::SetNoDataValue( double dfNoData )
310 :
311 : {
312 2 : IRISDataset *poGDS = (IRISDataset *) poDS;
313 : // if( poGDS->bNoDataSet && poGDS->dfNoDataValue == dfNoData )
314 : // return CE_None;
315 :
316 :
317 2 : poGDS->bNoDataSet = TRUE;
318 2 : poGDS->dfNoDataValue = dfNoData;
319 :
320 2 : return CE_None;
321 : }
322 :
323 : /************************************************************************/
324 : /* GetNoDataValue() */
325 : /************************************************************************/
326 :
327 0 : double IRISRasterBand::GetNoDataValue( int * pbSuccess )
328 :
329 : {
330 0 : IRISDataset *poGDS = (IRISDataset *) poDS;
331 :
332 :
333 0 : if( poGDS->bNoDataSet )
334 : {
335 0 : if( pbSuccess )
336 0 : *pbSuccess = TRUE;
337 :
338 0 : return poGDS->dfNoDataValue;
339 : }
340 :
341 0 : return GDALPamRasterBand::GetNoDataValue( pbSuccess );
342 : }
343 :
344 :
345 : /************************************************************************/
346 : /* ==================================================================== */
347 : /* IRISDataset */
348 : /* ==================================================================== */
349 : /************************************************************************/
350 :
351 : /************************************************************************/
352 : /* IRISDataset() */
353 : /************************************************************************/
354 :
355 2 : IRISDataset::IRISDataset()
356 :
357 : {
358 2 : bHasLoadedProjection = FALSE;
359 2 : fp = NULL;
360 2 : pszSRS_WKT = NULL;
361 2 : adfGeoTransform[0] = 0.0;
362 2 : adfGeoTransform[1] = 1.0;
363 2 : adfGeoTransform[2] = 0.0;
364 2 : adfGeoTransform[3] = 0.0;
365 2 : adfGeoTransform[4] = 0.0;
366 2 : adfGeoTransform[5] = 1.0;
367 2 : }
368 :
369 : /************************************************************************/
370 : /* ~IRISDataset() */
371 : /************************************************************************/
372 :
373 2 : IRISDataset::~IRISDataset()
374 :
375 : {
376 2 : FlushCache();
377 2 : if( fp != NULL )
378 2 : VSIFCloseL( fp );
379 2 : CPLFree( pszSRS_WKT );
380 2 : }
381 :
382 : /************************************************************************/
383 : /* Calculates the projection and Geotransform */
384 : /************************************************************************/
385 1 : void IRISDataset::LoadProjection()
386 : {
387 1 : bHasLoadedProjection = TRUE;
388 1 : float fEquatorialRadius = float( (CPL_LSBUINT32PTR (abyHeader+220+320+12)))/100; //They give it in cm
389 1 : float fInvFlattening = float( (CPL_LSBUINT32PTR (abyHeader+224+320+12)))/1000000; //Point 3.2.27 pag 3-15
390 : float fFlattening;
391 : float fPolarRadius;
392 :
393 1 : if(fEquatorialRadius == 0){ // if Radius is 0, change to 6371000 Point 3.2.27 pag 3-15 (old IRIS verions)
394 0 : fEquatorialRadius = 6371000;
395 0 : fPolarRadius = fEquatorialRadius;
396 0 : fInvFlattening = 0;
397 0 : fFlattening = 0;
398 : } else {
399 1 : if (fInvFlattening == 0){ //When inverse flattening is infinite, they use 0
400 1 : fFlattening = 0;
401 1 : fPolarRadius = fEquatorialRadius;
402 : } else {
403 0 : fFlattening = 1/fInvFlattening;
404 0 : fPolarRadius = fEquatorialRadius * (1-fFlattening);
405 : }
406 : }
407 :
408 1 : float fCenterLon = 360 * float((CPL_LSBUINT32PTR (abyHeader+112+320+12))) / 4294967295LL;
409 1 : float fCenterLat = 360 * float((CPL_LSBUINT32PTR (abyHeader+108+320+12))) / 4294967295LL;
410 :
411 1 : float fProjRefLon = 360 * float((CPL_LSBUINT32PTR (abyHeader+244+320+12))) / 4294967295LL;
412 1 : float fProjRefLat = 360 * float((CPL_LSBUINT32PTR (abyHeader+240+320+12))) / 4294967295LL;
413 :
414 : float fRadarLocX, fRadarLocY, fScaleX, fScaleY;
415 :
416 1 : fRadarLocX = float (CPL_LSBSINT32PTR (abyHeader + 112 + 12 )) / 1000;
417 1 : fRadarLocY = float (CPL_LSBSINT32PTR (abyHeader + 116 + 12 )) / 1000;
418 :
419 1 : fScaleX = float (CPL_LSBSINT32PTR (abyHeader + 88 + 12 )) / 100;
420 1 : fScaleY = float (CPL_LSBSINT32PTR (abyHeader + 92 + 12 )) / 100;
421 :
422 1 : OGRSpatialReference oSRSOut;
423 :
424 : ////MERCATOR PROJECTION
425 1 : if(EQUAL(aszProjections[nProjectionCode],"Mercator")){
426 1 : OGRCoordinateTransformation *poTransform = NULL;
427 1 : OGRSpatialReference oSRSLatLon;
428 :
429 : oSRSOut.SetGeogCS("unnamed ellipse",
430 : "unknown",
431 : "unnamed",
432 : fEquatorialRadius, fInvFlattening,
433 : "Greenwich", 0.0,
434 1 : "degree", 0.0174532925199433);
435 :
436 1 : oSRSOut.SetMercator(fProjRefLat,fProjRefLon,1,0,0);
437 1 : oSRSOut.exportToWkt(&pszSRS_WKT);
438 :
439 : //The center coordinates are given in LatLon on the defined ellipsoid. Necessary to calculate geotransform.
440 :
441 : oSRSLatLon.SetGeogCS("unnamed ellipse",
442 : "unknown",
443 : "unnamed",
444 : fEquatorialRadius, fInvFlattening,
445 : "Greenwich", 0.0,
446 1 : "degree", 0.0174532925199433);
447 :
448 : poTransform = OGRCreateCoordinateTransformation( &oSRSLatLon,
449 1 : &oSRSOut );
450 1 : std::pair <double,double> oPositionX2 = GeodesicCalculation(fCenterLat, fCenterLon, 90, fScaleX, fEquatorialRadius, fPolarRadius, fFlattening);
451 1 : std::pair <double,double> oPositionY2 = GeodesicCalculation(fCenterLat, fCenterLon, 0, fScaleY, fEquatorialRadius, fPolarRadius, fFlattening);
452 :
453 : double dfLon2, dfLat2;
454 1 : dfLon2 = oPositionX2.first;
455 1 : dfLat2 = oPositionY2.second;
456 : double dfX, dfY, dfX2, dfY2;
457 1 : dfX = fCenterLon ;
458 1 : dfY = fCenterLat ;
459 1 : dfX2 = dfLon2;
460 1 : dfY2 = dfLat2;
461 :
462 1 : if( poTransform == NULL || !poTransform->Transform( 1, &dfX, &dfY ) )
463 0 : CPLError( CE_Failure, CPLE_None, "Transformation Failed\n" );
464 :
465 1 : if( poTransform == NULL || !poTransform->Transform( 1, &dfX2, &dfY2 ) )
466 0 : CPLError( CE_Failure, CPLE_None, "Transformation Failed\n" );
467 :
468 1 : adfGeoTransform[0] = dfX - (fRadarLocX * (dfX2 - dfX));
469 1 : adfGeoTransform[1] = dfX2 - dfX;
470 1 : adfGeoTransform[2] = 0.0;
471 1 : adfGeoTransform[3] = dfY + (fRadarLocY * (dfY2 - dfY));
472 1 : adfGeoTransform[4] = 0.0;
473 1 : adfGeoTransform[5] = -1*(dfY2 - dfY);
474 :
475 1 : delete poTransform;
476 :
477 0 : }else if(EQUAL(aszProjections[nProjectionCode],"Azimutal equidistant")){
478 :
479 : oSRSOut.SetGeogCS("unnamed ellipse",
480 : "unknown",
481 : "unnamed",
482 : fEquatorialRadius, fInvFlattening,
483 : "Greenwich", 0.0,
484 0 : "degree", 0.0174532925199433);
485 0 : oSRSOut.SetAE(fProjRefLat,fProjRefLon,0,0);
486 0 : oSRSOut.exportToWkt(&pszSRS_WKT) ;
487 0 : adfGeoTransform[0] = -1*(fRadarLocX*fScaleX);
488 0 : adfGeoTransform[1] = fScaleX;
489 0 : adfGeoTransform[2] = 0.0;
490 0 : adfGeoTransform[3] = fRadarLocY*fScaleY;
491 0 : adfGeoTransform[4] = 0.0;
492 0 : adfGeoTransform[5] = -1*fScaleY;
493 : //When the projection is different from Mercator or Azimutal equidistant, we set a standard geotransform
494 : } else {
495 0 : adfGeoTransform[0] = -1*(fRadarLocX*fScaleX);
496 0 : adfGeoTransform[1] = fScaleX;
497 0 : adfGeoTransform[2] = 0.0;
498 0 : adfGeoTransform[3] = fRadarLocY*fScaleY;
499 0 : adfGeoTransform[4] = 0.0;
500 0 : adfGeoTransform[5] = -1*fScaleY;
501 1 : }
502 :
503 1 : }
504 :
505 : /******************************************************************************/
506 : /* The geotransform in Mercator projection must be calculated transforming */
507 : /* distance to degrees over the ellipsoid, using Vincenty's formula. */
508 : /* The following method is ported from a version for Javascript by Chris */
509 : /* Veness distributed under a CC-BY 3.0 licence, whose conditions is that the */
510 : /* following copyright notice is retained as well as the link to : */
511 : /* http://www.movable-type.co.uk/scripts/latlong-vincenty-direct.html */
512 : /******************************************************************************/
513 :
514 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
515 : /* Vincenty Direct Solution of Geodesics on the Ellipsoid (c) Chris Veness 2005-2012 */
516 : /* */
517 : /* from: Vincenty direct formula - T Vincenty, "Direct and Inverse Solutions of Geodesics on the */
518 : /* Ellipsoid with application of nested equations", Survey Review, vol XXII no 176, 1975 */
519 : /* http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf */
520 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
521 :
522 2 : std::pair <double,double> IRISDataset::GeodesicCalculation(float fLat, float fLon, float fAngle, float fDist, float fEquatorialRadius, float fPolarRadius, float fFlattening)
523 : {
524 2 : std::pair <double,double> oOutput;
525 2 : double dfAlpha1 = DEG2RAD * fAngle;
526 2 : double dfSinAlpha1 = sin(dfAlpha1);
527 2 : double dfCosAlpha1 = cos(dfAlpha1);
528 :
529 2 : double dfTanU1 = (1-fFlattening) * tan(fLat*DEG2RAD);
530 2 : double dfCosU1 = 1 / sqrt((1 + dfTanU1*dfTanU1));
531 2 : double dfSinU1 = dfTanU1*dfCosU1;
532 :
533 2 : double dfSigma1 = atan2(dfTanU1, dfCosAlpha1);
534 2 : double dfSinAlpha = dfCosU1 * dfSinAlpha1;
535 2 : double dfCosSqAlpha = 1 - dfSinAlpha*dfSinAlpha;
536 2 : double dfUSq = dfCosSqAlpha * (fEquatorialRadius*fEquatorialRadius - fPolarRadius*fPolarRadius) / (fPolarRadius*fPolarRadius);
537 2 : double dfA = 1 + dfUSq/16384*(4096+dfUSq*(-768+dfUSq*(320-175*dfUSq)));
538 2 : double dfB = dfUSq/1024 * (256+dfUSq*(-128+dfUSq*(74-47*dfUSq)));
539 :
540 2 : double dfSigma = fDist / (fPolarRadius*dfA);
541 2 : double dfSigmaP = 2*M_PI;
542 :
543 : double dfSinSigma;
544 : double dfCosSigma;
545 : double dfCos2SigmaM;
546 : double dfDeltaSigma;
547 :
548 6 : while (fabs(dfSigma-dfSigmaP) > 1e-12) {
549 2 : dfCos2SigmaM = cos(2*dfSigma1 + dfSigma);
550 2 : dfSinSigma = sin(dfSigma);
551 2 : dfCosSigma = cos(dfSigma);
552 : dfDeltaSigma = dfB*dfSinSigma*(dfCos2SigmaM+dfB/4*(dfCosSigma*(-1+2*dfCos2SigmaM*dfCos2SigmaM)-
553 2 : dfB/6*dfCos2SigmaM*(-3+4*dfSinSigma*dfSinSigma)*(-3+4*dfCos2SigmaM*dfCos2SigmaM)));
554 2 : dfSigmaP = dfSigma;
555 2 : dfSigma = fDist / (fPolarRadius*dfA) + dfDeltaSigma;
556 : }
557 :
558 2 : double dfTmp = dfSinU1*dfSinSigma - dfCosU1*dfCosSigma*dfCosAlpha1;
559 : double dfLat2 = atan2(dfSinU1*dfCosSigma + dfCosU1*dfSinSigma*dfCosAlpha1,
560 2 : (1-fFlattening)*sqrt(dfSinAlpha*dfSinAlpha + dfTmp*dfTmp));
561 2 : double dfLambda = atan2(dfSinSigma*dfSinAlpha1, dfCosU1*dfCosSigma - dfSinU1*dfSinSigma*dfCosAlpha1);
562 2 : double dfC = fFlattening/16*dfCosSqAlpha*(4+fFlattening*(4-3*dfCosSqAlpha));
563 : double dfL = dfLambda - (1-dfC) * fFlattening * dfSinAlpha *
564 2 : (dfSigma + dfC*dfSinSigma*(dfCos2SigmaM+dfC*dfCosSigma*(-1+2*dfCos2SigmaM*dfCos2SigmaM)));
565 2 : double dfLon2 = fLon*DEG2RAD+dfL;
566 2 : if (dfLon2 > M_PI)
567 0 : dfLon2 = dfLon2 - 2*M_PI;
568 2 : if (dfLon2 < -1*M_PI)
569 0 : dfLon2 = dfLon2 + 2*M_PI;
570 2 : oOutput.first = dfLon2*RAD2DEG;
571 2 : oOutput.second = dfLat2*RAD2DEG;
572 :
573 2 : return oOutput;
574 : }
575 :
576 : /************************************************************************/
577 : /* GetGeoTransform() */
578 : /************************************************************************/
579 :
580 1 : CPLErr IRISDataset::GetGeoTransform( double * padfTransform )
581 :
582 : {
583 1 : if (!bHasLoadedProjection)
584 0 : LoadProjection();
585 1 : memcpy( padfTransform, adfGeoTransform, sizeof(double)*6 );
586 1 : return CE_None;
587 : }
588 :
589 : /************************************************************************/
590 : /* GetProjectionRef() */
591 : /************************************************************************/
592 :
593 1 : const char *IRISDataset::GetProjectionRef(){
594 1 : if (!bHasLoadedProjection)
595 1 : LoadProjection();
596 1 : return pszSRS_WKT;
597 : }
598 :
599 : /************************************************************************/
600 : /* Identify() */
601 : /************************************************************************/
602 :
603 11380 : int IRISDataset::Identify( GDALOpenInfo * poOpenInfo )
604 :
605 : {
606 :
607 : /* -------------------------------------------------------------------- */
608 : /* Confirm that the file is an IRIS file */
609 : /* -------------------------------------------------------------------- */
610 : //Si no el posem, peta al fer el translate, quan s'obre Identify des de GDALIdentifyDriver
611 11380 : if( poOpenInfo->nHeaderBytes < 640 )
612 11258 : return FALSE;
613 :
614 :
615 122 : short nId1 = CPL_LSBSINT16PTR(poOpenInfo->pabyHeader);
616 122 : short nId2 = CPL_LSBSINT16PTR(poOpenInfo->pabyHeader+12);
617 122 : unsigned short nType = CPL_LSBUINT16PTR (poOpenInfo->pabyHeader+24);
618 :
619 : /*Check if the two headers are 27 (product hdr) & 26 (product configuration), and the product type is in the range 1 -> 34*/
620 122 : if( !(nId1 == 27 && nId2 == 26 && nType > 0 && nType < 35) )
621 120 : return FALSE;
622 :
623 2 : return TRUE;
624 : }
625 :
626 : /************************************************************************/
627 : /* FillString() */
628 : /************************************************************************/
629 :
630 14 : static void FillString(char* szBuffer, size_t nBufferSize, void* pSrcBuffer)
631 : {
632 190 : for(size_t i = 0; i < nBufferSize - 1; i++)
633 176 : szBuffer[i] = ((char*)pSrcBuffer)[i];
634 14 : szBuffer[nBufferSize-1] = '\0';
635 14 : }
636 :
637 : /************************************************************************/
638 : /* Open() */
639 : /************************************************************************/
640 :
641 1316 : GDALDataset *IRISDataset::Open( GDALOpenInfo * poOpenInfo )
642 :
643 : {
644 1316 : if (!Identify(poOpenInfo))
645 1314 : return NULL;
646 : /* -------------------------------------------------------------------- */
647 : /* Confirm the requested access is supported. */
648 : /* -------------------------------------------------------------------- */
649 2 : if( poOpenInfo->eAccess == GA_Update )
650 : {
651 : CPLError( CE_Failure, CPLE_NotSupported,
652 : "The IRIS driver does not support update access to existing"
653 0 : " datasets.\n" );
654 0 : return NULL;
655 : }
656 :
657 : /* -------------------------------------------------------------------- */
658 : /* Create a corresponding GDALDataset. */
659 : /* -------------------------------------------------------------------- */
660 : IRISDataset *poDS;
661 :
662 2 : poDS = new IRISDataset();
663 :
664 2 : poDS->fp = VSIFOpenL( poOpenInfo->pszFilename, "rb" );
665 2 : if (poDS->fp == NULL)
666 : {
667 0 : delete poDS;
668 0 : return NULL;
669 : }
670 :
671 : /* -------------------------------------------------------------------- */
672 : /* Read the header. */
673 : /* -------------------------------------------------------------------- */
674 2 : VSIFReadL( poDS->abyHeader, 1, 640, poDS->fp );
675 2 : int nXSize = CPL_LSBSINT32PTR(poDS->abyHeader+100+12);
676 2 : int nYSize = CPL_LSBSINT32PTR(poDS->abyHeader+104+12);
677 2 : int nNumBands = CPL_LSBSINT32PTR(poDS->abyHeader+108+12);
678 :
679 2 : poDS->nRasterXSize = nXSize;
680 :
681 2 : poDS->nRasterYSize = nYSize;
682 2 : if (poDS->nRasterXSize <= 0 || poDS->nRasterYSize <= 0 )
683 : {
684 : CPLError( CE_Failure, CPLE_AppDefined,
685 : "Invalid dimensions : %d x %d",
686 0 : poDS->nRasterXSize, poDS->nRasterYSize);
687 0 : delete poDS;
688 0 : return NULL;
689 : }
690 :
691 2 : if( !GDALCheckBandCount(nNumBands, TRUE) )
692 : {
693 0 : delete poDS;
694 0 : return NULL;
695 : }
696 :
697 : /* -------------------------------------------------------------------- */
698 : /* Create band information objects. */
699 : /* -------------------------------------------------------------------- */
700 8 : for (int iBandNum = 1; iBandNum <= nNumBands; iBandNum++) {
701 2 : poDS->SetBand( iBandNum, new IRISRasterBand( poDS, iBandNum ));
702 2 : poDS->GetRasterBand(iBandNum)->SetNoDataValue(-9999);
703 : }
704 :
705 : /* -------------------------------------------------------------------- */
706 : /* Setting the Metadata */
707 : /* -------------------------------------------------------------------- */
708 : //See point 3.2.26 at page 3.12 of the manual
709 2 : poDS->nProductCode = CPL_LSBUINT16PTR (poDS->abyHeader+12+12);
710 2 : poDS->SetMetadataItem( "PRODUCT_ID", CPLString().Printf("%d", poDS->nProductCode ));
711 2 : if( poDS->nProductCode >= ARRAY_ELEMENT_COUNT(poDS->aszProductNames) )
712 : {
713 0 : delete poDS;
714 0 : return NULL;
715 : }
716 :
717 2 : poDS->SetMetadataItem( "PRODUCT",poDS->aszProductNames[poDS->nProductCode]);
718 :
719 2 : poDS->nDataTypeCode = CPL_LSBUINT16PTR (poDS->abyHeader+130+12);
720 2 : if( poDS->nDataTypeCode >= ARRAY_ELEMENT_COUNT(poDS->aszDataTypeCodes) )
721 : {
722 0 : delete poDS;
723 0 : return NULL;
724 : }
725 2 : poDS->SetMetadataItem( "DATA_TYPE_CODE",poDS->aszDataTypeCodes[poDS->nDataTypeCode]);
726 :
727 2 : if( poDS->nDataTypeCode >= ARRAY_ELEMENT_COUNT(poDS->aszDataTypes) )
728 : {
729 0 : delete poDS;
730 0 : return NULL;
731 : }
732 2 : poDS->SetMetadataItem( "DATA_TYPE",poDS->aszDataTypes[poDS->nDataTypeCode]);
733 :
734 2 : unsigned short nDataTypeInputCode = CPL_LSBUINT16PTR (poDS->abyHeader+144+12);
735 2 : if( nDataTypeInputCode >= ARRAY_ELEMENT_COUNT(poDS->aszDataTypeCodes) )
736 : {
737 0 : delete poDS;
738 0 : return NULL;
739 : }
740 2 : poDS->SetMetadataItem( "DATA_TYPE_INPUT_CODE",poDS->aszDataTypeCodes[nDataTypeInputCode]);
741 :
742 2 : unsigned short nDataTypeInput = CPL_LSBUINT16PTR (poDS->abyHeader+144+12);
743 2 : if( nDataTypeInput >= ARRAY_ELEMENT_COUNT(poDS->aszDataTypes) )
744 : {
745 0 : delete poDS;
746 0 : return NULL;
747 : }
748 2 : poDS->SetMetadataItem( "DATA_TYPE_INPUT",poDS->aszDataTypes[nDataTypeInput]);
749 :
750 2 : poDS->nProjectionCode = * (unsigned char *) (poDS->abyHeader+146+12);
751 2 : if( poDS->nProjectionCode >= ARRAY_ELEMENT_COUNT(poDS->aszProjections) )
752 : {
753 0 : delete poDS;
754 0 : return NULL;
755 : }
756 :
757 : ////TIMES
758 2 : int nSeconds = CPL_LSBSINT32PTR(poDS->abyHeader+20+12);
759 :
760 2 : int nHour = (nSeconds - (nSeconds%3600)) /3600;
761 2 : int nMinute = ((nSeconds - nHour * 3600) - (nSeconds - nHour * 3600)%60)/ 60;
762 2 : int nSecond = nSeconds - nHour * 3600 - nMinute * 60;
763 :
764 2 : short nYear = CPL_LSBSINT16PTR(poDS->abyHeader+26+12);
765 2 : short nMonth = CPL_LSBSINT16PTR(poDS->abyHeader+28+12);
766 2 : short nDay = CPL_LSBSINT16PTR(poDS->abyHeader+30+12);
767 :
768 2 : poDS->SetMetadataItem( "TIME_PRODUCT_GENERATED", CPLString().Printf("%d-%02d-%02d %02d:%02d:%02d", nYear, nMonth, nDay, nHour, nMinute, nSecond ) );
769 :
770 :
771 2 : nSeconds = CPL_LSBSINT32PTR(poDS->abyHeader+32+12);
772 :
773 2 : nHour = (nSeconds - (nSeconds%3600)) /3600;
774 2 : nMinute = ((nSeconds - nHour * 3600) - (nSeconds - nHour * 3600)%60)/ 60;
775 2 : nSecond = nSeconds - nHour * 3600 - nMinute * 60;
776 :
777 2 : nYear = CPL_LSBSINT16PTR(poDS->abyHeader+26+12);
778 2 : nMonth = CPL_LSBSINT16PTR(poDS->abyHeader+28+12);
779 2 : nDay = CPL_LSBSINT16PTR(poDS->abyHeader+30+12);
780 :
781 4 : poDS->SetMetadataItem( "TIME_INPUT_INGEST_SWEEP", CPLString().Printf("%d-%02d-%02d %02d:%02d:%02d", nYear, nMonth, nDay, nHour, nMinute, nSecond ) );
782 :
783 : ///Site and task information
784 :
785 2 : char szSiteName[17] = ""; //Must have one extra char for string end!
786 2 : char szVersionName[9] = "";
787 :
788 2 : FillString(szSiteName, sizeof(szSiteName), poDS->abyHeader+320+12);
789 2 : FillString(szVersionName, sizeof(szVersionName), poDS->abyHeader+16+320+12);
790 2 : poDS->SetMetadataItem( "PRODUCT_SITE_NAME",szSiteName);
791 2 : poDS->SetMetadataItem( "PRODUCT_SITE_IRIS_VERSION",szVersionName);
792 :
793 2 : FillString(szSiteName, sizeof(szSiteName), poDS->abyHeader+90+320+12);
794 2 : FillString(szVersionName, sizeof(szVersionName), poDS->abyHeader+24+320+12);
795 2 : poDS->SetMetadataItem( "INGEST_SITE_NAME",szSiteName);
796 2 : poDS->SetMetadataItem( "INGEST_SITE_IRIS_VERSION",szVersionName);
797 :
798 2 : FillString(szSiteName, sizeof(szSiteName), poDS->abyHeader+74+320+12);
799 2 : poDS->SetMetadataItem( "INGEST_HARDWARE_NAME",szSiteName);
800 :
801 2 : char szConfigFile[13] = "";
802 2 : FillString(szConfigFile, sizeof(szConfigFile), poDS->abyHeader+62+12);
803 2 : poDS->SetMetadataItem( "PRODUCT_CONFIGURATION_NAME",szConfigFile);
804 :
805 2 : char szTaskName[13] = "";
806 2 : FillString(szTaskName, sizeof(szTaskName), poDS->abyHeader+74+12);
807 2 : poDS->SetMetadataItem( "TASK_NAME",szTaskName);
808 :
809 2 : short nRadarHeight = CPL_LSBSINT16PTR(poDS->abyHeader+284+320+12);
810 4 : poDS->SetMetadataItem( "RADAR_HEIGHT",CPLString().Printf("%d m",nRadarHeight));
811 2 : short nGroundHeight = CPL_LSBSINT16PTR(poDS->abyHeader+118+320+12);
812 4 : poDS->SetMetadataItem( "GROUND_HEIGHT",CPLString().Printf("%d m",nRadarHeight-nGroundHeight)); //Ground height over the sea level
813 :
814 2 : unsigned short nFlags = CPL_LSBUINT16PTR (poDS->abyHeader+86+12);
815 : //Get eleventh bit
816 2 : nFlags=nFlags<<4;
817 2 : nFlags=nFlags>>15;
818 2 : if (nFlags == 1){
819 1 : poDS->SetMetadataItem( "COMPOSITED_PRODUCT","YES");
820 1 : unsigned int compositedMask = CPL_LSBUINT32PTR (poDS->abyHeader+232+320+12);
821 1 : poDS->SetMetadataItem( "COMPOSITED_PRODUCT_MASK",CPLString().Printf("0x%08x",compositedMask));
822 : } else{
823 1 : poDS->SetMetadataItem( "COMPOSITED_PRODUCT","NO");
824 : }
825 :
826 : //Wave values
827 2 : poDS->SetMetadataItem( "PRF",CPLString().Printf("%d Hz",CPL_LSBSINT32PTR(poDS->abyHeader+120+320+12)));
828 4 : poDS->SetMetadataItem( "WAVELENGTH",CPLString().Printf("%4.2f cm",(float) CPL_LSBSINT32PTR(poDS->abyHeader+148+320+12)/100));
829 2 : unsigned short nPolarizationType = CPL_LSBUINT16PTR (poDS->abyHeader+172+320+12);
830 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
831 2 : if (nPolarizationType == 1)
832 0 : fNyquist = fNyquist * 2;
833 2 : else if(nPolarizationType == 2)
834 0 : fNyquist = fNyquist * 3;
835 2 : else if(nPolarizationType == 3)
836 0 : fNyquist = fNyquist * 4;
837 2 : poDS->fNyquistVelocity = fNyquist;
838 2 : poDS->SetMetadataItem( "NYQUIST_VELOCITY",CPLString().Printf("%.2f m/s",fNyquist));
839 :
840 : ///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
841 : //See point 3.2.25 at page 3.12 of the manual
842 2 : if (EQUAL(poDS->aszProductNames[poDS->nProductCode],"PPI")){
843 : //Degrees = 360 * (Binary Angle)*2^N
844 : //float fElevation = 360 * float((CPL_LSBUINT16PTR (poDS->abyHeader+164+12))) / 65536;
845 1 : float fElevation = 360 * float((CPL_LSBSINT16PTR (poDS->abyHeader+164+12))) / 65536;
846 :
847 1 : poDS->SetMetadataItem( "PPI_ELEVATION_ANGLE",CPLString().Printf("%f",fElevation));
848 1 : if (EQUAL(poDS->aszDataTypeCodes[poDS->nDataTypeCode],"dBZ"))
849 1 : poDS->SetMetadataItem( "DATA_TYPE_UNITS","dBZ");
850 : else
851 0 : poDS->SetMetadataItem( "DATA_TYPE_UNITS","m/s");
852 : //See point 3.2.2 at page 3.2 of the manual
853 1 : } else if (EQUAL(poDS->aszProductNames[poDS->nProductCode],"CAPPI")){
854 1 : float fElevation = ((float) CPL_LSBSINT32PTR(poDS->abyHeader+4+164+12))/100;
855 1 : poDS->SetMetadataItem( "CAPPI_HEIGHT",CPLString().Printf("%.1f m",fElevation));
856 1 : float fAzimuthSmoothingForShear = 360 * float((CPL_LSBUINT16PTR (poDS->abyHeader+10+164+12))) / 65536;
857 2 : poDS->SetMetadataItem( "AZIMUTH_SMOOTHING_FOR_SHEAR" ,CPLString().Printf("%.1f", fAzimuthSmoothingForShear));
858 1 : unsigned int nMaxAgeVVPCorrection = CPL_LSBUINT32PTR (poDS->abyHeader+24+164+12);
859 2 : poDS->SetMetadataItem( "MAX_AGE_FOR_SHEAR_VVP_CORRECTION" ,CPLString().Printf("%d s", nMaxAgeVVPCorrection));
860 1 : if (EQUAL(poDS->aszDataTypeCodes[poDS->nDataTypeCode],"dBZ"))
861 1 : poDS->SetMetadataItem( "DATA_TYPE_UNITS","dBZ");
862 : else
863 0 : poDS->SetMetadataItem( "DATA_TYPE_UNITS","m/s");
864 : //See point 3.2.32 at page 3.19 of the manual
865 0 : } else if (EQUAL(poDS->aszProductNames[poDS->nProductCode],"RAIN1") || EQUAL(poDS->aszProductNames[poDS->nProductCode],"RAINN")){
866 0 : short nNumProducts = CPL_LSBSINT16PTR(poDS->abyHeader+170+320+12);
867 0 : poDS->SetMetadataItem( "NUM_FILES_USED",CPLString().Printf("%d",nNumProducts));
868 :
869 0 : float fMinZAcum= (float)((CPL_LSBUINT32PTR (poDS->abyHeader+164+12))-32768)/1000;
870 0 : poDS->SetMetadataItem( "MINIMUM_Z_TO_ACUMULATE",CPLString().Printf("%f",fMinZAcum));
871 :
872 0 : unsigned short nSecondsOfAccumulation = CPL_LSBUINT16PTR (poDS->abyHeader+6+164+12);
873 0 : poDS->SetMetadataItem( "SECONDS_OF_ACCUMULATION",CPLString().Printf("%d s",nSecondsOfAccumulation));
874 :
875 0 : unsigned int nSpanInputFiles = CPL_LSBUINT32PTR (poDS->abyHeader+24+164+12);
876 0 : poDS->SetMetadataItem( "SPAN_OF_INPUT_FILES",CPLString().Printf("%d s",nSpanInputFiles));
877 0 : poDS->SetMetadataItem( "DATA_TYPE_UNITS","mm");
878 :
879 0 : char szInputProductName[13] = "";
880 0 : for(int k=0; k<12;k++)
881 0 : szInputProductName[k] = * (char *) (poDS->abyHeader+k+12+164+12);
882 0 : poDS->SetMetadataItem( "INPUT_PRODUCT_NAME",CPLString().Printf("%s",szInputProductName));
883 :
884 0 : if (EQUAL(poDS->aszProductNames[poDS->nProductCode],"RAINN"))
885 0 : poDS->SetMetadataItem( "NUM_HOURS_ACCUMULATE",CPLString().Printf("%d",CPL_LSBUINT16PTR (poDS->abyHeader+10+164+12)));
886 :
887 : //See point 3.2.73 at page 3.36 of the manual
888 0 : } else if (EQUAL(poDS->aszProductNames[poDS->nProductCode],"VIL")){
889 0 : float fBottomHeigthInterval = (float) CPL_LSBSINT32PTR(poDS->abyHeader+4+164+12) / 100;
890 0 : poDS->SetMetadataItem( "BOTTOM_OF_HEIGTH_INTERVAL",CPLString().Printf("%.1f m",fBottomHeigthInterval));
891 0 : float fTopHeigthInterval = (float) CPL_LSBSINT32PTR(poDS->abyHeader+8+164+12) / 100;
892 0 : poDS->SetMetadataItem( "TOP_OF_HEIGTH_INTERVAL",CPLString().Printf("%.1f m",fTopHeigthInterval));
893 0 : poDS->SetMetadataItem( "VIL_DENSITY_NOT_AVAILABLE_VALUE","-1");
894 0 : poDS->SetMetadataItem( "DATA_TYPE_UNITS","mm");
895 : //See point 3.2.68 at page 3.36 of the manual
896 0 : } else if (EQUAL(poDS->aszProductNames[poDS->nProductCode],"TOPS")){
897 0 : float fZThreshold = (float) CPL_LSBSINT16PTR(poDS->abyHeader+4+164+12) / 16;
898 0 : poDS->SetMetadataItem( "Z_THRESHOLD",CPLString().Printf("%.1f dBZ",fZThreshold));
899 0 : poDS->SetMetadataItem( "ECHO_TOPS_NOT_AVAILABLE_VALUE","-1");
900 0 : poDS->SetMetadataItem( "DATA_TYPE_UNITS","km");
901 : //See point 3.2.20 at page 3.10 of the manual
902 0 : } else if (EQUAL(poDS->aszProductNames[poDS->nProductCode],"MAX")){
903 0 : float fBottomInterval = (float) CPL_LSBSINT32PTR(poDS->abyHeader+4+164+12) / 100;
904 0 : poDS->SetMetadataItem( "BOTTOM_OF_INTERVAL",CPLString().Printf("%.1f m",fBottomInterval));
905 0 : float fTopInterval = (float) CPL_LSBSINT32PTR(poDS->abyHeader+8+164+12) / 100;
906 0 : poDS->SetMetadataItem( "TOP_OF_INTERVAL",CPLString().Printf("%.1f m",fTopInterval));
907 0 : int nNumPixelsSidePanels = CPL_LSBSINT32PTR(poDS->abyHeader+12+164+12);
908 0 : poDS->SetMetadataItem( "NUM_PIXELS_SIDE_PANELS",CPLString().Printf("%d",nNumPixelsSidePanels));
909 0 : short nHorizontalSmootherSidePanels = CPL_LSBSINT16PTR(poDS->abyHeader+16+164+12);
910 0 : poDS->SetMetadataItem( "HORIZONTAL_SMOOTHER_SIDE_PANELS",CPLString().Printf("%d",nHorizontalSmootherSidePanels));
911 0 : short nVerticalSmootherSidePanels = CPL_LSBSINT16PTR(poDS->abyHeader+18+164+12);
912 0 : poDS->SetMetadataItem( "VERTICAL_SMOOTHER_SIDE_PANELS",CPLString().Printf("%d",nVerticalSmootherSidePanels));
913 : }
914 :
915 : /* -------------------------------------------------------------------- */
916 : /* Initialize any PAM information. */
917 : /* -------------------------------------------------------------------- */
918 2 : poDS->SetDescription( poOpenInfo->pszFilename );
919 2 : poDS->TryLoadXML();
920 :
921 : /* -------------------------------------------------------------------- */
922 : /* Check for overviews. */
923 : /* -------------------------------------------------------------------- */
924 2 : poDS->oOvManager.Initialize( poDS, poOpenInfo->pszFilename );
925 :
926 2 : return( poDS );
927 : }
928 :
929 : /************************************************************************/
930 : /* GDALRegister_IRIS() */
931 : /************************************************************************/
932 :
933 582 : void GDALRegister_IRIS()
934 :
935 : {
936 : GDALDriver *poDriver;
937 :
938 582 : if( GDALGetDriverByName( "IRIS" ) == NULL )
939 : {
940 561 : poDriver = new GDALDriver();
941 :
942 561 : poDriver->SetDescription( "IRIS" );
943 : poDriver->SetMetadataItem( GDAL_DMD_LONGNAME,
944 561 : "IRIS data (.PPI, .CAPPi etc)" );
945 : poDriver->SetMetadataItem( GDAL_DMD_HELPTOPIC,
946 561 : "frmt_various.html#IRIS" );
947 561 : poDriver->SetMetadataItem( GDAL_DMD_EXTENSION, "ppi" );
948 561 : poDriver->SetMetadataItem( GDAL_DCAP_VIRTUALIO, "YES" );
949 :
950 561 : poDriver->pfnOpen = IRISDataset::Open;
951 561 : poDriver->pfnIdentify = IRISDataset::Identify;
952 :
953 561 : GetGDALDriverManager()->RegisterDriver( poDriver );
954 : }
955 582 : }
|
