make_nexrad.c

Go to the documentation of this file.

#define _POSIX_SOURCE
char fbuf[100];
char directory[100];
char database_directory[100];
char nexrad_directory[100];
 
#include "prototypes.h"
#define TRUE 1
#define FALSE 0
 
char *convertJulianDate();
 
main(int argc,char **argv)
   
{
 extern int errno;
 char tarbuf[1000];
 int irain;
 struct nexrad nexrad;
 int l;
 int maxi,maxj;
 int nrapx[131][131],nrapy[131][131];
 char dbuf[200];
 extern int errno;
 FILE *afosfile;
 FILE *fp;
 char ch, prevchar,pprevchar;
 int i,j,n,ihr,iminute,len,k,kk;
 char clat[9],clon[9],height[7],year[3],month[3],day[3];
 char hour[3],minute[3],radar[4],chour[3];
 unsigned char ncharrows;
 char param[7],t1[7],t2[7],t3[7];
 char maxval[7],bias[7],error[7],jdate[7],jmin[7];
 char outname[81],sjunk[50],lradar[4];
 char cdate[7];
 div_t time;
 unsigned char data[300];
 short  nruns,nbytes;
 short xrun,level;
 unsigned char xchar;
 double lat,lon;
 struct HRAP hrap;
 struct HRAP irap;
 float cx,cy;
 unsigned char value[131][131];
 unsigned char cbytes1,cbytes2; 
 FILE *fpout;
 short int maxvalue;
 char fname[100];
 char *p;
 int imonth,iyear,ihour,iday;
 char ibuf[100];
 char dname[100];
 float xlat,nlat,xlon,nlon;
 float hmaxx,hmaxy,hminx,hminy;
 int imaxx,iminx,imaxy,iminy;
 char kbuf[120];
 float val[500][500];
 char *nexname[20];
 FILE *fq;
 float pw,pnexrad[255];
 int ii,jj;
 char name[100];
 int ier;
          
strcpy(nexrad_directory,"/usr/mapper/nexrad");
 
afosfile = fopen(argv[1],"rb");
 
if(afosfile==NULL) {
 
                    printf("could not open %s\n",argv[1]);
 
                    exit(1);
 
          }
 
fread(data,sizeof(char),300,afosfile);
 
for(i=0;i<250;i++) {
     
          n=memcmp(&data[i],"HDP",3);
          if(n==0)
              break;
 
                  }
 
if(n != 0) {
 
           printf("could not find HDP\n");
           fclose(afosfile);
           exit(1);
 
     }
 
rewind(afosfile);
     
data[i+6]=0;
 
memcpy(radar,&data[i+3],4);
 
printf("decoding of %s started\n",argv[1]);
  
for (;;) {
 
    n = fscanf(afosfile,"%c",&ch);
    if (n==EOF) {
 
         printf("Error 1 %s\n",fname);
         fclose(afosfile);
         exit(1);
 
       }
 
    if (ch == ' ' && prevchar=='1' && pprevchar == '8') 
               break;
 
    pprevchar=prevchar;
    prevchar=ch;
 
     }
 
   /* read header block */
 
    n = fscanf(afosfile,"%7s%8s%6s%2s%2s%2s%2s%2s",clat,clon,height,year,month,
           day,hour,minute); 
  
    printf("%s %s %s %s ",radar,clat,clon,height);
 
    if (n==EOF) {
 
         printf("Error 2 %s\n",fname);
         fclose(afosfile);
         exit(1);
 
       }
 
 n = fscanf(afosfile,"%6s%6s%6s%6s",param,t1,t2,t3);
 if (n==EOF)  {
 
         printf("Error 3 %s\n",fname);
         fclose(afosfile);
         exit(1);
 
       }
 
 for (i=0;i<13;i++)
    fscanf(afosfile,"%6s",param);
 
 n = fscanf(afosfile,"%6s%6s%6s%6s%6s%6s%6s",maxval,bias,error,
       jdate,jmin,param,param);
 
 if (n==EOF){
 
         printf("Error 3 %s\n",fname);
         fclose(afosfile);
         exit(1);
 
       }
 
 strcpy(cdate,convertJulianDate((short)atoi(jdate)));
 
  /* determine product valid time, if not H+00 product quit */
 
 time = div(atoi(jmin),60);
 ihr = time.quot;
 iminute = time.rem;
 
 printf("%s %02d:%02d max=%6s bias=%6s err=%6s\n",
        cdate, ihr, iminute, maxval,bias,error);
 
 if (iminute > 10) {
 
                 printf("%s too old - aborting\n",argv[1]);
 
                 fclose(afosfile);
                 exit(1);
 
           }
 
 sprintf(chour,"%02d",ihr);
 
 fread(sjunk,sizeof(char),11,afosfile);
 
 maxvalue=0;
 
 n = fread(&ncharrows, sizeof(char),1,afosfile);
 
 for (i=0; i<ncharrows; i++)
 
     {
 
      n = fread(&cbytes1, sizeof(char),1,afosfile);
      n = fread(&cbytes2, sizeof(char),1,afosfile);
      nbytes=cbytes2;
      
      n = fread(data, sizeof(unsigned char), nbytes-2,afosfile);
 
      if(n != nbytes-2)
 
      {
        printf(" unexpected EOF encoutered -- ");
        printf(" number of bytes expected = %d",nbytes-2);
        printf(" number of bytes found =%d \n",n);
 
        fclose(afosfile); 
        exit(1);
 
      }
 
      nruns = (short)((float)(nbytes-2)/2.);
      kk=0;
 
      for (j=0;j<nruns; j++)
     {
      xchar = (unsigned char)data[j*2];
      xrun = (short) xchar;
      level = (short) data[1+j*2];       
 
          for(k=0;k<xrun;k++) {
 
                        if(level > 0) {
                   
                          if(i > 130 || kk > 130) {
 
                               printf("Error 3 %s\n",fname);
                               fclose(afosfile);
                               exit(1);
 
                                           }
 
                          value[130-i][kk]=(char)level;
 
                          if(value[130-i][kk] != 255 &&
                             value[130-i][kk] > maxvalue)
                                            maxvalue=value[130-i][kk] ;
                          
                               }
 
                        else
                             value[130-i][kk]=0;
                   
          kk++;
 
         }
 
    }
 
    }
 
  fclose(afosfile);
 
  iyear=atoi(year);
  imonth=atoi(month);
  iday=atoi(day);
  ihour=atoi(chour);
 
  for(i=0;i<3;i++)
             lradar[i]=tolower(radar[i]);
  
  lradar[3]=0;
  
  sprintf(outname,"%s/raw/%s.%02d-%02d-%02d-%02d",nexrad_directory,
          lradar,iyear,imonth,iday,ihour);
 
  fpout=fopen(outname,"wb");
  if(fpout==NULL) {
 
                  printf("could not open %s\n",outname);
 
                  exit(1);
 
        }
 
  fwrite(value,sizeof(char),131*131,fpout);
 
  fclose(fpout);
 
  printf("max value is %d\n",maxvalue);
 
 /* get max min hrap bins */
 
  iyear=atoi(year);
  imonth=atoi(month);
  iday=atoi(day);
  ihour=atoi(chour);
 
  sprintf(dname,"%s/%02d/%02d",nexrad_directory,iyear,imonth);
 
  errno=access(dname,R_OK);
 
  if(errno != 0) {
  
         sprintf(ibuf,"mkdir -p %s",dname);
         system(ibuf);
 
         sprintf(ibuf,"chmod -R 777 /usr/mapper/nexrad");
         system(ibuf);
                           
         }
     
 sprintf(dbuf,"/usr/mapper/NEXRAD_radars");
 
 fp=fopen(dbuf,"r");
 
 if(fp==NULL) {
 
              printf("could not open /usr/mapper/NEXRAD_radars\n");
              exit(1);
 
              }
 
 p=fgets(kbuf,80,fp);
 if(p==NULL) {
 
              printf("error in format /usr/mapper/NEXRAD_radars\n");
              exit(1);
 
           }
 
 p=pars_line(kbuf,"x_lat=",dbuf);
 
 if(p==NULL) {
 
              printf("no x_lat\n");
              exit(1);
              
              }
              
 xlat=atof(dbuf);
   
 p=pars_line(kbuf,"n_lat=",dbuf);
 
 if(p==NULL) {
 
              printf("no n_lat\n");
              exit(1);
              
              }
              
 nlat=atof(dbuf);
 
 p=pars_line(kbuf,"x_lon=",dbuf);
 
 if(p==NULL) {
 
              printf("no x_lon\n");
              exit(1);
              
              }
              
 xlon=atof(dbuf);
   
 p=pars_line(kbuf,"n_lon=",dbuf);
 
 if(p==NULL) {
 
              printf("no n_lon\n");
              exit(1);
              
              }
              
 nlon=atof(dbuf);
  
 
 /* initialize val */
 
 for(i=0;i<500;i++) {
 for(j=0;j<500;j++) {
 
           val[i][j]=-9999;
           
           }
           
           }
 
 /* initialize name */
 
 for(i=0;i<20;i++)
        nexname[i]=calloc(20,1);
 
/* initialize pnexrad */
 
 for(i=0;i<255;i++) {
 
      pw=(-6) +( (float)(i-1))*.125;
      pw=pw/10;
      pw=.0397 * pow(10,pw);
 
      if(i==0)
             pnexrad[i]=0;
 
      else
             pnexrad[i]=pw;
 
    }
 
 hmaxx=-1;hmaxy=-1;
 hminx=9999;hminy=9999;
  
 hrap=LatLongToHrap(xlat,xlon);
 
 if(hrap.x > hmaxx)
             hmaxx=hrap.x;
 
 
 if(hrap.x < hminx)
             hminx=hrap.x;
 
 
 if(hrap.y > hmaxy)
             hmaxy=hrap.y;
 
 
 if(hrap.y < hminy)
             hminy=hrap.y; 
 
 hrap=LatLongToHrap(xlat,nlon);
 
 if(hrap.x > hmaxx)
             hmaxx=hrap.x;
 
 
 if(hrap.x < hminx)
             hminx=hrap.x;
 
 
 if(hrap.y > hmaxy)
             hmaxy=hrap.y;
 
 
 if(hrap.y < hminy)
             hminy=hrap.y; 
 
 hrap=LatLongToHrap(nlat,xlon);
 
 if(hrap.x > hmaxx)
             hmaxx=hrap.x;
 
 
 if(hrap.x < hminx)
             hminx=hrap.x;
 
 
 if(hrap.y > hmaxy)
             hmaxy=hrap.y;
 
 
 if(hrap.y < hminy)
             hminy=hrap.y; 
 
 hrap=LatLongToHrap(nlat,nlon);
 
 if(hrap.x > hmaxx)
             hmaxx=hrap.x;
 
 
 if(hrap.x < hminx)
             hminx=hrap.x;
 
 
 if(hrap.y > hmaxy)
             hmaxy=hrap.y;
 
 
 if(hrap.y < hminy)
             hminy=hrap.y; 
 
 imaxx=hmaxx+.5;
 iminx=hminx+.5;
 
 imaxy=hmaxy+.5;
 iminy=hminy+.5;
 
 maxi=imaxx-iminx;
 maxj=imaxy-iminy;
   
 /* need times to get for each file */
 
 for(l=0;l<20;l++) {
 
     p=fgets(kbuf,80,fp);
     if(p==NULL) {
 
                 fclose(fp);
                 break;
 
           }
 
     p=pars_line(kbuf,"radar=",dbuf);
 
     if(p==NULL)
               continue;
                   
     strcpy(name,dbuf); 
 
     for(i=0;i<3;i++)
              name[i]=tolower(name[i]);
 
     name[3]=0;
 
     strcpy(nexname[l],name);
 
     p=pars_line(kbuf,"latitude=",dbuf);
 
     if(p==NULL)
               continue;
 
     lat=atof(p);
 
     p=pars_line(kbuf,"longitude=",dbuf);
 
     if(p==NULL)
               continue;
 
     lon=atof(p);
 
     if(lon > 0)
              lon=-lon; 
    
     /* get hrap coordinates of radar */
 
     hrap=LatLongToHrap(lat,-lon);
 
     cx=(int)hrap.x+1;
     cy=(int)hrap.y;
 
     for (i=0; i < 131; i++){     
     for (k=0; k < 131 ;k++){
 
 
        irap.x=cx + (float)(k - 65);
        irap.y=cy + (float)(i-65);
 
        nrapx[i][k]=irap.x;
        nrapy[i][k]=irap.y;        
 
         }
 
        }
 
     sprintf(fbuf,"/usr/mapper/nexrad/raw/%s.%02d-%02d-%02d-%02d",
             name,iyear,imonth,iday,ihour);
 
     fq=fopen(fbuf,"rb");
 
             /* what if missing ?? */
 
     if(fq==NULL) {
 
                        printf("could not find %s\n",fbuf);
                        continue;
 
              }
 
     fread(&nexrad,sizeof(char),131*131,fq);
 
     fclose(fq);
             
     for (i=0; i < 131; i++) {
     for (j=0; j < 131 ;j++) {
 
                     ii=nrapx[i][j]-iminx;
                     jj=nrapy[i][j]-iminy;
 
                     if(ii < 0 || ii > maxi ||
                        jj < 0 || jj > maxj)
                                       continue;
                     
                     if(nexrad.value[i][j] == 255)
                                             continue;
 
                     if(pnexrad[nexrad.value[i][j]] > 
                        val[ii][jj])
                        val[ii][jj]=pnexrad[nexrad.value[i][j]];
 
           }
 
       }
 
       }
 
 
/* now add up all numbers and divide by total contributions */
 
fclose(fp);
 
sprintf(fbuf,"/usr/mapper/nexrad/%02d/%02d/ngrid.%02d-%02d-%02d-%02d",
iyear,imonth,iyear,imonth,iday,ihour);
 
printf("fbuf is %s\n",fbuf);
 
fp=fopen(fbuf,"w");
 
fprintf(fp,"%d %d %d %d 1\n",iminx,iminy,maxi,maxj);
 
for(i=0;i<maxi;i++) {
 
for(j=0;j<maxj;j++) { 
                       
 
                    if(val[i][j] >=0)
                      irain=val[i][j]* 100;
                   
                    else
                      irain=-9999;
                   
                    fprintf(fp," %5d",irain);
 
          }
 
       }
 
fclose(fp);
 
sprintf(tarbuf,"/usr/local/bin/gzip -f %s",fbuf);
printf("tarbuf %s\n",tarbuf);
ier=system(tarbuf);
printf("ieris %d %d\n",ier,errno);
 
}
 
 
char *convertJulianDate(jdate)
      short jdate;
{
 char cdate[7];
 int days[]={31,28,31,30,31,30,31,31,30,31,30,31} ;
 int year,month,ndays,date,total,nyear,leap_year;
 total=0;
 for (year=1970;year<2010;year++)
    {
      leap_year=FALSE;
      if (year==1972 || year==1976 || year==1980 ||
           year==1984 || year==1988 || year==1992 || year==1996 ||
           year==2000 || year==2004 || year==2008)
          leap_year=TRUE;
     for (month=0;month<12;month++)
     {
      total = total+days[month];
      if (month==1 && leap_year==TRUE)
          total++;
      if (total >= jdate)
        {
         ndays=days[month];
         if (month==1 && leap_year==TRUE) ndays++;
         date=ndays - (total-jdate);
         month=month+1;
         nyear=year-1900;
         if (year > 2000) nyear=year-2000;
         sprintf(cdate,"%02d%02d%02d",month,date,nyear);
         return cdate;
        }
     }
    }
}