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 * Space Applications Corporation
35 #define RAD2DEG 57.29578 /* radian to degree conversion */
36 #define MAXRAYS 512 /* loop safety valve when traversing a sweep */
38 extern int radar_verbose_flag;
41 /*************************************************************/
45 /*************************************************************/
46 void RSL_free_carpi(Carpi *carpi)
48 /* Frees memory allocated to a carpi structure, and associated
49 pointer and data arrays.
54 if (carpi->data != NULL)
56 if (carpi->data[0] != NULL)
57 free(carpi->data[0]); /* Free the 2D data array. */
58 free(carpi->data); /* Free the vector of pointers. */
60 free(carpi); /* Free the carpi structure. */
64 /*************************************************************/
68 /*************************************************************/
69 Carpi *RSL_new_carpi(int nrows, int ncols)
71 /* Allocates memory for a carpi structure, and associated pointer
78 /* Allocate a carpi. */
79 c = (Carpi *)calloc(1, sizeof(Carpi));
80 if (c == NULL) perror("RSL_new_carpi");
81 /* Allocate a vector of 'nrows' pointers. */
82 c->data = (Carpi_value **)calloc(nrows, sizeof(Carpi_value *));
83 if (c->data == NULL) perror("RSL_new_carpi");
84 /* Allocate a 2 dimensional array for the actual data values. */
85 data = (Carpi_value *)calloc(nrows*ncols, sizeof(Carpi_value));
86 if (data == NULL) perror("RSL_new_carpi");
87 /* Fill all the 'nrows' pointer slots created above. */
88 for (row=0; row<nrows; row++)
89 c->data[row] = data + row*ncols;
94 /*************************************************************/
96 /* RSL_volume_to_carpi */
98 /*************************************************************/
99 Carpi *RSL_volume_to_carpi(Volume *v, float h, float grnd_r,
100 float dx, float dy, int nx, int ny,
101 int radar_x, int radar_y, float lat, float lon)
103 * Creates a carpi from a volume structure.
105 * +--------------------+ ^
108 * | cell size dx,dy | ny
111 * +--------------------+ V
113 * <-------- nx -------->
115 * Radar centered at radar_x, radar_y.
122 if (v == NULL) return NULL;
123 cappi = RSL_cappi_at_h(v, h, grnd_r);
126 cappi->interp_method = 0;
128 carpi = RSL_cappi_to_carpi(cappi, dx, dy, lat, lon, nx, ny,
130 if (carpi == NULL) return NULL;
131 RSL_free_cappi(cappi);
135 /*************************************************************/
137 /* RSL_find_rng_azm */
139 /*************************************************************/
140 void RSL_find_rng_azm(float *r, float *ang, float x, float y)
142 /* Convert Cartesian coords (x,y) to polar (rng,angle); 0<angle<360 */
143 *r = sqrt((double)x*x + (double)y*y);
147 *ang = RAD2DEG*atan((double)y/x);
155 if (y >= 0.0) *ang = 0.0;
160 /*************************************************************/
162 /* RSL_cappi_to_carpi */
164 /*************************************************************/
165 Carpi *RSL_cappi_to_carpi(Cappi *cappi, float dx, float dy, float lat,
166 float lon, int nx, int ny, int radar_x, int radar_y)
167 /****** Simple and straightforward algorithm:
168 Divide each of the nx*ny carpi cells into scx*scy subcells.
169 Find the data value for each subcell from the cappi rays.
170 Average the subcell data values over a cell to obtain the cell value.
171 Store the cell value into the 2_D carpi array.
176 int row, col, j, k, m, n, scx, scy, valid_subcells;
177 float x, y, rng, azm, cell, cell_diag, gate_size;
178 float carpi_max_rng, radar_max_rng;
179 float subcell[3][3]; /* Maximum of 9 subcells per carpi cell*/
181 if (cappi == NULL) return NULL;
182 if (cappi->sweep == NULL) return NULL;
183 first_ray = RSL_get_first_ray_of_sweep(cappi->sweep);
184 if (first_ray == NULL) return NULL; /* No data. */
186 if (radar_verbose_flag) fprintf(stderr,"\nCreating carpi...\n");
188 /* Allcate space for a carpi, and fill in its values. */
189 carpi = RSL_new_carpi(ny, nx);
190 carpi->month = cappi->month;
191 carpi->day = cappi->day;
192 carpi->year = cappi->year;
193 carpi->hour = cappi->hour;
194 carpi->minute = cappi->minute;
195 carpi->sec = cappi->sec;
200 carpi->radar_x = radar_x;
201 carpi->radar_y = radar_y;
202 carpi->height = cappi->height;
205 strncpy(carpi->radar_type, cappi->radar_type, sizeof(cappi->radar_type));
206 carpi->field_type = cappi->field_type;
207 carpi->interp_method = cappi->interp_method;
208 carpi->f = first_ray->h.f;
209 carpi->invf = first_ray->h.invf;
211 gate_size = first_ray->h.gate_size / 1000.0; /* km */
212 cell_diag = sqrt(dx*dx + dy*dy);
213 /* How many subcells per carpi cell? The larger the carpi cell relative
214 to gate_size, the more subcells we want. Will have scx*scy subcells
215 per carpi cell. Note below that the max no. of subcells is 9 . */
216 if ((dy - gate_size) < 0.0)
218 else if ((dy - gate_size) <= 1.0)
223 if ((dx - gate_size) < 0.0)
225 else if ((dx - gate_size) <= 1.0)
230 /* Max range of the radar data */
231 radar_max_rng = first_ray->h.nbins * gate_size; /* km */
232 /* Max range of desired carpi is max of x and y directions. */
233 carpi_max_rng = nx / 2.0 * dx; /* km */
234 if ( (ny / 2.0 * dy) > carpi_max_rng )
235 carpi_max_rng = ny / 2.0 * dy;
236 /* carpi_max_rng is smaller of (radar_max_rng, carpi_max_rng) */
237 if (radar_max_rng < carpi_max_rng)
238 carpi_max_rng = radar_max_rng;
240 if (radar_verbose_flag)
241 fprintf(stderr,"carpi_max_rng:%.1f(km) beam_width:%.1f gate_size:%d(m)\n",
242 carpi_max_rng, cappi->sweep->h.beam_width, first_ray->h.gate_size);
244 /* For each cell (row,col) comprising the carpi...*/
245 for (row=0; row<ny; row++)
247 for (col=0; col<nx; col++)
249 /* For each subcell (there are scx*scy subcells per cell)
250 of carpi cell (row,col)...*/
251 for (n=0; n<scy; n++)
253 /* Convert carpi coords (row,col) to km coords (x,y) relative to
255 y = (row + n/(float)scy - radar_y)*dy;
256 for (m=0; m<scx; m++)
258 x = (col + m/(float)scx - radar_x)*dx;
259 /* Find range and azimuth of subcell (relative to radar) */
260 RSL_find_rng_azm(&rng, &azm, x, y);
261 /* Check if carpi cell outside of data range by noting
262 location of lower left corner of cell */
266 if ((rng - cell_diag) > carpi_max_rng)
268 /* Totality of carpi cell lies outside data range. */
269 cell = (float) BADVAL;
272 } /* end if (m == 0) */
273 /* cell lies in data range. Get data value from cappi. */
274 subcell[n][m] = RSL_get_value_from_cappi(cappi,rng,azm);
275 } /* end for (m=... */
276 } /* end for (n=... */
278 /* All subcell values have now been determined. Average them
279 to determine the cell value. */
280 valid_subcells = 0; /* number subcells having valid data values.*/
282 for (j=0; j<scy; j++)
284 for (k=0; k<scx; k++)
285 if ((subcell[j][k] <= BADVAL) && (subcell[j][k] >= NOECHO))
287 else /* valid data value. */
289 cell = cell + subcell[j][k];
292 } /* end for (j=0; j<scy; j++) */
293 if (valid_subcells != 0) cell = cell/valid_subcells;
294 else cell = (float) BADVAL;
297 carpi->data[row][col] = (Carpi_value) carpi->invf(cell);
298 } /* end for (col=... */
299 } /* end for (row=... */