Daily_QC
get_hrap_coord.c
Go to the documentation of this file.
1 #define _POSIX_SOURCE
2 #include "prototypes_new.h"
3 
4 void get_hrap_coord(double maximum_latitude,double minimum_latitude,double center_longitude,int smonth,int emonth)
5 
6 {
7  extern int maxmin_used;
8  extern struct topo *topo;
9  extern char station_list_file[100];
10  extern char hrap_gage_file[100];
11  extern int isohyets_used;
12  extern struct station station[3000];
13  extern Widget scrollbar;
14  extern Widget top_level;
15  extern int posit;
16  extern int max_stations;
17  extern Widget s_text;
18  int i,k,kk,maxi,maxj,l,m,h;
19  struct HRAP hrap;
20  struct HRAP irap;
21  float cx,cy;
22  extern struct dval dval;
23  char dbuf[200];
24  extern struct hrap_grid *hrap_grid;
25  extern struct isoh *isoh;
26  extern struct maxmin *maxmin;
27  extern struct pcp *pcp;
28  extern struct pcp *spf;
29  extern struct pcp *tpf;
30  extern int pcp_in_use[500];
31  long r,s;
32  float conv=.0174,totlat,sorted[30];
33  float dist,dist1,dist2,distance,value;
34  float maxvalue,minvalue,mindistance,maxdistance;
35  float isohlat,isohlon;
36  int jj,latindex,lonindex;
37  int mini,minj;
38  int slider_size,scrollm,ivalue,increment,page_increment;
39  struct stat statbuf;
40  double omaximum_latitude,ominimum_latitude,ocenter_longitude;
41  int ier,newflag,mm,mer;
42  FILE *fz;
43  time_t ost_mtime;
44  char *p,kbuf[200];
45  struct tm *tm;
46  int ix,iy;
47  float lat,lon;
48  float elevation;
49  float maxminlat,maxminlon;
50 
51  hrap_grid=(struct hrap_grid *) calloc(1,sizeof(struct hrap_grid));
52  if(hrap_grid==NULL) {
53 
54  printf("could not allocate hrap_grid space\n");
55  exit(1);
56 
57  }
58 
59  pcp=(struct pcp *) calloc(1,sizeof(struct pcp));
60  if(pcp==NULL) {
61 
62  printf("could not allocate pcp space\n");
63  exit(1);
64 
65  }
66 
67  spf=(struct pcp *) calloc(1,sizeof(struct pcp));
68  if(spf==NULL) {
69 
70  printf("could not allocate pcp space\n");
71  exit(1);
72 
73  }
74 
75  tpf=(struct pcp *) calloc(1,sizeof(struct pcp));
76  if(tpf==NULL) {
77 
78  printf("could not allocate pcp space\n");
79  exit(1);
80 
81  }
82  for(i=0;i<500;i++)
83  pcp_in_use[i]=-1;
84 
85  maxi=-1;maxj=-1;mini=999999;minj=99999;
86 
87  strcpy(dbuf,"HRAP Coordinates\n");
88  posit=posit+strlen(dbuf);
89  XmTextInsert(s_text,posit,dbuf);
90 
91  XtVaGetValues(scrollbar,XmNmaximum,&scrollm,NULL);
92 
93  XmScrollBarGetValues(scrollbar,&ivalue,&slider_size,
94  &increment,&page_increment);
95 
96  XmScrollBarSetValues(scrollbar,scrollm-slider_size,
97  slider_size, increment,page_increment,True);
98 
99  XmUpdateDisplay(top_level);
100 
101  /* need longitude of 4 corners of display */
102 
103  totlat=maximum_latitude-minimum_latitude+2.0;
104 
105  hrap=LatLongToHrap((float)maximum_latitude,(float)center_longitude-7);
106 
107  cx=(int)hrap.x+1;
108  cy=(int)hrap.y;
109 
110  if(cx > maxi)
111  maxi=cx;
112 
113  if(cy > maxj)
114  maxj=cy;
115 
116  if(cx < mini)
117  mini=cx;
118 
119  if(cy < minj)
120  minj=cy;
121 
122  hrap=LatLongToHrap((float)maximum_latitude,(float)center_longitude+7);
123 
124  cx=(int)hrap.x+1;
125  cy=(int)hrap.y;
126 
127  if(cx > maxi)
128  maxi=cx;
129 
130  if(cy > maxj)
131  maxj=cy;
132 
133  if(cx < mini)
134  mini=cx;
135 
136  if(cy < minj)
137  minj=cy;
138 
139  hrap=LatLongToHrap((float)minimum_latitude,(float)center_longitude-7);
140  cx=(int)hrap.x+1;
141  cy=(int)hrap.y;
142 
143  if(cx > maxi)
144  maxi=cx;
145 
146  if(cy > maxj)
147  maxj=cy;
148 
149  if(cx < mini)
150  mini=cx;
151 
152  if(cy < minj)
153  minj=cy;
154 
155  hrap=LatLongToHrap((float)minimum_latitude,(float)center_longitude+7);
156 
157  cx=(int)hrap.x+1;
158  cy=(int)hrap.y;
159 
160  if(cx > maxi)
161  maxi=cx;
162 
163  if(cy > maxj)
164  maxj=cy;
165 
166  if(cx < mini)
167  mini=cx;
168 
169  if(cy < minj)
170  minj=cy;
171 
172 /*allocate necessary space */
173 
174 hrap_grid->hrap_minx=mini;
175 hrap_grid->hrap_miny=minj;
176 hrap_grid->maxi=maxi-mini;
177 hrap_grid->maxj=maxj-minj;
178 
179 maxi=hrap_grid->maxi;
180 maxj=hrap_grid->maxj;
181 
182 printf("maxi %d maxj %d\n",maxi,maxj);
183 
184 pcp->value=(short int **) calloc(maxi,sizeof(short int *));
185 if(pcp->value==NULL) {
186 
187  printf("could not allocate pcp space\n");
188  exit(1);
189 
190  }
191 
192 spf->value=(short int **) calloc(maxi,sizeof(short int *));
193 if(spf->value==NULL) {
194 
195  printf("could not allocate pcp space\n");
196  exit(1);
197 
198  }
199 
200 tpf->value=(short int **) calloc(maxi+1000,sizeof(short int *));
201 if(tpf->value==NULL) {
202 
203  printf("could not allocate pcp space\n");
204  exit(1);
205 
206  }
207 
208 for(i=0;i<maxi;i++) {
209 
210  pcp->value[i]=(short int *) calloc(maxj,sizeof(short int));
211 
212  if(pcp->value[i]==NULL) {
213 
214  printf("no memory for pcp array\n");
215  exit(1);
216 
217  }
218 
219  }
220 
221 for(i=0;i<maxi;i++) {
222 
223  spf->value[i]=(short int *) calloc(maxj,sizeof(short int));
224 
225  if(spf->value[i]==NULL) {
226 
227  printf("no memory for pcp array\n");
228  exit(1);
229 
230  }
231 
232  }
233 
234 for(i=0;i<maxi+1000;i++) {
235 
236  tpf->value[i]=(short int *) calloc(maxj,sizeof(short int));
237 
238  if(tpf->value[i]==NULL) {
239 
240  printf("no memory for pcp array\n");
241  exit(1);
242 
243  }
244 
245  }
246 
247 hrap_grid->coord=(struct coord**) calloc(maxi,sizeof(struct coord *));
248 hrap_grid->isoh=(short int***) calloc(maxi,sizeof(short int **));
249 hrap_grid->gage=(struct gage**) calloc(maxi,sizeof(struct gage *));
250 hrap_grid->owner=(short int**) calloc(maxi,sizeof(short int **));
251 hrap_grid->max=(short int***) calloc(maxi,sizeof(short int **));
252 hrap_grid->min=(short int***) calloc(maxi,sizeof(short int **));
253 hrap_grid->elev=(short int**) calloc(maxi,sizeof(short int **));
254 
255 if(hrap_grid->coord==NULL ||
256  hrap_grid->isoh==NULL ||
257  hrap_grid->gage==NULL ||
258  hrap_grid->owner==NULL ||
259  hrap_grid->max==NULL ||
260  hrap_grid->min==NULL ||
261  hrap_grid->elev==NULL ) {
262 
263  printf("no memory for hrap array\n");
264  exit(1);
265 
266  }
267 
268 for(k=0;k<12;k++) {
269 
270  ier=is_good(k,smonth,emonth);
271 
272  if(ier==-1)
273  continue;
274 
275  hrap_grid->isoh[k]=(short int **)calloc(maxi,sizeof(short int *));
276  hrap_grid->max[k]=(short int **)calloc(maxi,sizeof(short int *));
277  hrap_grid->min[k]=(short int **)calloc(maxi,sizeof(short int *));
278 
279  if(hrap_grid->isoh[k]==NULL ||
280  hrap_grid->max[k]==NULL ||
281  hrap_grid->min[k]==NULL) {
282 
283  printf("no memory for hrap array\n");
284  exit(1);
285 
286  }
287 
288  }
289 
290 for(i=0;i<maxi;i++) {
291 
292  hrap_grid->coord[i]=(struct coord *) calloc(maxj,sizeof(struct coord));
293  hrap_grid->gage[i]=(struct gage *) calloc(maxj,sizeof(struct gage));
294  hrap_grid->owner[i]=(short int *) calloc(maxj,sizeof(short int));
295  hrap_grid->elev[i]=(short int *) calloc(maxj,sizeof(short int));
296 
297  if(hrap_grid->coord[i]==NULL ||
298  hrap_grid->gage[i]== NULL ||
299  hrap_grid->owner[i]==NULL ||
300  hrap_grid->elev[i]==NULL) {
301 
302  printf("no memory for hrap array\n");
303  exit(1);
304 
305  }
306 
307 
308  for(k=0;k<12;k++) {
309 
310  ier=is_good(k,smonth,emonth);
311 
312  if(ier==-1)
313  continue;
314 
315  hrap_grid->isoh[k][i]=(short int *) calloc(maxj,sizeof(short int));
316  hrap_grid->max[k][i]=(short int *) calloc(maxj,sizeof(short int));
317  hrap_grid->min[k][i]=(short int *) calloc(maxj,sizeof(short int));
318 
319  if(hrap_grid->isoh[k][i]==NULL ||
320  hrap_grid->max[k][i]==NULL ||
321  hrap_grid->min[k][i]==NULL ) {
322 
323  printf("no memory for hrap array\n");
324  exit(1);
325 
326  }
327 
328  }
329 
330  }
331 
332 newflag=0;
333 
334 if(hrap_gage_file[0]==0)
335  strcpy(hrap_gage_file,"map_qc.grid");
336 
337 ier=stat(hrap_gage_file,&statbuf);
338 
339 /* get creation time of station_list file */
340 
341 mer=stat(station_list_file,&statbuf);
342 
343 /* grid file does not exist */
344 
345 if(ier==-1)
346  newflag=2;
347 
348 else {
349 
350  fz=fopen(hrap_gage_file,"r");
351 
352  p=fgets(dbuf,100,fz);
353 
354  ier=sscanf(dbuf,"%ld %lf %lf %lf",&ost_mtime,&omaximum_latitude,
355  &ominimum_latitude,&ocenter_longitude);
356 
357  tm=gmtime(&statbuf.st_mtime);
358 
359  if(ost_mtime != statbuf.st_mtime ||
360  omaximum_latitude != maximum_latitude ||
361  ominimum_latitude != minimum_latitude ||
362  ocenter_longitude != center_longitude)
363  newflag=1;
364 
365  fclose(fz);
366 
367  }
368 
369 if(newflag==0) {
370 
371  fz=fopen(hrap_gage_file,"r");
372  p=fgets(dbuf,100,fz);
373 
374  }
375 
376 else if(newflag==1)
377  fz=fopen(hrap_gage_file,"w+");
378 
379 else
380  fz=fopen(hrap_gage_file,"w");
381 
382 if(newflag!=0) {
383 
384  strcpy(dbuf,"Gage to grid recalculation\n");
385  posit=posit+strlen(dbuf);
386  XmTextInsert(s_text,posit,dbuf);
387 
388  XtVaGetValues(scrollbar,XmNmaximum,&scrollm,NULL);
389 
390  XmScrollBarGetValues(scrollbar,&ivalue,&slider_size,
391  &increment,&page_increment);
392 
393  XmScrollBarSetValues(scrollbar,scrollm-slider_size,
394  slider_size, increment,page_increment,True);
395 
396  XmUpdateDisplay(top_level);
397 
398  fprintf(fz,"%ld %f %f %f\n",statbuf.st_mtime,maximum_latitude,
399  minimum_latitude,center_longitude);
400 
401  }
402 
403 for (i=0; i < hrap_grid->maxi; i++){
404 
405 for (k=0; k < hrap_grid->maxj ;k++){
406 
407  irap.x=mini+i;
408  irap.y=minj+k;
409 
410  hrap=HrapToLatLong(irap);
411 
412  hrap_grid->coord[i][k].lat=hrap.y;
413  hrap_grid->coord[i][k].lon=hrap.x;
414 
415  r=dval.a * cos(hrap.y*conv)/(1+sin(hrap.y*conv))
416  * cos((hrap.x-dval.lo-90)*conv) + dval.xo +.5;
417 
418  hrap_grid->coord[i][k].x=r;
419 
420  s=dval.a * cos(hrap.y*conv)/(1+sin(hrap.y*conv))
421  * sin((hrap.x-dval.lo-90)*conv) + dval.yo + .5;
422 
423  hrap_grid->coord[i][k].y=s;
424 
425  iy=(topo->max_lat-hrap.y)/topo->delta_lat;
426  ix=(topo->max_lon-hrap.x)/topo->delta_lon;
427 
428  if(topo==NULL || ix <= 0 || iy <= 0 || ix >= topo->maxi-1 || iy >=
429  topo->maxj-1) {
430 
431  hrap_grid->elev[i][k]=-1;
432 
433  }
434 
435 
436  else {
437 
438  distance=0.0;
439  value=0.0;
440 
441  for(kk=ix;kk<=ix+1;kk++) {
442  for(jj=iy;jj<=iy+1;jj++) {
443 
444  lat=topo->max_lat-(float)jj * topo->delta_lat;
445  lon=topo->max_lon-(float)kk * topo->delta_lon;
446 
447  dist2=(hrap.x-lon)*cos((lat+hrap.y)/2*conv);
448 
449  dist1=hrap.y-lat;
450 
451  dist=pow(dist1,2)+pow(dist2,2);
452 
453  if(dist<=0.000001)
454  dist=0.000001;
455 
456  dist=1/dist;
457 
458  value=value+dist*(float)topo->value[kk][jj];
459  distance=distance+dist;
460 
461  }
462  }
463 
464  elevation=(value*32.1)/distance;
465  hrap_grid->elev[i][k]=(int)elevation;
466 
467  }
468 
469  /* we have lat long of point we need seasonal isohyet */
470 
471  for(l=0;l<30;l++)
472  sorted[l]=9999999;
473 
474  if(newflag==0){
475 
476  m=0;
477  fread(kbuf,sizeof(char),150,fz);
478  kbuf[150]=0;
479  for(mm=0;mm<30;mm++) {
480 
481  hrap_grid->gage[i][k].index[mm]=atoi(&kbuf[m]);
482 
483  m=m+5;
484 
485 
486  }
487 
488 
489  }
490 
491  else {
492 
493  for(m=0;m<max_stations;m++) {
494 
495  dist1=hrap.y-station[m].lat;
496  dist2=(hrap.x-station[m].lon)*cos((hrap.y+station[m].lat)/2*conv);
497 
498  dist=pow(dist1,2)+pow(dist2,2);
499 
500  for(l=0;l<30;l++) {
501 
502  if(dist < sorted[l]) {
503 
504  for(h=29; h > l; h--) {
505 
506  sorted[h]=sorted[h-1];
507 
508  hrap_grid->gage[i][k].index[h]=
509  hrap_grid->gage[i][k].index[h-1];
510 
511  }
512 
513  sorted[l]=dist;
514 
515  hrap_grid->gage[i][k].index[l]=m;
516 
517  break;
518 
519  }
520 
521 
522  }
523 
524  }
525 
526  for(l=0;l<30;l++)
527  fprintf(fz,"%04d ",hrap_grid->gage[i][k].index[l]);
528 
529  }
530 
531  for(h=0;h<12;h++) {
532 
533  ier=is_good(h,smonth,emonth);
534 
535  if(ier==-1)
536  continue;
537 
538  hrap_grid->isoh[h][i][k]=-1;
539  hrap_grid->max[h][i][k]=-9999;
540  hrap_grid->min[h][i][k]=-9999;
541 
542  if(isohyets_used!=-0){
543 
544  latindex=(isoh->max_lat-hrap.y)/isoh->delta_lat;
545  lonindex=(isoh->max_lon-hrap.x)/isoh->delta_lon;
546 
547  /* next get average isohyet for hrap point */
548 
549  if(latindex >= 0 && lonindex >= 0 &&
550  latindex < isoh->maxj-1 && lonindex < isoh->maxi-1) {
551 
552  distance=0.0;
553  value=0.0;
554 
555  for(kk=latindex;kk<=latindex+1;kk++) {
556 
557  isohlat=isoh->max_lat - (float)kk*isoh->delta_lat;
558 
559  for(jj=lonindex;jj<=lonindex+1;jj++) {
560 
561  isohlon=isoh->max_lon - (float)jj*isoh->delta_lon;
562 
563  dist1=isohlat-irap.y;
564  dist2=(isohlon-irap.x)*cos((isohlat+irap.y)/2*conv);
565 
566  dist=pow(dist1,2)+pow(dist2,2);
567 
568  if(isoh->value[h][jj][kk] <=0)
569  continue;
570 
571  value=value+dist*isoh->value[h][jj][kk];
572  distance=distance+dist;
573 
574  }
575  }
576 
577  if(distance != 0)
578  hrap_grid->isoh[h][i][k]=value/distance;
579 
580 
581  }
582 
583  }
584 
585  if(maxmin_used!=0) {
586 
587  latindex=(maxmin->max_lat-hrap.y)/maxmin->delta_lat;
588  lonindex=(maxmin->max_lon-hrap.x)/maxmin->delta_lon;
589 
590  if(latindex >= 0 && lonindex >= 0 &&
591  latindex < maxmin->maxj-1 && lonindex < maxmin->maxi-1) {
592 
593  maxdistance=0.0;
594  mindistance=0.0;
595  maxvalue=0.0;
596  minvalue=0.0;
597 
598  for(kk=latindex;kk<=latindex+1;kk++) {
599 
600  maxminlat=maxmin->max_lat - (float)kk*maxmin->delta_lat;
601 
602  for(jj=lonindex;jj<=lonindex+1;jj++) {
603 
604  maxminlon=maxmin->max_lon - (float)jj*maxmin->delta_lon;
605 
606  dist1=maxminlat-irap.y;
607  dist2=(maxminlon-irap.x)*cos((maxminlat+irap.y)/2*conv);
608 
609  dist=pow(dist1,2)+pow(dist2,2);
610 
611  if(maxmin->maxvalue[h][jj][kk] >-500) {
612 
613  maxvalue=maxvalue+dist*maxmin->maxvalue[h][jj][kk];
614  maxdistance=maxdistance+dist;
615 
616  }
617 
618 
619  if(maxmin->minvalue[h][jj][kk] >-500) {
620 
621  minvalue=minvalue+dist*maxmin->minvalue[h][jj][kk];
622  mindistance=mindistance+dist;
623 
624  }
625 
626  }
627  }
628 
629  if(maxdistance != 0.0)
630  hrap_grid->max[h][i][k]=maxvalue/maxdistance;
631 
632  if(mindistance != 0.0)
633  hrap_grid->min[h][i][k]=minvalue/mindistance;
634 
635  }
636 
637  }
638 
639  }
640 
641  }
642 
643  }
644 
645 fclose(fz);
646 return;
647 
648 }
649 
650 
651 
652 
653 
654 
655 
656 
657 
658 
659 
660 
max_stations
int max_stations
Definition: daily_qc.c:199
s
char s
Definition: build_list.c:122
tpf
struct pcp * tpf
Definition: daily_qc.c:279
pcp_in_use
int pcp_in_use[500]
Definition: daily_qc.c:243
hrap_gage_file
char hrap_gage_file[1000]
Definition: daily_qc.c:186
top_level
Widget top_level
Definition: daily_qc.c:214
pcp
struct pcp * pcp
Definition: daily_qc.c:277
hrap_grid
struct hrap_grid * hrap_grid
Definition: daily_qc.c:285
scrollbar
Widget scrollbar
Definition: daily_qc.c:215
isohyets_used
int isohyets_used
Definition: daily_qc.c:192
posit
int posit
Definition: daily_qc.c:255
emonth
int emonth
Definition: daily_qc.c:142
spf
struct pcp * spf
Definition: daily_qc.c:278
station_list_file
char station_list_file[1000]
Definition: daily_qc.c:173
topo
struct topo * topo
Definition: daily_qc.c:286
isoh
struct isoh * isoh
Definition: daily_qc.c:287
s_text
Widget s_text
Definition: daily_qc.c:214
maxmin
struct maxmin * maxmin
Definition: daily_qc.c:22
maxmin_used
int maxmin_used
Definition: daily_qc.c:17
smonth
int smonth
Definition: daily_qc.c:142
maxj
int maxj
Definition: get_frzlvl_data.c:3
maxi
int maxi
Definition: get_frzlvl_data.c:3
get_hrap_coord
void get_hrap_coord(double maximum_latitude, double minimum_latitude, double center_longitude, int smonth, int emonth)
Definition: get_hrap_coord.c:4
printf
printf("pcp %d\n", pcp_in_use[103])
exit
exit(1)
LatLongToHrap
struct HRAP LatLongToHrap(float, float)
HrapToLatLong
struct HRAP HrapToLatLong(struct HRAP)
prototypes_new.h
HRAP
Definition: misc.h:407
HRAP::y
float y
Definition: misc.h:410
HRAP::x
float x
Definition: misc.h:409
coord
Definition: misc.h:422
coord::lat
float lat
Definition: misc_new.h:673
coord::y
short int y
Definition: misc.h:425
coord::lon
float lon
Definition: misc_new.h:674
coord::x
short int x
Definition: misc.h:424
dval
Definition: misc.h:397
dval::lo
double lo
Definition: misc.h:402
dval::a
double a
Definition: misc.h:399
dval::yo
double yo
Definition: misc.h:401
dval::xo
double xo
Definition: misc.h:400
gage
Definition: misc.h:517
gage::index
short int index[20]
Definition: misc.h:519
hrap_grid
Definition: misc.h:505
hrap_grid::hrap_miny
short int hrap_miny
Definition: misc.h:510
hrap_grid::min
short int *** min
Definition: misc_new.h:814
hrap_grid::hrap_minx
short int hrap_minx
Definition: misc.h:509
hrap_grid::maxi
short int maxi
Definition: misc.h:507
hrap_grid::maxj
short int maxj
Definition: misc.h:508
hrap_grid::coord
struct coord ** coord
Definition: misc.h:511
hrap_grid::gage
struct gage ** gage
Definition: misc.h:512
hrap_grid::isoh
short int ** isoh
Definition: misc.h:513
hrap_grid::elev
short int ** elev
Definition: misc_new.h:815
hrap_grid::max
short int *** max
Definition: misc_new.h:813
hrap_grid::owner
short int ** owner
Definition: misc_new.h:811
isoh
Definition: misc.h:465
isoh::delta_lon
float delta_lon
Definition: misc.h:476
isoh::max_lon
float max_lon
Definition: misc.h:472
isoh::max_lat
float max_lat
Definition: misc.h:471
isoh::value
short int ** value
Definition: misc.h:468
isoh::delta_lat
float delta_lat
Definition: misc.h:475
maxmin::max_lat
float max_lat
Definition: misc_new.h:769
maxmin::delta_lat
float delta_lat
Definition: misc_new.h:773
maxmin::delta_lon
float delta_lon
Definition: misc_new.h:774
maxmin::max_lon
float max_lon
Definition: misc_new.h:770
maxmin::maxvalue
short int *** maxvalue
Definition: misc_new.h:765
maxmin::minvalue
short int *** minvalue
Definition: misc_new.h:766
pcp
Definition: misc.h:523
pcp::value
short int ** value
Definition: misc.h:526
station
Definition: misc.h:232
station::lat
float lat
Definition: misc.h:239
station::lon
float lon
Definition: misc.h:240
topo
Definition: misc.h:450
topo::delta_lat
float delta_lat
Definition: misc.h:460
topo::max_lat
float max_lat
Definition: misc.h:456
topo::maxi
int maxi
Definition: misc.h:454
topo::delta_lon
float delta_lon
Definition: misc.h:461
topo::value
short int ** value
Definition: misc.h:453
topo::max_lon
float max_lon
Definition: misc.h:457
topo::maxj
int maxj
Definition: misc.h:455
Generated by doxygen 1.9.1