render_pcp.c

Go to the documentation of this file.

/* Converts point precipitation data in struct pdata to hrap grid.  The
   20 closest stations are precalculated for each grid point - this
   speeds computations for data sets with many precipitation points.
   If there are no good precipitation points for a grid point, then
   a recalculation is made using all precipitation points. 1/R**2
   interpolation is used.  If requested, the final interpolation
   is scaled using seasonal isohyets. The grid is saved as a disk
   file and used later for a raster or vector (HRAP) plot. */
 
#include "prototypes_new.h"
 
void render_pcp(int pcpn_day,int pcpn_time,int pcpn_time_step)
 
{
 
extern int isom;
extern Widget top_level;
extern int pcp_in_use[500];
extern Cursor watch_cursor;
extern int method;
extern struct pdata pdata[10];
extern struct station station[3000];
int i,j,h,hh,time_pos,htotal;
double distance,dist1,dist2,dist,value;
double temp;
extern struct hrap_grid *hrap_grid;
extern struct pcp *pcp;
extern Window window;
extern Display *display;
extern int max_stations;
int totals[5];
 
 
printf("in render_pcp\n");
 
for(i=0;i<5;i++)
        totals[i]=0;
 
XDefineCursor(display,window,watch_cursor); 
XmUpdateDisplay(top_level);
 
if(pcpn_time_step==0)
              time_pos=pcpn_time;
 
else
              time_pos=4;
 
for(i=0;i<hrap_grid->maxi;i++) {
 
for(j=0;j<hrap_grid->maxj;j++) {
 
           if(hrap_grid->owner[i][j]==-1) {
           
                                 pcp->value[i][j]=0;
                                 continue;
                                 
                                 }
 
           value=0.0;
           distance=0.0;
           htotal=0;
 
           for(h=0;h<30;h++) {
 
               hh=hrap_grid->gage[i][j].index[h];
 
               if(pdata[pcpn_day].stn[hh].frain[time_pos].data < 0 ){
            /*   printf("no data\n");*/
                                                 continue;
};
 
               if(pdata[pcpn_day].stn[hh].frain[time_pos].qual != 0 &&
                  pdata[pcpn_day].stn[hh].frain[time_pos].qual != 8 &&
                  pdata[pcpn_day].stn[hh].frain[time_pos].qual != 6 &&
                  pdata[pcpn_day].stn[hh].frain[time_pos].qual != 3 &&
                  pdata[pcpn_day].stn[hh].frain[time_pos].qual != 4 &&
          pdata[pcpn_day].stn[hh].frain[time_pos].qual != 2)
                          continue;
 
               dist1=(double)(hrap_grid->coord[i][j].x-
                      station[hh].x);    
               dist2=(double)(hrap_grid->coord[i][j].y-
                      station[hh].y);  
 
               dist=pow(dist1,2)+pow(dist2,2);
 
               if(dist < .00001)
                        dist=.00001;
               
               dist=1/dist;
 
               temp=pdata[pcpn_day].stn[hh].frain[time_pos].data*dist;
 
               if(method==2 && station[hh].isoh[isom] > 0) 
                     temp=temp*(float)hrap_grid->isoh[isom][i][j]/
                          (station[hh].isoh[isom] * 25.4);
 
               if(method==2 && station[hh].isoh[isom] <= 0)
              
                                 continue;
 
               value=value+temp;
 
               distance=distance+dist;                    
 
               htotal++;
               
               if(htotal ==10)
                            break;
 
         }
                
 
           if(htotal < 4) {
          
           value=0.0;
           distance=0.0;
           htotal=0;
 
           for(h=0;h<max_stations;h++) {
 
               if(pdata[pcpn_day].stn[h].frain[time_pos].data < 0 )
                                                 continue;
 
               if(pdata[pcpn_day].stn[h].frain[time_pos].qual != 0 &&
                  pdata[pcpn_day].stn[h].frain[time_pos].qual != 8 &&
                  pdata[pcpn_day].stn[h].frain[time_pos].qual != 6 &&
                  pdata[pcpn_day].stn[h].frain[time_pos].qual != 3 &&
                  pdata[pcpn_day].stn[h].frain[time_pos].qual != 4 &&
          pdata[pcpn_day].stn[h].frain[time_pos].qual != 2)
                            continue;
 
               dist1=(double)(hrap_grid->coord[i][j].x-
                      station[h].x);    
               dist2=(double)(hrap_grid->coord[i][j].y-
                      station[h].y);  
 
               dist=pow(dist1,2)+pow(dist2,2);
 
               if(dist < .00001)
                        dist=.00001;
               
               dist=1/dist;
 
               temp=pdata[pcpn_day].stn[h].frain[time_pos].data*dist;
 
               if(method==2 && station[h].isoh[isom] <= 0)
                                 continue;
 
               if(method==2 && station[h].isoh[isom] > 0)
                     temp=temp*(float)hrap_grid->isoh[isom][i][j]/
                          (station[h].isoh[isom]*25.4);
 
               value=value+temp;
 
               distance=distance+dist;                    
 
               htotal++;
 
         }
 
}
 
           if(htotal == 0) 
               pcp->value[i][j]=0;
 
           else
               pcp->value[i][j]=(int)(value/distance*100);
  
           if(pcp->value[i][j] < .01)
                    pcp->value[i][j]=0;
     
     }
 
}
 
if(pcpn_time_step==0)
              time_pos=pcpn_day*4+3-pcpn_time;
 
else
              time_pos=40+pcpn_day;
 
pcp_in_use[time_pos]=1;
 
write_file("pcp",time_pos,pcp);
 
}