Pileus Git - ~andy/rsl/blob - carpi.c

   1 /*
   2     NASA/TRMM, Code 910.1.
   3     This is the TRMM Office Radar Software Library.
   4     Copyright (C) 1996, 1997
   5             Mike Kolander
   6             Space Applications Corporation
   7             Vienna, Virginia
   8
   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.
  13
  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.
  18
  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.
  22 */
  23 /* Mike Kolander
  24  * Space Applications Corporation
  25  * NASA/Goddard 910.1
  26  */
  27
  28 #include <stdio.h>
  29 #include <math.h>
  30 #include <stdlib.h>
  31 #include <string.h>
  32
  33 #include "rsl.h"
  34
  35 #define RAD2DEG 57.29578 /* radian to degree conversion */
  36 #define MAXRAYS 512      /* loop safety valve when traversing a sweep */
  37
  38 extern int radar_verbose_flag;
  39
  40
  41 /*************************************************************/
  42 /*                                                           */
  43 /*                     RSL_free_carpi                        */
  44 /*                                                           */
  45 /*************************************************************/
  46 void RSL_free_carpi(Carpi *carpi)
  47 {
  48         /* Frees memory allocated to a carpi structure, and associated
  49                  pointer and data arrays.
  50   */
  51
  52         if (carpi != NULL)
  53         {
  54           if (carpi->data != NULL)
  55                 {
  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. */
  59                 }
  60                 free(carpi);   /* Free the carpi structure. */
  61         }
  62 }
  63
  64 /*************************************************************/
  65 /*                                                           */
  66 /*                     RSL_new_carpi                         */
  67 /*                                                           */
  68 /*************************************************************/
  69 Carpi *RSL_new_carpi(int nrows, int ncols)
  70 {
  71         /* Allocates memory for a carpi structure, and associated pointer
  72                  and data arrays.
  73         */
  74         Carpi *c;
  75         Carpi_value *data;
  76         int row;
  77
  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;
  90
  91         return(c);
  92 }
  93
  94 /*************************************************************/
  95 /*                                                           */
  96 /*                     RSL_volume_to_carpi                   */
  97 /*                                                           */
  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)
 102 /*
 103  * Creates a carpi from a volume structure.
 104  *
 105  *   +--------------------+   ^
 106  *   |                    |   |
 107  *   |                    |   |
 108  *   |  cell size dx,dy   |   ny
 109  *   |                    |   |
 110  *   |                    |   |
 111  *   +--------------------+   V
 112  * lat,lon
 113  *   <-------- nx -------->
 114  *
 115  *  Radar centered at radar_x, radar_y.
 116  */
 117
 118 {
 119   Cappi *cappi;
 120   Carpi *carpi;
 121
 122   if (v == NULL) return NULL;
 123   cappi = RSL_cappi_at_h(v, h, grnd_r);
 124   cappi->lat = lat;
 125   cappi->lon = lon;
 126   cappi->interp_method = 0;
 127
 128   carpi = RSL_cappi_to_carpi(cappi, dx, dy, lat, lon, nx, ny,
 129                                                                                                                  radar_x, radar_y);
 130   if (carpi == NULL) return NULL;
 131   RSL_free_cappi(cappi);
 132   return carpi;
 133 }
 134
 135 /*************************************************************/
 136 /*                                                           */
 137 /*                      RSL_find_rng_azm                     */
 138 /*                                                           */
 139 /*************************************************************/
 140 void RSL_find_rng_azm(float *r, float *ang, float x, float y)
 141 {
 142 /* Convert Cartesian coords (x,y) to polar (rng,angle); 0<angle<360 */
 143         *r = sqrt((double)x*x + (double)y*y);
 144
 145         if (x != 0.0)
 146         {
 147                 *ang = RAD2DEG*atan((double)y/x);
 148                 if (x > 0.0)
 149                   *ang = 90.0 - *ang;
 150                 else
 151                   *ang = 270.0 - *ang;
 152         }
 153         else   /* x is 0.0 */
 154         {
 155                 if (y >= 0.0) *ang = 0.0;
 156                 else *ang = 180.0;
 157         }
 158 }
 159
 160 /*************************************************************/
 161 /*                                                           */
 162 /*                     RSL_cappi_to_carpi                    */
 163 /*                                                           */
 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.
 172           ********/
 173 {
 174         Carpi *carpi;
 175         Ray *first_ray;
 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*/
 180
 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. */
 185
 186         if (radar_verbose_flag) fprintf(stderr,"\nCreating carpi...\n");
 187
 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;
 196         carpi->dx = dx;
 197         carpi->dy = dy;
 198         carpi->nx = nx;
 199         carpi->ny = ny;
 200         carpi->radar_x = radar_x;
 201         carpi->radar_y = radar_y;
 202         carpi->height = cappi->height;
 203         carpi->lat = lat;
 204         carpi->lon = lon;
 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;
 210
 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)
 217           scy = 1;
 218         else if ((dy - gate_size) <= 1.0)
 219           scy = 2;
 220         else
 221           scy = 3;
 222
 223         if ((dx - gate_size) < 0.0)
 224           scx = 1;
 225         else if ((dx - gate_size) <= 1.0)
 226           scx = 2;
 227         else
 228           scx = 3;
 229
 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;
 239
 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);
 243
 244         /* For each cell (row,col) comprising the carpi...*/
 245         for (row=0; row<ny; row++)
 246         {
 247                 for (col=0; col<nx; col++)
 248                 {
 249                         /* For each subcell (there are scx*scy subcells per cell)
 250                                  of carpi cell (row,col)...*/
 251                         for (n=0; n<scy; n++)
 252                         {
 253                                 /* Convert carpi coords (row,col) to km coords (x,y) relative to
 254                                          radar. */
 255                                 y = (row + n/(float)scy - radar_y)*dy;
 256                                 for (m=0; m<scx; m++)
 257                                 {
 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 */
 263                                         if (m == 0)
 264                                         {
 265                                           if (n == 0)
 266                                                         if ((rng - cell_diag) > carpi_max_rng)
 267                                                         {
 268                                                                 /* Totality of carpi cell lies outside data range. */
 269                                                                 cell = (float) BADVAL;
 270                                                                 goto escape;
 271                                                         }
 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=... */
 277
 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.*/
 281                         cell = 0.0;
 282                         for (j=0; j<scy; j++)
 283                         {
 284                           for (k=0; k<scx; k++)
 285                                   if ((subcell[j][k] <= BADVAL) && (subcell[j][k] >= NOECHO))
 286                                     continue;
 287                                         else  /* valid data value. */
 288                                         {
 289                                                 cell = cell + subcell[j][k];
 290                                                 valid_subcells++;
 291                                         }
 292                         }  /* end for (j=0; j<scy; j++) */
 293                         if (valid_subcells != 0) cell = cell/valid_subcells;
 294                         else cell = (float) BADVAL;
 295
 296 escape:
 297                         carpi->data[row][col] = (Carpi_value) carpi->invf(cell);
 298                 } /* end for (col=... */
 299         }  /* end for (row=...  */
 300
 301         return(carpi);
 302 }