3 This is the TRMM Office Radar Software Library.
4 Copyright (C) 1996, 1997
6 Space Applications Corporation
9 This library is free software; you can redistribute it and/or
10 modify it under the terms of the GNU Library General Public
11 License as published by the Free Software Foundation; either
12 version 2 of the License, or (at your option) any later version.
14 This library is distributed in the hope that it will be useful,
15 but WITHOUT ANY WARRANTY; without even the implied warranty of
16 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 Library General Public License for more details.
19 You should have received a copy of the GNU Library General Public
20 License along with this library; if not, write to the Free
21 Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
24 /******************************************************************
25 * Opens a Mcgill format radar data file, and reads and
26 * decompresses the contents.
28 * All Mcgill data structures referenced herein are
29 * defined in the file "mcgill.h".
31 * There exist 3 functions in this file which are
32 * accessible to application routines:
33 * -----------------------------------------------------
34 * mcgFile_t *mcgFileOpen(int *code, char *filename);
35 * Opens and verifies the content of a Mcgill format
36 * data file, reads in the file header record, and
37 * initializes the mcgFile_t structure.
38 * -------------------------------------------------------
39 * int mcgFileClose(mcgFile_t *file);
40 * Closes the Mcgill data file.
41 * -------------------------------------------------------
42 * int mcgRayBuild(mcgRay_t *ray, mcgFile_t *file);
43 * Returns in the mcgRay_t structure the reflectivity values from
44 * one ray of data from a Mcgill format radar data file.
45 * Successive calls return data from successive rays found
50 *******************************************************************/
62 /* Fixed elevation angles for Mcgill sweeps. */
63 static float elev_angle[3][24] =
65 /* PAFB volume scan format types 3 and 4 (12 sweeps) */
67 0.6, 1.0, 1.6, 2.5, 3.6, 5.2, 7.5, 10.5, 14.7, 20.2,
70 /* PAFB volume scan format types 1 and 2 (24 sweeps) */
72 0.6, 0.79, 1.01, 1.27, 1.57, 1.93, 2.35, 2.84,
73 3.41, 4.08, 4.86, 5.78, 6.85, 8.09, 9.55, 11.23,
74 13.18, 15.42, 17.98, 20.89, 24.16, 27.78, 31.73, 35.97
76 /* Sao Paulo radar?? (20 sweeps) */
78 1.0, 1.3, 1.6, 2.0, 2.4, 2.9, 3.6, 4.2, 5.1,
79 6.1, 7.3, 8.7, 10.3, 12.2, 14.4, 16.9, 19.8,
84 /* Fixed number_of_bins for Mcgill sweeps. */
85 static int num_bins_normal[24] =
87 /* PAFB volume scan formats 1 and 3 (normal mode) */
88 240, 180, 180, 180, 180, 180, 180, 180, 180, 180,
89 180, 180, 180, 180, 180, 180, 180, 180, 180, 180,
92 static int num_bins_compressed[24] =
94 /* PAFB volume scan formats 2 and 4 (compressed mode) */
95 180, 120, 120, 120, 120, 120, 120, 120, 120, 120,
96 120, 120, 120, 120, 120, 120, 120, 120, 120, 120,
100 FILE *uncompress_pipe (FILE *fp);
102 /**********************************************************************/
103 mcgFile_t *mcgFileOpen(int *code, char *filename)
104 /**********************************************************************/
106 /* Open a Mcgill format radar data file, and read and verify
107 the content of the first (header) record from the file.
108 This function returns one of the following coded integer values:
109 MCG_OK: Normal return.
110 MCG_OPEN_FILE_ERR: Error occurred while attempting to open
112 MCG_EOF: Unexpected End_of_File occurred reading from file.
113 MCG_READ_ERR: Error occurred while reading from file.
114 MCG_FORMAT_ERR: Format error encountered reading from file.
115 SYS_NO_SPACE: Error attempting to allocate memory space.
118 char buffer[MCG_RECORD];
121 /* Allocate space for the mcgFile structure. */
122 if ((file = (mcgFile_t *)malloc(sizeof(mcgFile_t))) == NULL)
124 *code = SYS_NO_SPACE;
128 /* Open Mcgill data file for reading */
129 if ((file->fp = fopen(filename, "r")) == NULL)
131 *code = MCG_OPEN_FILE_ERR;
134 file->fp = uncompress_pipe(file->fp); /* Transparently, use gunzip. */
135 /* Read first (header) record from data file into buffer */
136 if (fread(buffer, sizeof(char), MCG_RECORD, file->fp) < MCG_RECORD)
138 *code = MCG_FORMAT_ERR;
142 /* Move header values from buffer to the mcgHeader structure.*/
143 memcpy(&file->head, buffer, sizeof(mcgHeader_t));
145 /* Check the leading header characters for the site ID string. */
146 if ((strncmp((char *)&file->head, "PAFB RADAR VOLUME SCAN AT", 25) == 0) ||
147 (strncmp((char *)&file->head, "P A B RADAR VOLUME SCAN AT", 27) == 0))
149 file->site = MCG_PAFB_1_2;
151 /* Check if file from Sao Paulo.
152 else if (strncmp((char *)&file->head, SAOP_ID_STRING) == 0)
154 file->site = MCG_SAOP;
157 else /* Can't identify site/format. */
159 *code = MCG_FORMAT_ERR;
163 /* Check validity of some file header values. If
164 not valid, we assume the file is garbled beyond legibility. */
165 if ((file->head.vol_scan_format < 1) || (file->head.vol_scan_format > 4)
166 || (file->head.day < 1) || (file->head.day > 31)
167 || (file->head.month < 1) || (file->head.month > 12))
169 /* Invalid/unexpected file format */
170 *code = MCG_FORMAT_ERR;
174 /* Set the proper number_of_bins, depending on scan format. */
175 switch (file->head.vol_scan_format)
179 file->num_bins = num_bins_normal;
183 file->num_bins = num_bins_compressed;
189 /* I've never seen scan_formats 3 or 4 from PAFB, but I make
190 allowance for them here. */
191 if (file->site == MCG_PAFB_1_2)
193 if ((file->head.vol_scan_format == 3) ||
194 (file->head.vol_scan_format == 4))
196 file->site = MCG_PAFB_3_4;
200 /* If a Mcgill Csp block is contained in this file, it's next.
201 Skip past it. It's useless for our purposes. */
202 if (file->head.csp_rec == 1)
204 /* if (fseek(file->fp, MCG_CSP, SEEK_CUR) != 0)
206 *code = MCG_READ_ERR;
210 /* Allocate buffer space for the csp block, read the block
211 from the disk file into the buffer, then free the buffer.
212 An ugly way to discard the block, but the fseek solution
213 above won't work with stdin capabilities.
215 if ((csp_buffer=(char *)malloc(MCG_CSP)) == NULL)
217 *code = SYS_NO_SPACE;
220 if (fread(csp_buffer, sizeof(char), MCG_CSP, file->fp) < MCG_CSP)
222 *code = MCG_READ_ERR;
228 /* File is open and properly initialized. */
234 int rsl_pclose( FILE *stream);
235 /**********************************************************************/
236 int mcgFileClose(mcgFile_t *file)
237 /**********************************************************************/
239 /* Close a Mcgill format radar data file. Returns one of the
240 following coded integer values:
241 MCG_OK: File closed successfully.
242 MCG_CLOSE_FILE_ERR: System error occurred during file close
247 if (rsl_pclose(file->fp) == 0)
253 return(MCG_CLOSE_FILE_ERR);
258 /************************************************************************/
259 int mcgRecordRead(mcgRecord_t *record, mcgFile_t *file)
260 /************************************************************************/
262 /* Reads a Mcgill logical record (2048 bytes) from the data file
263 into the mcgRecord_t structure. Returns one of the following coded
265 MCG_OK: Successfully read a record.
266 MCG_EOF: End_of_File condition occurred while reading a record.
267 MCG_READ_ERR: Read error occurred.
270 /* Read data from file directly into the record structure. */
271 if (fread(record, sizeof(char), MCG_RECORD, file->fp) < MCG_RECORD)
276 return(MCG_READ_ERR);
284 /************************************************************************/
285 mcgSegmentID mcgSegmentKeyIdentify(int key)
286 /************************************************************************/
288 /* Decodes the first byte from a Mcgill file segment.
289 The first byte of the segment identifies the type of the segment.
290 Mcgill segment types are data, elevation, and end_of_data.
292 if ((key > 0) && (key < 16))
293 return(MCG_DATA_SEG);
294 else if ((key > 31) && (key < 63))
295 return(MCG_ELEV_SEG);
299 return(MCG_UNDEFINED_SEG);
304 /*************************************************************************/
305 int mcgRayBuild(mcgRay_t *ray, mcgFile_t *file)
306 /*************************************************************************/
308 /* Returns in the ray structure the reflectivity values from one
309 ray of data from a Mcgill format radar data file. Successive calls
310 to this routine return data from successive rays found in the file.
311 This function returns the following coded integer values:
312 MCG_OK : Ray structure successfully filled with data values.
313 MCG_EOD: Empty ray structure, No more data in file.
314 MCG_EOF: Premature End_of_File condition encountered.
315 MCG_READ_ERR: System error occurred while reading from file.
316 MCG_FORMAT_ERR: Format error encountered while reading from file.
319 /* This function is typically called about 9000 times while reading a
320 Mcgill data file. The following static variables must retain their
321 values between successive calls. */
322 static int seg_num=MCG_MAX_SEG_NUM;
323 static int eod_found = FALSE;
324 static int sweep_num=0;
325 static float elev=0.0;
328 mcgSegmentID seg_type;
329 static mcgRecord_t record;
332 /* If we've previously found the end_of_data (eod) data segment in
333 the last Mcgill record, we're done. Do this for robustness, in case
334 the calling application routine screws up. */
340 /* Initialize the ray data values. */
343 ray->sweep_num = sweep_num;
344 memset(ray->data, 0, 240);
345 ray->num_bins = file->num_bins[sweep_num];
347 /* Loop to fill the ray bins with data values from consecutive
348 Mcgill data segments. */
351 /* Read in from the Mcgill file a new record, if necessary. */
352 if (seg_num == MCG_MAX_SEG_NUM)
354 if ((n=mcgRecordRead(&record, file)) < MCG_OK)
358 /* record_empty = FALSE; */
363 /* The first byte of segment is used to identify segment type. */
364 seg_type = mcgSegmentKeyIdentify(record.segment[seg_num][0]);
368 /* Determine the azimuth angle for this data segment */
369 azm = record.segment[seg_num][1] * 64.0 +
370 record.segment[seg_num][2] - 1.0;
371 /* Check if this data segment contains data belonging to
372 the current ray. If so, add the segment bin values to
373 the ray structure. */
374 if ((azm == ray->azm) || (ray->azm == EMPTY))
377 /* Compute the ray bin base_address for the 16 bin values
378 from this data segment. */
379 base = (record.segment[seg_num][0] - 1)*16;
380 /* Move the 16 bin values from segment to ray structure */
383 ray->data[base + j] = record.segment[seg_num][3+j];
386 else /* This data segment belongs to the next ray. Return the
387 current ray. We'll start with this segment the next
388 time mcgRayBuild() is called. */
395 sweep_num = record.segment[seg_num][0] - 31 - 1;
396 elev = elev_angle[file->site][sweep_num];
397 /* If ray structure is not empty, return it, since
398 we've found a Mcgill elevation segment, which initiates
400 if (ray->azm != EMPTY)
403 if (seg_num == MCG_MAX_SEG_NUM)
404 /* record_empty = TRUE; */
407 else /* Ray structure is empty. Continue on to read the
408 next Mcgill segment.*/
411 ray->sweep_num = sweep_num;
417 if (ray->azm == EMPTY)
419 else /* Return the ray structure we've filled. */
424 return(MCG_FORMAT_ERR);
429 } /* end while (1) */