/* (c) copyright 2003 David G. Long, Brigham Young University */ /***************************************************************** Filename: slice.c This brief programs illustrates how to read a QuikSCAT L1B file and generate 3 dB contour slice shapes using BYU-generated slice shape files. written by DGL at BYU 6/23/2003 revised by DGL at BYU 2/02/2007 + fixed slice longitude scaling note: some sample code to read HDF files is included, though the main purpose of the code is to illustrate how to use the shape tables. This code should be compiled to include the HDF header directory ******************************************************************/ #include #include #include #include /* HDF include (local installation) */ #define abs(a) (((a) > 0.0) ? (a) : -(a)) #define slice_start 0 /* slices to look at */ #define slice_stop 8 /* define enum and structure type to store hdf variable arrays efficently */ typedef enum {flt=0, uint4=1, uint2=2, uint1=3} storesize; typedef struct{ int rank,dim[3]; storesize size; union { float *flt; unsigned int *uint4; unsigned short int *uint2; unsigned char *uint1; } data; } qsdata; void free_qsdata(qsdata p) /* free allocated structure memory */ { free(p.data.flt); p.data.flt=NULL; } /* function prototypes for HDF reading routines */ void *get_mem(int32 dt,int32 numval); void getTextAfterNewline(); qsdata extract_sds(int32 sd_id, char *in_var, int irec, int slab_size, storesize size); void read_attrib_byname(int32 sd_id, char *name, int *ntype, int *nval, char *val); void set_masks(unsigned long *frame_mask, unsigned short *frame_qual_mask, unsigned long *inst_mask, unsigned short *sigma_mask, unsigned short *qual_mask, unsigned long *slice_mask); /* slice shape routine function prototypes */ int load_inst_mode_shape(int inst_mode); void get_shape(float time, float az, int ib, int islice, float cen_lat, float cen_lon, int *nc, float *cor_lat, float *cor_lon, float orb_period); /* global variables */ unsigned short frame_qual_mask; /* frame_qual_flag */ unsigned short sigma_mask; /* sigma0_mode_flag */ unsigned short qual_mask; /* sigma0_qual_flag */ unsigned long frame_mask; /* frame_err_status */ unsigned long inst_mask; /* frame_inst_status */ unsigned long slice_mask; /* slice_qual_flag */ int old_inst_mode = -1; /* instrument gate width mode, initially flagged as unset */ /* begin main program */ int main(int argc,char *argv[]) { char *filename, fn[500], *outpath, op[300], *c, *c2, outfile[500]; char line[1025], attr_name[512], attr_data2[512], attr_chk[512]; char values[10240], string[1024]; int32 retn, attr_index, adata_type, count, sd_id; int8 *attr_data; int nframes; qsdata orbit_time, sigma0_mode_flag, slice_qual_flag; qsdata frame_inst_status, frame_err_status, frame_qual_flag, num_pulses; qsdata cell_lat, cell_lon, antenna_azimuth, z_vel; qsdata slice_lat, slice_lon, slice_sigma0; qsdata slice_azimuth, slice_incidence; int ir, nframe, rev_num; int i, errflag; int ntype, nval; int iirev, iyears, idays; float rev_orbit_period; double orb_period; FILE *fid; int iframe, npulses, ii, j, ib, ibb, nc, ascend; float cen_lat, cen_lon, ant_az; float time_orb, dcen_lat, dcen_lon, dcor_lat, dcor_lon; int inst_mode; float nlat[8], nlon[8], cx, cy, a[8], b[8], w[32], amin, amax, sigma; int isgn, ix, iy, icnt=0; if (argc < 2) { fprintf(stdout,"\nusage: %s L1B_input_name \n",argv[0]); fprintf(stdout," input parameters:\n"); fprintf(stdout," L1B_input_name = input L1B hdf file (full path)\n"); return(0); } /* Read the input filename */ if (argv[1] == NULL){ printf("Enter name of input L1B file:"); fgets(line,sizeof(line),stdin); sscanf(line,"%s",fn); filename=fn; } else filename=argv[1]; printf("\nL1B Filename: %s\n",filename); /* set global data acceptance flags */ set_masks(&frame_mask, &frame_qual_mask, &inst_mask, &sigma_mask, &qual_mask, &slice_mask); /* Open the HDF input file and initiate the SD interface */ sd_id = SDstart(filename, DFACC_RDONLY); if (sd_id < 0) { printf("*** SDstart error \n"); exit(-1); } /* Make sure that the file is a QuikSCAT or SeaWinds Level 2B file */ attr_index=SDfindattr(sd_id,"ShortName"); retn=SDattrinfo(sd_id,attr_index,attr_name, &adata_type,&count); if (retn < 0) printf("*** error reading ShortName SDattrinfo\n"); attr_data=(int8 *)malloc(count * DFKNTsize(adata_type)); retn=SDreadattr(sd_id,attr_index,attr_data); if (retn < 0) printf("*** error reading ShortName SDreadattr\n"); getTextAfterNewline(attr_data,2,attr_data2); free(attr_data); printf("HDF ShortName: %s\n",attr_data2); if (strcmp(attr_data2,"QSCATL1B") != 0 && strcmp(attr_data2,"SWSL1B") != 0){ printf("ERROR: Input file is not a QuikSCAT/SeaWinds Level 1B file\n"); return(0); } /* read selected vdata header info into strings and decode into variables as needed*/ read_attrib_byname(sd_id,"rev_number",&ntype,&nval,values); getTextAfterNewline(values,2,string); iirev=atoi(string); rev_num=iirev; printf("L1B rev_number: %d\n",iirev); read_attrib_byname(sd_id,"l1b_actual_frames",&ntype,&nval,values); getTextAfterNewline(values,2,string); nframes=atoi(string); read_attrib_byname(sd_id,"RangeBeginningDate",&ntype,&nval,values); getTextAfterNewline(values,2,string); printf("Begin date: %s\n",string); sscanf(string,"%4d-%3d",&iyears,&idays); read_attrib_byname(sd_id,"RangeEndingDate",&ntype,&nval,values); getTextAfterNewline(values,2,string); printf("End date: %s\n",string); read_attrib_byname(sd_id,"RangeBeginningTime",&ntype,&nval,values); getTextAfterNewline(values,2,string); printf("Begin time: %s\n",string); read_attrib_byname(sd_id,"RangeEndingTime",&ntype,&nval,values); getTextAfterNewline(values,2,string); printf("End time: %s\n",string); read_attrib_byname(sd_id,"rev_orbit_period",&ntype,&nval,values); getTextAfterNewline(values,2,string); printf("rev_orbit_period: %s\n",string); rev_orbit_period=atof(string); sscanf(string,"%lf",&orb_period); if (rev_orbit_period < 100.0) rev_orbit_period=6061.0; printf("rev_orbit_period: %s %f\n",string,rev_orbit_period); printf("\nL1B file: %s\n",filename); printf("Year: %4d JD: %3d Rev: %d L1B frames: %d\n",iyears,idays,iirev,nframes); /* Read selected SDSs of the L1B file from frame */ ir=0; /* start frame */ nframe=nframes; /* end frame to read*/ printf("\nBegin reading HDF SDSs Number of frames: %d\n",nframes); orbit_time=extract_sds(sd_id,"orbit_time",ir,nframe,uint4); frame_inst_status=extract_sds(sd_id,"frame_inst_status",ir,nframe,uint4); frame_err_status=extract_sds(sd_id,"frame_err_status",ir,nframe,uint4); frame_qual_flag=extract_sds(sd_id,"frame_qual_flag",ir,nframe,uint2); num_pulses=extract_sds(sd_id,"num_pulses",ir,nframe,uint1); z_vel=extract_sds(sd_id,"z_vel",ir,nframe,flt); cell_lat=extract_sds(sd_id,"cell_lat",ir,nframe,flt); cell_lon=extract_sds(sd_id,"cell_lon",ir,nframe,flt); sigma0_mode_flag=extract_sds(sd_id,"sigma0_mode_flag",ir,nframe,uint2); slice_qual_flag=extract_sds(sd_id,"slice_qual_flag",ir,nframe,uint4); antenna_azimuth=extract_sds(sd_id,"antenna_azimuth",ir,nframe,flt); slice_lat=extract_sds(sd_id,"slice_lat",ir,nframe,flt); slice_lon=extract_sds(sd_id,"slice_lon",ir,nframe,flt); slice_sigma0=extract_sds(sd_id,"slice_sigma0",ir,nframe,flt); slice_azimuth=extract_sds(sd_id,"slice_azimuth",ir,nframe,flt); slice_incidence=extract_sds(sd_id,"slice_incidence",ir,nframe,flt); printf("End reading SDSs\n"); /* close access to HDF file */ SDend(sd_id); /* set old instrument mode value to force change when first valid value available */ old_inst_mode=-3; /* open output file */ fid=fopen("slice.DAT","w"); /* for each frame */ for (iframe=ir; iframe < nframe; iframe++) { if ((iframe % 100) == 0) printf("Processing frame %d\n",iframe); /* only keep valid data */ if ((frame_err_status.data.uint4[iframe] & frame_mask) == 0 && (frame_inst_status.data.uint4[iframe] & inst_mask) == 0 && (frame_qual_flag.data.uint2[iframe] & frame_qual_mask) == 0) { /* get instrument gate width mode at start of frame and see if appropriate table already loaded */ inst_mode=((frame_inst_status.data.uint4[iframe] >> 4) % 8); if (inst_mode != old_inst_mode) { /* load new shape table */ printf(" Changing Instrument Gate Width Mode to %d\n",inst_mode); old_inst_mode=inst_mode; if (load_inst_mode_shape(inst_mode) < 0) { printf("*** could not load shape table for mode %d\n",inst_mode); exit(-1); } } if (z_vel.data.flt[iframe] > 0.0) ascend=1; else ascend=0; /* for each pulse in frame */ npulses=num_pulses.data.uint1[iframe]; for (i=0; i 90.0 && ant_az < 270.0) ibb++; /* aft-looking, otherwise forware-looking */ /* for demo, keep only measurements near equator */ if (abs(cen_lat) < 1.0) { /* for each slice */ for (ii=slice_start; ii nn ) break; } outputtext[ii] = '\0'; return; } /* utility routines */ void set_masks(unsigned long *frame_mask, unsigned short *frame_qual_mask, unsigned long *inst_mask, unsigned short *sigma_mask, unsigned short *qual_mask, unsigned long *slice_mask) { int i; /* set to ones all bits which we will reject if set */ /* frame_err_status */ *frame_mask = ~(1<<23 | 1<<24); /* frame_qual_flag */ *frame_qual_mask = ~(1 | 1<<1 | 1<<2 | 1<<3); /* frame_inst_status */ *inst_mask = ~(1<<2 | 1<<4 | 1<<5 | 1<<6 | 1<<7 | 1<<8 | 1<<9 | 1<<10 | 1<<11 | 1<<23); /* sigma0_mode_flag */ *sigma_mask = ~(1<<2 | 1<<3 | 1<<6 | 1<<7 | 1<<8 | 1<<10); /* sigma0_qual_flag */ *qual_mask = ~(1<<1 | 1<<2 | 1<<3); /* slice_qual_flag */ for (i=0; i<8; i++) *slice_mask = *slice_mask | (1<<(i*4) | 1<<(i*4+1) | 1<<(i*4+2)); *slice_mask = ~*slice_mask; printf("Masks %X %X %X %X %X %X\n",*frame_mask, *frame_qual_mask, *inst_mask, *sigma_mask, *qual_mask, *slice_mask); /* c frame_err_status c frame_mask=ibclr(frame_mask,12) ! clear SES Data Overrun Flag #?# c frame_mask=ibclr(frame_mask,13) ! clear SES Data Underrun Flag #?# c frame_mask=ibclr(frame_mask,14) ! clear Payload Bus Interface (PBI) Flag #?# frame_mask=ibclr(frame_mask,23) ! clear Attitude Data Flag #?# frame_mask=ibclr(frame_mask,24) ! clear Ephemeris Data Flag #?# c frame_qual_flag frame_qual_mask=ibclr(frame_qual_mask,0) ! clear frame filler bit 1 frame_qual_mask=ibclr(frame_qual_mask,1) ! clear frame filler bit 2 frame_qual_mask=ibclr(frame_qual_mask,2) ! clear frame CRC flag bit 1 frame_qual_mask=ibclr(frame_qual_mask,3) ! clear frame CRC flag bit 2 c frame_inst_status inst_mask=ibclr(inst_mask,2) ! clear first pulse count bit inst_mask=ibclr(inst_mask,4) ! clear gate width bits inst_mask=ibclr(inst_mask,5) ! clear gate width bits inst_mask=ibclr(inst_mask,6) ! clear gate width bits inst_mask=ibclr(inst_mask,7) ! clear hires/egg acquisition flag inst_mask=ibclr(inst_mask,8) ! clear begin calibration bit inst_mask=ibclr(inst_mask,9) ! clear SES select bit inst_mask=ibclr(inst_mask,10) ! clear SAS select bit inst_mask=ibclr(inst_mask,11) ! clear TWTA select bit inst_mask=ibclr(inst_mask,23) ! clear TWTA monitor enable bit #?# c sigma0_mode_flag sigma_mask=ibclr(sigma_mask,2) ! clear beam flag bit sigma_mask=ibclr(sigma_mask,3) ! clear location bit sigma_mask=ibclr(sigma_mask,6) ! clear gate width bits sigma_mask=ibclr(sigma_mask,7) ! clear gate width bits sigma_mask=ibclr(sigma_mask,8) ! clear gate width bits sigma_mask=ibclr(sigma_mask,10) ! clear SES select bit c sigma_mask=ibclr(sigma_mask,11) ! clear spin rate bit c sigma0_qual_flag c qual_mask=ibclr(qual_mask,0) ! clear sigma0 useable bit #?# qual_mask=ibclr(qual_mask,1) ! clear low snr bit #?# qual_mask=ibclr(qual_mask,2) ! clear negative sigma-0 bit qual_mask=ibclr(qual_mask,3) ! clear sigma0 out of range bit #?# c slice_qual_flag do i=1,8 j=(i-1)*4 slice_mask=ibclr(slice_mask,j+1) ! clear sign flag bit slice_mask=ibclr(slice_mask,j+0) ! clear gain threshold bit #?# slice_mask=ibclr(slice_mask,j+2) ! clear low SNR flag bit #?# c slice_mask=ibclr(slice_mask,j+3) ! clear location bit #?# end do */ } /*********************************************************************************/ /* slice shape code */ float rinterp(float xx, float yy, float r1, float r2, float r3, float r4) { return((r2-r1)*yy+r1+((r4-r3)*yy+r3-(r2-r1)*yy-r1)*xx); } /* shape table storage area */ static struct { int norbs,nazs; float orbtimes[36],az_angs[8*36*36*2],corners[16*8*36*36*2]; } shape; /* load slice shape table into memory */ int load_inst_mode_shape(int mode) { char *loc_path, loc_file[512], loc_file_string[]="."; FILE *fid; int i,ib,iaz,iorb,islice; /* get path to shape table locations from environmental variable BYU_SHAPES_DIR set to to path to local location prior to call */ if ((loc_path =(char *) getenv("BYU_SHAPES_DIR"))==NULL) { loc_path = loc_file_string; printf("*** The environment variable BYU_SHAPES_DIR is not set/\n"); } sprintf(loc_file,"%s/shape_in%d.tab",loc_path,mode+1); printf(" Reading shape table file %s\n",loc_file); fid=fopen(loc_file,"r"); if (fid==NULL) { printf("*** could not open shape table file %s\n",loc_file); return(-1); } ib=0; fscanf(fid,"%d %d",&shape.norbs,&shape.nazs); for (i=0;i= table_orbit_period) t=fmod(t,table_orbit_period); for (i=0; i shape.norbs) { j=j-shape.norbs; if (t >= shape.orbtimes[i] && t < shape.orbtimes[j]+table_orbit_period) { io1=i; io2=j; goto label_10; } } else { if (t >= shape.orbtimes[i] && t < shape.orbtimes[j]) { io1=i; io2=j; goto label_10; } } } printf("*** orbit time error in get_shape %f %f %f %d\n",time,t,orb_period,shape.norbs); exit(-1); /* error return */ label_10: delo=shape.orbtimes[io2]-shape.orbtimes[io1]; if (io2 < io1) delo=delo+orb_period; delo=(t-shape.orbtimes[io1])/delo; a=fmod(az+720.0,360.0)*0.1; ia1=(int) a; ia2=ia1+1; if (ia2 >= shape.nazs) ia2=ia2-shape.nazs; dela=a-ia1; /* bilinearly interpolate tables, first the azimuth angle */ ind1=islice+ia1*8+io1*8*shape.nazs+ib*8*shape.nazs*shape.norbs; ind2=islice+ia1*8+io2*8*shape.nazs+ib*8*shape.nazs*shape.norbs; ind3=islice+ia2*8+io1*8*shape.nazs+ib*8*shape.nazs*shape.norbs; ind4=islice+ia2*8+io2*8*shape.nazs+ib*8*shape.nazs*shape.norbs; cd1=cos(DTR*shape.az_angs[ind1]); cd2=cos(DTR*shape.az_angs[ind2]); cd3=cos(DTR*shape.az_angs[ind3]); cd4=cos(DTR*shape.az_angs[ind4]); sd1=sin(DTR*shape.az_angs[ind1]); sd2=sin(DTR*shape.az_angs[ind2]); sd3=sin(DTR*shape.az_angs[ind3]); sd4=sin(DTR*shape.az_angs[ind4]); cx=rinterp(dela,delo,cd1,cd2,cd3,cd4); sx=rinterp(dela,delo,sd1,sd2,sd3,sd4); /* fix azimuth angle interpolation */ cd1=sqrt(cx*cx+sx*sx); cx=cx/cd1; sx=sx/cd1; /* for each slice corner, interpolate and rotate */ for (i=0; i<8; i++) { ind1=i+islice*16+ia1*8*16+io1*8*16*shape.nazs+ib*8*16*shape.nazs*shape.norbs; ind2=i+islice*16+ia1*8*16+io2*8*16*shape.nazs+ib*8*16*shape.nazs*shape.norbs; ind3=i+islice*16+ia2*8*16+io1*8*16*shape.nazs+ib*8*16*shape.nazs*shape.norbs; ind4=i+islice*16+ia2*8*16+io2*8*16*shape.nazs+ib*8*16*shape.nazs*shape.norbs; xx=rinterp(dela,delo,shape.corners[ind1], shape.corners[ind2], shape.corners[ind3], shape.corners[ind4]); yy=rinterp(dela,delo,shape.corners[ind1+8], shape.corners[ind2+8], shape.corners[ind3+8], shape.corners[ind4+8]); x[i] = cx*xx + sx*yy; y[i] = -sx*xx + cx*yy; } /* convert x,y corner to lat,lon */ xx=sin(DTR*cen_lat); r=(1.0-FLAT*xx*xx)*AEARTH; rsq=r*cos(DTR*cen_lat); for (i=0; i<8; i++) { temp=x[i]/rsq; if (temp > 1.0) temp = 1.0; if (temp < -1.0) temp = -1.0; temp=RTD*asin(temp); cor_lat[i]=RTD*asin((y[i]-(1.0-cos(DTR*temp))*sin(DTR*cen_lat)*rsq)/r)+cen_lat; cor_lon[i]=temp+cen_lon; } *nc=8; return; }