/*************************************************************************** * * * ----------- * * ssmitdrta 1.2 * * ----------- * * * * Read NESDIS SSM/I TDR file, output antenna temperatures. Also * * output satellite ID #, time, revolution #, latitudes, and * * longitudes. * * * * To convert the antenna temperatures to brightness temperatures, * * run the Fortran program ssmitdrtb. It will read the Ta file * * produced by this program and output Tb's. * * * *************************************************************************** * * * Usage: ssmitdrta in_file out_file * * * * Two parameters are given on the command line; * * (1) input file name - the NESDIS TDR data set, and * * (2) output file name - the Ta output file. * * * *************************************************************************** * * * Written 11/94 by Don Moss, University of Alabama in Huntsville * * * * Revised 3/95 by Don Moss for DPR 95-092 - Release 1.1 * * Define fseek parameters for non-standard compilers, i.e. not * * compliant with ISO 9899:1990 sect. 7.9.9.2. * * Report statistics at end. * * * * Revised 5/97 by Don Moss for DOR 97-239 - Release 1.2 * * Enable reading of new NESDIS "NPR" format. * * * *************************************************************************** * * * Output File: * * * * Satellite ID #, time, revolution #, latitudes, longitudes, * * surface type codes, and Ta's are output to a formatted * * (text) file. There are 16 records for each A-B scan pair. * * The first record contains satellite ID number (for example, * * 13 corresponds to DMSP F13), time (B-scan start time in * * seconds since the beginning of 1987), and approximate * * revolution count since launch. The second record contains * * 128 A-scan latitudes. The third record contains 128 B-scan * * latitudes. Then follow two records with the corresponding * * longitudes. (Note: South latitudes are represented as * * negative numbers. Longitudes are degrees east of the prime * * meridian; thus east longitudes are between 0 and 180, and * * west longitudes are between 180 and 359.99.) The next two * * records contain 128 A-scan and 128 B-scan surface types, * * respectively. Next are five records containing Ta's for * * the low-resolution channels, 19V, 19H, 22V, 37V, and 37H, * * in that order. The last four records of the set contain * * the high-resolution channel data: A-scan 85V, B-scan 85V, * * A-scan 85H, and B-scan 85H. These 16 records are output * * for each input record. The records are formatted as * * follows: * * * * Rec Format Length Contents * * --------------------------------------------- * * 1 I3,F16.4,I9 28 Sat ID, time, rev * * 2 128 F7.2 896 A-scan latitudes * * 3 128 F7.2 896 B-scan latitudes * * 4 128 F7.2 896 A-scan longitudes * * 5 128 F7.2 896 B-scan longitudes * * 6 128 I2 256 A-scan surface types * * 7 128 I2 256 B-scan surface types * * 8 64 F7.2 448 19V * * 9 64 F7.2 448 19H * * 10 64 F7.2 448 22V * * 11 64 F7.2 448 37V * * 12 64 F7.2 448 37H * * 13 128 F7.2 896 85V A-scan * * 14 128 F7.2 896 85V B-scan * * 15 128 F7.2 896 85H A-scan * * 16 128 F7.2 896 85H B-scan * * --------------------------------------------- * * * ***************************************************************************/ #include #include #include #include #define TRUE 1 #define FALSE 0 #ifndef SEEK_SET #define SEEK_SET 0 #define SEEK_CUR 1 #define SEEK_END 2 #endif #define ERROUT stdout #define BUFSIZE 4444 #define SPN_FRAME 12798 #define NPR_FRAME 3604 #define NPR_HDRS 2158 #define SPN_FILLER 0xA5 #define EOF_LEN 6 #define LORES 64 #define HIRES 128 /*-----return codes-----*/ #define BAD_LEN -1 #define PAST_EOF -2 #define OK 0 #define END_FILE 1 #define CANT_OPEN 2 #define BAD_HDRS 3 #define BAD_EOF 4 #define WRITE_ERR 6 #define FATAL_ERR 8 /*-----TDR block names-----*/ #define ProdID 0 #define DataSeq 1 #define RevHdrDD 2 #define ScanHdr1DD 3 #define ScanHdr2DD 4 #define ScanDD 5 #define RevHdr 6 #define ScanHdr1 7 #define ScanHdr2 8 #define Scan 9 /*-----macros-----*/ /* * V1 evaluates a 1-byte binary value as an unsigned number * V2 evaluates a 2-byte binary value as an unsigned number * V4 evaluates a 4-byte binary value as an unsigned number * VLat gets latitude and converts its range to [-90:+90] * VLon gets longitude * FPRTVAL writes a formatted value to the output file * FPRTNL writes a newline character to the output file */ #define V1(p) ((long) *(buf+p)) #define V2(p) (256L * (long) *(buf+p) + (long) *(buf+p+1)) #define V4(p) (256L * (256L * (256L * (long) *(buf+p) + (long) *(buf+p+1)) + (long) *(buf+p+2)) + (long) *(buf+p+3)) #define VLat(p) (V2(p) - 9000L) #define VLon(p) (V2(p)) #define FPRTVAL(fmt, val) if (fprintf(outfp, fmt, val) < 0) return TRUE #define FPRTNL if (fprintf(outfp, "\n") < 0) return TRUE /*-----program name and version-----*/ #define PROG_NAME "ssmitdrta" /* name of this program */ /***************************************************/ #define RELEASE_LVL "1.2" /* NOTE: change release level with every revision */ /***************************************************/ /*-----types-----*/ typedef char bool; #if INT_MAX < 0x7FFFFFFF typedef long int32; #else typedef int int32; #endif #if USHRT_MAX < 0xFFFF typedef unsigned int uint16; typedef int int16; #else typedef unsigned short uint16; typedef short int16; #endif /*-----functions-----*/ int32 GetBlock(void); /* read a block from TDR file */ int32 Checksum(int32 len); /* calculate block checksum (should be 0) */ int32 ReadValid(int32 blktype); /* read and validate a header block */ void FillArrays(void); /* fill scan arrays */ bool WriteValues(void); /* write values to output file */ void ShowStat(int nominal_in, int processed); /* report statistics */ /*-----header info-----*/ int32 cyr, cmon, cday, /* file creation date */ chr, cmin, /* file creation time */ scans, /* number of scans in file (from DataSeq) */ scid, /* spacecraft ID */ rev, /* nominal rev # */ bjld, bhr, bmin, bsec, /* begin day of year, time */ ejld, ehr, emin, esec, /* ending day of year, time */ ajld, ahr, amin, asec, /* ascending node doy, time */ lsat; /* logical satellite ID */ /*-----scan data-----*/ int32 scann; /* scan number (from ScanHdr1) */ int32 bst; /* B-scan start time (sec) scaled by 10000 */ double xtime; /* bst as seconds since 0z 1 Jan 1987 */ uint16 v19[LORES], /* Ta's */ h19[LORES], v22[LORES], v37[LORES], h37[LORES], v85[HIRES][2], h85[HIRES][2], lon[HIRES][2]; /* longitudes */ int16 lat[HIRES][2]; /* latitudes */ char sft[HIRES][2]; /* surface types */ /*-----other global variables-----*/ FILE *infp; /* input file pointer */ FILE *outfp; /* output file pointer */ unsigned char buf[BUFSIZE]; /* input file buffer */ int32 blklen[] = {28, 32, 190, 370, 1138, 370, 30, 76, 194, 3334}; char *blknam[10] = {"ProdID", "DataSeq", "RevHdrDD", "ScanHdr1DD", "ScanHdr2DD", "ScanDD", "RevHdr", "ScanHdr1", "ScanHdr2", "Scan"}; char pgm[] = PROG_NAME; /* this program's name */ int32 spotL, spotH; /* low- & high-res spot index */ bool SPN; /* true if reading a NESDIS SPN file */ bool NPR; /* true if reading a NESDIS NPR file */ int framesize; /* frame size of SPN or NPR file */ bool verbose = TRUE; /* verbose mode */ int scans_read = 0; /* number of scans successfully read */ /*================================ main ================================*/ int main(int argc, char *argv[]) { char *infn; /* input file name pointer */ char *outfn; /* output file name pointer */ char check[8] = {0, 14, 1, 1, 'F', 'N', 'O', 'C'}; /* the data above is what's supposed to be in the Product ID Block */ int32 len; /* block length */ int32 len0301 = 2 * 0x0301; /* quasi-length of certain erroneous SPN blocks */ int32 len0681 = 2 * 0x0681; /* quasi-length of certain erroneous SPN blocks */ int32 offset; /* offset to next frame, for skipping bad blocks */ int32 status; /* error condition */ int32 fcyr, fcmon, fcday; /* TDR file creation date */ int32 fchr, fcmin; /* and time */ int32 year; /* year of data */ int32 iyr; /* intermediate variable */ int32 ttime; /* offset from beg. 1987 till beg. current day */ int32 satid; /* satellite ID */ int32 bst0 = 0; /* previous bst */ /*-------------------------------------------------------------------------*/ /*-----check arguments-----*/ if (argc != 3) { fprintf(ERROUT, "\n%s %s\n", pgm, RELEASE_LVL); fprintf(ERROUT, "Usage: %s in_file out_file\n", pgm); return FATAL_ERR; } /*-----open input file-----*/ infn = argv[1]; if (!(infp = fopen(infn, "rb"))) { fprintf(ERROUT, "%s: can't open input file %s\n", pgm, infn); return CANT_OPEN; } else { if (verbose) printf("-----opening input file %s-----\n", infn); } /*-----open output file-----*/ outfn = argv[2]; if (!(outfp = fopen(outfn, "w"))) { fprintf(ERROUT, "%s: can't open output file %s\n", pgm, outfn); return CANT_OPEN; } else { if (verbose) printf("-----opening output file %s-----\n", outfn); } /*-----read first NPR frame-----*/ fseek(infp, NPR_HDRS, SEEK_SET); /* skip NPR headers */ if (fread(buf, 2, 1, infp) != 1) { /* read next two bytes */ fprintf(ERROUT, "\n%s: ***** Error reading first frame *****\n", pgm); return BAD_HDRS; } /*-----check for NPR filler-----*/ NPR = (buf[0] == 0 && buf[1] == 0); if (NPR) { SPN = FALSE; framesize = NPR_FRAME; } else { /*-----check for SPN filler-----*/ fseek(infp, SPN_FRAME - 2, SEEK_SET); /* skip to end of first SPN frame */ if (fread(buf, 2, 1, infp) != 1) { /* read last two bytes of it */ fprintf(ERROUT, "\n%s: ***** Error reading first frame *****\n", pgm); return BAD_HDRS; } SPN = (buf[0] == SPN_FILLER && buf[1] == SPN_FILLER); framesize = SPN_FRAME; /* make it nonzero in any case */ } /*-----reset file-----*/ fseek(infp, 0, SEEK_SET); /* return to beginning of file */ /*-------------------------------------------------------------------------*/ /*-----process non-repeating header blocks-----*/ /*---Product ID Block---*/ if (ReadValid(ProdID )) return BAD_HDRS; if (memcmp(buf, check, 8)) { fprintf(ERROUT, "%s: invalid ProdID block\n", pgm); return BAD_HDRS; } satid = 10 * (V1(18) - 48) + V1(19) - 48; fcyr = V2(20); fcmon = V1(22); fcday = V1(23); fchr = V1(24); fcmin = V1(25); if (verbose) printf("SatID %d File date %d/%02d/%02d time %02d:%02d\n", satid, fcyr, fcmon, fcday, fchr, fcmin); /*---Data Sequence Block---*/ if (ReadValid(DataSeq )) return BAD_HDRS; scans = V2(14); if (verbose) printf("Scans %d\n", scans); /*---Data Description Blocks---*/ if (ReadValid(RevHdrDD )) return BAD_HDRS; if (ReadValid(ScanHdr1DD)) return BAD_HDRS; if (ReadValid(ScanHdr2DD)) return BAD_HDRS; if (ReadValid(ScanDD )) return BAD_HDRS; /*---Rev Header Block---*/ if (ReadValid(RevHdr )) return BAD_HDRS; scid = V4(4); rev = V4(8); bjld = V2(12); bhr = V1(14); bmin = V1(15); bsec = V1(16); ejld = V2(17); ehr = V1(19); emin = V1(20); esec = V1(21); ajld = V2(22); ahr = V1(24); amin = V1(25); asec = V1(26); lsat = V1(27); if (verbose) { printf("SCID %d REV# %d\n", scid, rev); printf("BJLD %03d %02d:%02d:%02d\n", bjld, bhr, bmin, bsec); printf("EJLD %03d %02d:%02d:%02d\n", ejld, ehr, emin, esec); printf("AJLD %03d %02d:%02d:%02d LSAT %d\n", ajld, ahr, amin, asec, lsat); } /*---compute ttime (seconds from beg. 1987 till beg. current day---*/ year = fcyr; ttime = 0; /* initialize offset */ if (bjld > 330 && fcmon <= 2) year--; /* adjust year if necessary */ for (iyr = 1987; iyr < year; iyr++) /* sum days in */ ttime += (iyr % 4) ? 365 : 366; /* previous years */ ttime += (bjld - 1); /* add days through yesterday */ ttime *= 86400; /* convert to seconds */ /*-----if NPR, position to end of block-----*/ if (NPR) fseek(infp, framesize, SEEK_SET); /*-------------------------------------------------------------------------*/ /*-----read and process scans-----*/ for (;;) { /*---read Scan Header Block 1---*/ if ((len = GetBlock()) == BAD_LEN) return FATAL_ERR; else if (len == EOF_LEN || len == 0 || len == PAST_EOF) { status = END_FILE; break; /* exit for */ } else if (len != blklen[ScanHdr1]) { printf("%s: unexpected block length %d, now @ %08X\n", pgm, len, ftell(infp)); if (SPN && len == len0301) { printf(" -- skipping to next frame\n"); offset = framesize - (ftell(infp) % framesize); fseek(infp, offset, SEEK_CUR); } continue; } else if (Checksum(len)) { printf("\n%s: invalid hdr1 block, now @ %08X\n", pgm, ftell(infp)); continue; } /*---extract parameters from Scan Header Block 1---*/ scann = V2(4); /* not used */ bst = V4(6); /*---convert time to seconds since beginning of 1987---*/ if (bst < bst0) /* if time reverses (crosses midnight), */ ttime += 86400; /* add a day's worth of seconds to offset */ bst0 = bst; /* keep up with previous bst */ xtime = ttime + (double) bst / 10000.0; /*---read Scan Header Block 2---*/ if ((len = GetBlock()) == BAD_LEN) return FATAL_ERR; else if (len == EOF_LEN || len == 0) { status = BAD_EOF; break; /* exit for */ } else if (len != blklen[ScanHdr2]) { printf("%s: unexpected block length %d, now @ %08X\n", pgm, len, ftell(infp)); continue; } else if (Checksum(len)) { printf("\n%s: invalid hdr2 block, now @ %d\n", pgm, ftell(infp)); continue; } /*---read Scan Data Block---*/ if ((len = GetBlock()) == BAD_LEN) return FATAL_ERR; else if (len == EOF_LEN || len == 0) { status = BAD_EOF; break; /* exit for */ } else if (len != blklen[Scan]) { printf("\n%s: unexpected block length %d, now @ %08X\n", pgm, len, ftell(infp)); if (SPN && len == len0681) { printf(" -- skipping to next frame\n"); offset = framesize - (ftell(infp) % framesize); fseek(infp, offset, SEEK_CUR); } continue; } else if (Checksum(len)) { printf("\n%s: invalid data block, cksm %04X, now @ %08X\n", pgm, Checksum(len), ftell(infp)); continue; } scans_read++; /*---extract parameters from Scan Block---*/ FillArrays(); /*-----write values-----*/ if (WriteValues()) { fprintf(ERROUT, "%s: error writing output file\n"); return WRITE_ERR; } } /*-----end of for loop to read & process scans-----*/ /*-------------------------------------------------------------------------*/ /*-----check status, end program-----*/ switch (status) { case END_FILE: fclose(outfp); ShowStat(scans, scans_read); return OK; case BAD_EOF: fprintf(ERROUT, "%s: unexpected end-of-file encountered\n", pgm); ShowStat(scans, scans_read); return BAD_EOF; } } /*================================ ReadValid ================================*/ int32 ReadValid(int32 blktype) /* * Read and validate a header block. */ { int32 length = blklen[blktype]; /*-----read-----*/ if (fread(buf, length, 1, infp) != 1) { fprintf(ERROUT, "%s: can't read %s block, now @ %d\n", pgm, blknam[blktype], ftell(infp)); return BAD_HDRS; } /*-----validate-----*/ if (Checksum(length)) { fprintf(ERROUT, "%s: invalid %s block, now @ %d\n", pgm, blknam[blktype], ftell(infp)); return BAD_HDRS; } /*-----return-----*/ return 0; } /*================================ GetBlock =================================*/ int32 GetBlock() /* * Get a block if possible, and return the block length. Otherwise, if initial * read is past EOF, return PAST_EOF. Otherwise, return BAD_LEN. */ { int32 len, /* block length in bytes */ len2, /* intermediate variable */ offset; /* offset to next frame, for skipping bad blocks */ bool bad, /* true if block length is bad */ skip; /* true if reading SPN filler */ /*-----loop until a good block is read-----*/ do { do { /*-----read block length (i.e. number of two-byte words)-----*/ if (fread(buf, 2, 1, infp) != 1) { if (feof(infp)) return PAST_EOF; else { fprintf(ERROUT, "%s: can't read next block length code\n", pgm); return BAD_LEN; } } /*-----if an SPN filler was read, skip to next block and read again-----*/ skip = (SPN && buf[0] == SPN_FILLER && buf[1] == SPN_FILLER); if (skip) { offset = framesize - (ftell(infp) % framesize); fseek(infp, offset, SEEK_CUR); } } while (skip); /*-----convert length to bytes-----*/ len = V2(0) * 2; /*-----bytes remaining to be read-----*/ len2 = len - 2; /*-----is it valid?-----*/ bad = (len != 0 && (len < 4 || len > BUFSIZE)); if (bad) { printf("%s: ***** invalid block length %d, now @ %08X *****\n", pgm, len, ftell(infp)); /*-----no, skip to next frame if this is a NESDIS SPN or NPR file-----*/ if (SPN || NPR) { printf(" -- skipping to next frame\n"); offset = framesize - (ftell(infp) % framesize); fseek(infp, offset, SEEK_CUR); continue; } else { printf(" -- aborting\n"); return BAD_LEN; } } /*-----read the rest of the block-----*/ if (len2 > 0) { if (fread(buf + 2, len2, 1, infp) != 1) { fprintf(ERROUT, "%s: can't read next %d bytes\n", pgm, len2); return BAD_LEN; } } } while (bad); /*-----return the block length-----*/ return len; } /*================================ Checksum =================================*/ int32 Checksum(int32 len) { int32 i; int32 cksm = 0; for (i = 0; i < len; i+=2) cksm += V2(i); return cksm & 0xFFFF; } /*=============================== FillArrays ================================*/ void FillArrays(void) /* * Fill scan arrays. */ { int32 bb; /* offset into buffer */ for (spotL = 0; spotL < LORES; spotL++) { bb = spotL * 52; spotH = spotL * 2; v19[spotL] = V2(bb + 10); h19[spotL] = V2(bb + 12); v22[spotL] = V2(bb + 14); v37[spotL] = V2(bb + 16); h37[spotL] = V2(bb + 18); lat[spotH][0] = VLat(bb + 6); lon[spotH][0] = VLon(bb + 8); v85[spotH][0] = V2(bb + 20); h85[spotH][0] = V2(bb + 22); sft[spotH][0] = V1(bb + 24); lat[spotH][1] = VLat(bb + 26); lon[spotH][1] = VLon(bb + 28); v85[spotH][1] = V2(bb + 30); h85[spotH][1] = V2(bb + 32); sft[spotH][1] = V1(bb + 34); lat[spotH + 1][0] = VLat(bb + 36); lon[spotH + 1][0] = VLon(bb + 38); v85[spotH + 1][0] = V2(bb + 40); h85[spotH + 1][0] = V2(bb + 42); sft[spotH + 1][0] = V1(bb + 44); lat[spotH + 1][1] = VLat(bb + 46); lon[spotH + 1][1] = VLon(bb + 48); v85[spotH + 1][1] = V2(bb + 50); h85[spotH + 1][1] = V2(bb + 52); sft[spotH + 1][1] = V1(bb + 54); } } /*=============================== WriteValues ===============================*/ bool WriteValues(void) /* * Write values to output file. Return TRUE in the event of a write error; * otherwise return FALSE. */ { /*-----satellite ID, B-scan start time, rev #-----*/ if (fprintf(outfp, "%3d%16.4f%9d\n", scid, xtime, rev) < 0) return TRUE; /*-----A-scan latitudes-----*/ for (spotH = 0; spotH < HIRES; spotH++) FPRTVAL("%7.2f", (float) lat[spotH][0] / 100.0); FPRTNL; /*-----B-scan latitudes-----*/ for (spotH = 0; spotH < HIRES; spotH++) FPRTVAL("%7.2f", (float) lat[spotH][1] / 100.0); FPRTNL; /*-----A-scan longitudes-----*/ for (spotH = 0; spotH < HIRES; spotH++) FPRTVAL("%7.2f", (float) lon[spotH][0] / 100.0); FPRTNL; /*-----B-scan longitudes-----*/ for (spotH = 0; spotH < HIRES; spotH++) FPRTVAL("%7.2f", (float) lon[spotH][1] / 100.0); FPRTNL; /*-----A-scan surface types-----*/ for (spotH = 0; spotH < HIRES; spotH++) FPRTVAL("%2i", sft[spotH][0]); FPRTNL; /*-----B-scan surface types-----*/ for (spotH = 0; spotH < HIRES; spotH++) FPRTVAL("%2i", sft[spotH][1]); FPRTNL; /*-----19V-----*/ for (spotL = 0; spotL < LORES; spotL++) FPRTVAL("%7.2f", (float) v19[spotL] / 100.0); FPRTNL; /*-----19H-----*/ for (spotL = 0; spotL < LORES; spotL++) FPRTVAL("%7.2f", (float) h19[spotL] / 100.0); FPRTNL; /*-----22V-----*/ for (spotL = 0; spotL < LORES; spotL++) FPRTVAL("%7.2f", (float) v22[spotL] / 100.0); FPRTNL; /*-----37V-----*/ for (spotL = 0; spotL < LORES; spotL++) FPRTVAL("%7.2f", (float) v37[spotL] / 100.0); FPRTNL; /*-----37H-----*/ for (spotL = 0; spotL < LORES; spotL++) FPRTVAL("%7.2f", (float) h37[spotL] / 100.0); FPRTNL; /*-----A-scan 85V-----*/ for (spotH = 0; spotH < HIRES; spotH++) FPRTVAL("%7.2f", (float) v85[spotH][0] / 100.0); FPRTNL; /*-----B-scan 85V-----*/ for (spotH = 0; spotH < HIRES; spotH++) FPRTVAL("%7.2f", (float) v85[spotH][1] / 100.0); FPRTNL; /*-----A-scan 85H-----*/ for (spotH = 0; spotH < HIRES; spotH++) FPRTVAL("%7.2f", (float) h85[spotH][0] / 100.0); FPRTNL; /*-----B-scan 85H-----*/ for (spotH = 0; spotH < HIRES; spotH++) FPRTVAL("%7.2f", (float) h85[spotH][1] / 100.0); FPRTNL; /*-----return-----*/ return FALSE; /* indicate no error */ } /*=============================== ShowStat ==============================*/ void ShowStat(int nominal_in, int processed) { printf("\n\n========================================\n"); printf(" Nominal input scans . . . %5d\n", nominal_in); printf(" Scans processed . . . . . %5d\n", processed); printf("========================================\n"); }