There is a library called libsi.a that contains a number of routines for reading and writing satellite image files. The purpose of this library is to permit programmers to use files of this format without worrying about the specifics of the format. Functions currently exist to open an existing file and return a header structure, read a scanline, write a header structure, and write a scan line. Whether or not latitude/longitude data is contained in the image file or in a separate file should remain transparent to the caller. The following is a discussion of the available routines.
Satellite Image Files
There are library routines available for reading and writing satellite image files, as well as operating on satellite image file headers. Before reading data from a file, it is necessary to open it. Similarly, before writing data to a file, it must be created. In actuality, the open and create functions read and write only the header. There are also functions available for creating and copying headers.
Reading a Satellite Image File
There are two operations involved with reading a satellite image file:
The function si_open() opens a satellite image file (and its latitude/longitude file if it is separate) and returns a pointer to a SatImageHdr structure. This pointer may be used simply as a handle for subsequent read operations on the file, without ever accessing the members of the structure, if desired. The C binding is:
C Binding:
SatImageHdr *si_open(file,rwmode) char *file; int rwmode;
where file is the name of the file to open and rwmode is the read/write mode (0 = read only, 1 = write only, 2 = read & write). Embedded within the SatImageHdr structure is the file descriptor for the open file, as well as a file descriptor for the latitude/longitude file. Thus subsequent calls to si_read_scan_line() can access the open file(s) through the header structure.
Fortran Binding:
si_open(file,rwmode,hdr) character file*(*) integer rwmode record /SatImageHdr/ hdr
where file is the name of the file to open and rwmode is the read/write mode (0 = read only, 1 = write only, 2 = read & write). hdr is the SatImageHeader structure to be filled in. This header is used as a handle to the file for all subsequent read operations. The file sif.h must be included in the fortran source file because it contains the definition of the SatImageHdr record.
si_read_scanline()
After opening a satellite image file, a program will normally want to read the data. The data is accessed one scan line at a time through a series of calls to si_read_scanline(), usually in a loop:
C Binding:
si_read_scanline(hdr, par, lat, lon, num_pts, time) SatImageHdr *hdr; float *par, *lat, *lon; int *num_pts; float *time;
In this routine, hdr is the pointer that was handed back by si_open(). par, lat and lon are pointers to arrays of floats. Since no memory is allocated by si_read_scanline(), it is the calling routine's responsibility to make certain these pointers point to adequate storage. When si_read_scanline() returns, *num_pts contains the number of samples in the scanline just read, and *time contains the time of the scanline, if known (see the discussion of the time_flag in Appendix A). If the calling routine is not interested in one or more of the returned variables par, lat, lon, num_pts, or time, it may pass in a NULL pointer for that variable. si_read_scanline() returns: -1 for a failed read, 0 for a successful read, or 1 for a successful read with a discontinuity in the data (i.e., there is a gap between the scan line that was read and the one that preceded it.)
Fortran Binding:
si_read_scanline(hdr, par, lat, lon, num_pts, time) record /SatImageHdr/ hdr real par(*), lat(*), lon(*) integer*4 num_pts real time
The explanation is the same as for the C binding, except, obviously, NULL cannot be passed instead of a valid variable address.
si_read_raw_scanline()
si_read_raw_scanline() reads a single raw scan line of a satellite image from a file previously opened by si_open( ).
C Binding:
#include /usr/local/include/si.h si_read_raw_scanline(hdr, par, lat, lon, num_pts, time) SatImageHdr *hdr; char *par; float *lat, *lon; int *num_pts; float *time;hdr is the pointer to the SatImageHdr structure returned by si_open(). The behavior of si_read_raw_scanline() is identical to si_read_scanline(), except si_read_scanline() always returns FLOAT data in par, but si_read_raw_scanline() returns the data in whatever form it exists in the data file. Allowable types are float, 8-bit values, 16-bit values, and 32-bit values. The 8 and 16 bit values may be actual pixel intensities, or indices into a lookup table, if the file has a lookup table. par, lat, and lon are pointers to previously allocated arrays to hold the data read in. par will receive the data values, while lat and lon will receive the location information. In other words, the pixel intensity at lat[i],lon[i] is par[i]. The number of points in the scan line is returned in num_pts. If each scanline is time-stamped, the time for the scan line is returned in time. If any of the returned information is not required, a NULL pointer may be passed. It is the responsibility of the programmer to account for the different possible data types. In most cases, it is wiser to use only si_read_scanline(). si_read_raw_scanline() returns -1 for error (or end of file), 0 for a successful read, and 1 for a successful read with a discontinuity at the current scan line (i.e., there is a gap between the previous line and this one).
si_read_data()
All of the data in a satellite image file with a constant number of samples per scan line can be read in a single call to si_read_data().
C Binding:
si_read_data(hdr,par,lat,lon) SatImageHdr *hdr; float *par,*lat,*lon;
where hdr is a pointer to the header structure of the file being read, and par, lat and lon are either NULL (in which case that information is not returned) or pointers to pre-allocated arrays of a size large enough to accommodate all the data in the file. par receives the data values, and lat and lon receive the corresponding latitudes and longitudes.
Writing Satellite Image Files
There are several routines available to write out a satellite image file. These routines are analogous to the functions for reading a file.
si_write_header()
To create a satellite image file, call si_write_header(). This routine accepts a filled in header structure and creates a file with the appropriate header data. In some cases, a latitude/longitude file referenced in the header will already exist. In other cases, you might want to create it. si_write_header() permits either case.
C Binding:
si_write_header(file,hdr,latlon_flag) char *file; SatImageHdr *hdr; int latlon_flag; /* we create latlon file if TRUE */
in this function, file is the name of the file to be created. hdr is a pointer to a filled in header structure, and latlon_flag is 0 if no latitude/longitude file should be created, and 1 otherwise. The routine si_create_header() can be called prior to invoking si_write_header(). It will create a header structure and fill it in. The pointer it returns can be handed directly to si_write_header(). See the discussion of si_create_header() below.
Fortran Binding:
si_write_header(file,hdr,latlon_flag) character file*(*) record /SatImageHdr/ hdr integer*4 latlon_flagsi_write_scanline()
si_write_scanline() is analogous to si_read_scanline(). C Binding:
si_write_scanline(hdr, par, lat, lon, num_pts, time) SatImageHdr *hdr; float *par, *lat, *lon; int num_pts; float time;
hdr is a pointer to the header passed to si_write_scanline() which filled in the file descriptors elements of the header. par,lat and lon are pointers to arrays of dimension num_pts. par is the satellite image data. The lat and lon arrays must be included if:
si_write_scanline (hdr, par, lat, lon, num_pts, time) record /SatImageHdr/ hdr real*4 par(*), lat(*), lon(*) integer*4 num_pts real*4 timesi_write_raw_scanline
This routine exists, but it is a synonym for si_write_scanline, i.e., they do the same thing and are called exactly the same way.
Closing Satellite Data Files
When you are done reading or writing a file, close it with a call to si_close().
C Binding:
si_close(hdr) SatImageHdr *hdr;
hdr is a pointer to the header structure used to open the file. All open file descriptors associated with that header structure are closed.
Fortran Binding:
si_close(hdr) record /SatImageHdr/ hdrUtility Routines
Several utility routines exist for operating on SatImageHdr structures and satellite image files.
Copying Headers
C Binding:
SatImageHdr *si_dup_header(hdr,copy) SatImageHdr *hdr,*copy;
si_dup_header() copies a header structure and hands back a pointer to the new structure. This function is handy if you have several satellite image files to create with similar headers. Create the first one with si_create_header(), and subsequent ones with si_dup_header(). If copy is NULL, space for a new header is allocated. Otherwise, the contents of hdr are copied. Neither of these routines opens or creates a file. Only si_write_header() does that.
Fortran Binding:
si_dup_header_(src_hdr,dest_hdr) record /SatImageHdr/ src_hdr, dest_hdrCreating Headers
si_create_header() exists to permit the programmer to create and fill in a satellite image header structure with a single function call. A pointer to a SatImageHdr structure is returned. This pointer can then be passed to si_write_header().
C Binding:
SatImageHdr * si_create_header(sat_id, param, min, max, bad, ll_file, scans, samps, comment, year, month, day, time) int sat_id, param, scans, samps, year, month, day; char *ll_file, *comment; float min, max, bad, time;Fortran Binding:
si_create_header(hdr, sat_id, param, min, max, bad, ll_file, scans, samps, comment,year, month, day, time) integer*4 sat_id, param, scans integer*4 samps, year, month, day character ll_file*(*), comment*(*) real*4 min, max, bad, time
Called from Fortran, si_create_header() fills in an existing SatImageHdr structure.
Sometimes it is necessary to change a file header that has already been written. If you change an element of a satellite image header structure that has already been written out to a file, you may rewrite the updated header with a call to si_update_header().
C Binding:
si_update_header(hdr) SatImageHdr *hdr;Fortran Binding:
si_update_header(hdr) record /SatImageHdr/ hdrSkipping Scan Lines
When reading a satellite image file, some scan lines may not be of any interest. si_skip_scan_lines() permits the program to jump over an arbitrary number of scan lines in a file without reading them.
C Binding:
si_skip_scan_lines(hdr,num) SatImageHdr *hdr; int num;Fortran Binding:
si_skip_scan_lines(hdr,num) record /SatImageHdr/ hdr integer*4 num
hdr is a pointer to a satellite image file opened for reading. num is the number of scan lines to skip.si_skip_scanline() returns: -1 for failure, 0 for success, or 1 for success with a discontinuity in the data (i.e., there is a gap between the scan line that was read and the one that preceded it.)
Size of Header
To determine the current size of a header structure, use si_sizeof_header()
C Binding:
int si_sizeof_header(hdr) SatImageHdr *hdr;Fortran Binding:
integer function si_sizeof_header(hdr) record /SatImageHdr/ hdrTo Set Comment String (Fortran)
To set the value of the comment string header element from Fortran, use si_set_comment():
Fortran Binding:
si_set_comment(hdr, comment) record /SatImageHdr/ hdr character comment*(*)Compute Boundary of Satellite Image File
Call si_compute_boundary() (from C only) to return the vertices of a polygon that encloses the data (coverage) of a SI file.
C Binding:
si_compute_boundary(file,x,y,res) char *file; float **x; float **y; int res;
This routine computes the boundary of a single, continuous image. Arrays x,y contain coordinates of vertices comprising the boundary. The number of vertices is returned. res is the distance between adjacent input points/scans used to compute the boundary.
Overlay Files
There are library routines available for reading and writing overlay files, as well as operating on overlay data-set headers. Before reading data from a file, it is necessary to open it. Similarly, before writing data to a file, it must be created. The functions available for reading and writing overlay files are analogous to the Unix file system commands for manipulating files. In general, a file should be opened before it is read, and it should be created before it is written.
Reading an Overlay File
There are two operations involved with reading an overlay file:
The function ov_open() opens a file and checks that it is an overlay data file. It then returns a file descriptor for the open file. Subsequent calls to ov_read() will pass the file descriptor. The calling program may use the file descriptor to do normal Unix file system I/O operations, but there is really no good reason to do so.
C Binding:
ov_open(file,mode) char *file; int mode;
file is the name of the file to be opened. mode is 0 for read-only, 1 for write-only, and 2 for read-write permission.
Fortran Binding:
ov_open(file,fd,mode) character*80 file integer*4 fd,mode
file is opened, and the file descriptor is passed back in fd. fd is passed as a handle to subsequent calls to ov_read(). mode is 0 for read-only, 1 for write-only, and 2 for read-write permission.
Reading an Overlay File
Once an overlay file is open, a data-set may be read in by a call to ov_read().
C Binding:
OverlayHdr* ov_read(fd,read_data) int fd; int read_data;
fd is the file descriptor returned by ov_open(). ov_read() returns a pointer to an OverlayHdr structure. If read_data is 0, the data is not read and the only the header information is returned. The data is skipped, and the file pointer is positioned at the beginning of the next header (if any) in the file. If read_data is 1, the data is read and a pointer to it is included in the returned header structure. If read_data is 0, the data is not read and the only the header information is returned. The data is skipped, and the file pointer is positioned at the beginning of the next header (if any) in the file. If read_data is 1, the data is read and a pointer to it is included in the returned header structure.
Fortran Binding:
ov_read(fd,ohdr,comment_len,comment,title_len,title, units_len,units,param_desc_len param_desc,private_size, private,data_size,data,vdata_size,vdata,data_flag,lat, lon) integer*4 fd record /OverlayHdr/ ohdr integer*4 comment_len character comment*(*) integer*4 title_len character title*(*) integer*4 units_len character units*(*) integer*4 param_desc_len character param_desc*(*) integer*4 private_size character private*(*) integer*4 data_size real data(*) integer*4 vdata_size real vdata(*) integer*4 data_flag real lat(*) real lon(*)
If an overlay data file contains several data-sets, subsequent calls to ov_read() will return consecutive data sets. Currently, there is no good way to determine the number of data sets in a file, or to skip over datasets in a file.
Writing an Overlay File
There are two operations involved with writing an overlay file:
To create an overlay file for output, a program should call ov_create().
C Binding:
ov_create(file) char *file;
file is the name of the output file. ov_create() returns a file descriptor that should be passed to ov_write() to write data sets to the file.
Fortran Binding:
ov_create(file,fd) character file*(*) integer*4 fd
file is created, and the file descriptor is passed back in fd. fd is passed as a handle to subsequent calls to ov_write().
Writing an Overlay File
After creating an output file, any number of data-sets can be written to it. Each data-set is written by a call to ov_write().
C Binding:
ov_write(fd,hdr) int fd; OverlayHdr *hdr;
fd is the file descriptor returned by ov_create(), and hdr is a pointer to the header structure that defines the data set to be written.
Fortran Binding:
ov_write(fd,ohdr,comment,title units,param_desc, private,data,vdata) integer*4 fd OverlayHdr ohdr character comment*(*) character title*(*) character units*(*) character param_desc*(*) character private*(*) real data(*) real vdata(*)
Note that the proper dimensions of these arrays must be set by the calling routine in the OverlayHdr structure.
Closing Overlay Files
C Binding:
Close overlay files with a call to the Unix file system function close().Fortran Binding:
Files created or opened from Fortran can be closed with a call to ov_close(fd) integer*4 fd
fd is the file descriptor passed back by ov_create() or ov_open().
Utility Routines
There are a few library routines available that operate on overlay header structures.
Printing an Overlay Header
To print the contents of an OverlayHdr structure to the standard output, call ov_print_hdr().
C Binding:
ov_print_hdr(hdr) OverlayHdr *hdr;
hdr is a pointer to the structure to be printed.
Fortran Binding:
ov_print_hdr(hdr) record /OverlayHdr/ ohdrFreeing Overlay Data
Call ov_free() to free all of the data that was allocated dynamically for a given overlay data set. This call frees the comment string, the private data, and the contour or vector data associated with the OverlayHdr structure. This is a handy routine if you need to reuse or destroy an OverlayHdr structure.
C Binding:
ov_free(hdr) OverlayHdr *hdr;Fortran Binding:
ov_free(hdr) record /OverlayHdr/ ohdrCreating an Overlay Header
ov_create_hdr() creates an OverlayHdr structure and returns a pointer to the structure. It fills in the structure with the minimum subset of values to make the structure usable.
C Binding:
OverlayHdr * ov_create_hdr(type,num_rows,num_columns,data,vdata) int type; int num_rows,num_columns; float *data,*vdata;
type may be OV_CONTOUR, OV_VECTOR, or OV_OUTLINE. num_rows and num_columns define the dimensions of the data arrays. data and vdata are pointers to the primary and secondary data arrays.
For types OV_OUTLINE, OV_UNGRIDDED_VECTOR, and OV_UNGRIDDED_SCALAR, num_rows contains the number of data points in the file, and num_columns is ignored (set it to zero).
EXAMPLES
Reading and Writing Satellite Image Files from C
The following program gives an example of how to read and write satellite image files using the C callable routines in the library. The program reads in two input files, combines them in some fashion and writes an output file. Note that the function used to combine the files can be changed by simply changing the line that defines the macro f(a,b). The line that defines the macro fdesc(f1,f2) should also be altered so that the comment field of the output header reflects the contents of the file.
/* this program applies a function point-by-point on two input files which share a lat-lon file. It outputs a new file. Each point in the output file is the result of the function applied to the corresponding points in the input files. i.e., out[i,j] = f(in1[i,j],in2[i,j]) where i is the scan #, and j is the sample # within the scan. Note this program will not work on input files with a variable number of samples in each scan, but could be modified to do so. We also assume that the two input files share a lat-lon file.*/ #includeReading and Writing Satellite Image Files from Fortran#include "si.h" /*define the function to be applied on each point here */ #define f(a,b) (a+b)/2.0/* average the two values */ #define fdesc(f1,f2) sprintf(desc,"average of %s and %s",f1,f2); #define TRUE 1 #define FALSE 0 char desc[200]; /*arguments are in_file1, in_file2, out_file */ main(argc,argv) int argc; char *argv[]; { SatImageHdr *in1_hdr,*in2_hdr,*out_hdr; int i,j; char *input_file1,*input_file2,*output_file; float *data1, *data2, *outdata; /*pointers to float arrays */ if(argc < 3) { print_usage(argv[0]); exit(); } input_file1 = argv[1]; input_file2 = argv[2]; output_file = argv[3]; /*open the input files */ if((in1_hdr = si_open(input_file1,0)) == (SatImageHdr *)-1) { fprintf("Error opening `%s'.\n",input_file1); exit(); } if((in2_hdr = si_open(input_file2,0)) == (SatImageHdr *)-1) { fprintf("Error opening `%s'.\n",input_file2); exit(); } /*the output hdr will look a lot like either of the input headers, so just copy one */ out_hdr = si_dup_header(in1_hdr,NULL); /*change the comment field to reflect the contents of the output file */ /*fdesc is a macro at the beginning of this file that can be easily changed*/ /* (note that sprintf returns pointer to output string) */ out_hdr->comment = fdesc(input_file1,input_file2); /*no matter the input file data_type, the output will be FLOAT*/ out_hdr->data_type = FLOAT; /*create the output file & write out its header. We won't create a lat-lon file since we can share the one used by the input files */ si_write_header(output_file,out_hdr,FALSE); /*allocate enough space to hold 1 scan of data from each file */ if((data1 = (float *)malloc(in1_hdr->samps_per_scan * sizeof(float))) < 0) { fprintf("malloc failed\n"); exit(); } if((data2 = (float *)malloc(in2_hdr->samps_per_scan * sizeof(float))) < 0) { fprintf("malloc failed\n"); exit(); } if((outdata = (float *)malloc(in2_hdr->samps_per_scan * sizeof(float))) < 0) { fprintf("malloc failed\n"); exit(); } /*read a scanline at a time from each input file */ for(i=0; i num_scans; i++) { /* we don't care about lat,lon,num_pts or time so pass NULLs */ si_read_scanline(in1_hdr,data1,NULL,NULL,NULL,NULL); si_read_scanline(in2_hdr,data2,NULL,NULL,NULL,NULL); /* perform the `function' on each point in the scanline. put the result in outdata */ for(j=0; j samps_per_scan; j++) { outdata[j] = f(data1[j],data2[j]); } /*write the scanline of data */ si_write_scanline(out_hdr,outdata,NULL,NULL,NULL,NULL); } /*close the files */ si_close(in1_hdr); si_close(in2_hdr); si_close(out_hdr); } print_usage(program_name) char *program_name; { fprintf(stderr,"Usage: %s input_file1 input_file2 output_file\n", program_name); }
This example program performs the same function as the C program above.
c This program performs a function point-by-point on two input c files c which share a lat-lon file. It outputs a new file. Each point c in the output file is the result of the function applied to the c corresponding points in the input files. c i.e., out[i,j] = f(in1[i,j],in2[i,j]) c where i is the scan #, and j is the sample # within the scan. c Note this program will not work on input files with a variable c number of samples in each scan, but could be modified to do so. c We also assume that the two input files share a lat-lon file. program si_combine include "/eos/katz/si/sif.h" character inputfile1*150,inputfile2*150,outputfile*150 record /SatImageHdr/ in1_hdr, in2_hdr, out_hdr integer*4 i,j,num_pts real*4 time real*4 data1(2000), data2(2000), outdata(2000) real*4 lat(2000), lon(2000) c prompt the user for the names of the files print *,'Enter name of first input file:' read(*,*)inputfile1 print *,'Enter name of second input file:' read(*,*)inputfile2 print *,'Enter name of output file:' read(*,*)outputfile c open the input files call si_open(inputfile1,0,in1_hdr) call si_open(inputfile2,0,in2_hdr) c the output hdr will look a lot like either of the input headers, c so just copy one call si_dup_header(in1_hdr,out_hdr) c no matter the input file data_type, the output will be FLOAT out_hdr.data_type = FLOAT c change the comment field to reflect the contents of the output file c ALTER THIS LINE IF YOU CHANGE CFUNC() call si_set_comment(out_hdr, `Averaged Data') c create the output file & write out its header. We won't create a lat-lon c file since we can share the one used by the input files call si_write_header(outputfile,out_hdr,0) c loop. read a scanline at a time from each input file. the lat-lon info c can be ignored. Likewise, num_pts & time are not used. if(in1_hdr.num_scans .gt. 2000)then print *,'Warning: too many samples per scan line' end if do 100 i=1,in1_hdr.num_scans call si_read_scanline(in1_hdr,data1,lat,lon,num_pts,time) call si_read_scanline(in2_hdr,data2,lat,lon,num_pts,time) c perform the `function' on each point in the scan line do 200 j=1,in1_hdr.samps_per_scan call cfunc(data1(j),data2(j),outdata(j)) 200 continue c write out the new scanline to the output file c since the latlon flag of the call to si_write_header() was zero, c the lat & lon arrays are not used by si_write_scanline(). call si_write_scanline(out_hdr,outdata,lat,lon,num_pts,time) 100 continue c close the files call si_close(in1_hdr) call si_close(in2_hdr) call si_close(out_hdr) stop end c change this subroutine to perform any function on the input c points. this function defines how the two files are combined. c currently, the two input values are simply averaged. subroutine cfunc(in1, in2, out) real*4 in1, in2, out out = (in1+in2)/2.0 return endReading and Writing Overlay Files in Fortran
This program reads in an overlay file, divides each value by 100, then writes a new overlay file as output.
program ov_test include "ovf.h" character infile*200, outfile*200 character comment*200,title*200,units*20,param_desc*200 real*4 data(100000), vdata(100000), private record /OverlayHdr/ hdr c make sure input and output files are named on the command line if(iargc() .ne. 2)then write(*,*) `Usage: ov_testf infile outfile' stop end if c get the input and output filenames from the command line call getarg(1, infile) call getarg(2, outfile) c open and read the input file call ov_open(infile,fd,0) call ov_read(fd,hdr,100,comment,200,title,20,units,200, & param_desc,0,private,100000,data,100000,vdata,1) c print out the file header structure call ov_print_hdr(hdr) call ov_close(fd) c perform the operation on the data call adjust_data(hdr,data) c if it is vector data, do the vdata as well if(hdr.type .eq. 2) call adjust_data(hdr,vdata) c write out the file call ov_create(outfile,fd) call ov_write(fd,hdr,comment,title,units,param_desc, & private,data,vdata) call ov_close(fd) stop end c divide each data value by 100.0 subroutine adjust_data(hdr,data) include "ovf.h" record /OverlayHdr/ hdr integer*4 i,j real data(*) do 200 i= 0, hdr.num_rows-1, 1 do 200 j= 0, hdr.num_columns-1, 1 index = i*hdr.num_columns + j data(index) = data(index)/100.0 200 continue return end