MODIS Ocean Level 3 Data读取程序.doc

一、MODIS Ocean Level 3 Binned Data读取程序IDL语言程序1、 从HDF文件中读取SDSFunction Get_SDS, file_name, SDS_Name If (N_Elements( file_name ) eq 0) then Message, Must specify filename. Id = HDF_SD_Start( file_name, /READ ) Index = HDF_SD_NameToIndex( Id, SDS_Name ) SDS = HDF_SD_Select( Id, Index ) HDF_SD_GetData, SDS, Data HDF_SD_End, Id Return, DataEnd2、 从一个byte中提取1个或多个bit的资料; The IDL function, B, extracts bits from an array of byte values.; i = array of byte values; j = starting bit position; k = number of bits to be extracted; Kevin Turpie, SAIC GSOFunction Bit, i, j, k If (N_Elements( k ) eq 0) then k = 1 k = Fix( k ) If (N_Elements( j ) eq 0) then j = 1 j = Fix( j ) Return, IShft ( IShft ( Long (i), 32-j-k+1 ), -(31-k+1) )End3、 将二进制资料转换成经纬度坐标; This code converts bin number to (lat,lon) ; It is written in both C and IDL.; The IDL routine, , calls a C handler, binloc.c, which in turn; calls the bin2ll routine from a shareable library, created from bin_csub.c; (originally from SeaDAS). Compilation instructions in the binloc.c file; work for SGI. references a .so library file created using those; instructions. The user must be careful to correctly specify the path and; name of this library file for their system in .;- IDL Interface: -;-; BinLoc; This IDL procedure takes an array of L3 bin numbers (b) and returns the; corresponding latitude and longitude values.; ARGUMENTS:; b = corresponding array of bin numbers; Lat = lat/lon grid containing sum of values; Lon = lat/lon grid containing number of values; KEYWORDS:; NROWS = number of rows in input L3 bin grid.; default is 4320 for MODIS. for SeaWiFS, set to 2160.; AUTHOR: Kevin Turpie, 08 May 2000; SAIC General Sciences Corporation; ;- PRO BinLoc, b, Lat, Lon, NROWS=nrows;On_Error, 2 If (N_Elements( nrows ) eq 0) then nrows = 4320L nb = N_Elements( b ) b = Long( b ) nbins = nb nl = N_Elements( Lat ) If (nl eq 0) then begin Lat = FltArr( nb ) nl = nb EndIf If (nl ne nb) then Message, Bin and Lat array must be the same size. Lat = Float( Lat ) nl = N_Elements( Lon ) If (nl eq 0) then begin Lon = FltArr( nb ) nl = nb EndIf If (nl ne nb) then Message, Bin and Lon array must be the same size. Lon = Float( Lon ) iret = Call_External( binloc.so, binloc, $ b, Long( nbins ), Long( nrows), Lat, Lon ) END#include /* binloc (private)!C-* Description: binloc is a routine that can be called from IDL to convert L3 bin numbers into latitude and longitude coordinates. In IDL, call_external is used to pass information via the arguments argc and argv (see below). An IDL routine called was created to handle this routine. Input Parameters: argc number of arguments. argv array of input and output parameters passed by value by call_external in IDL. The elements are described below. Element I/O Description argv0 I corresponding array of bin numbers argv1 I total number of bins argv2 I number of rows in bin grid argv3 O array containing latitude values argv4 O array containing longitude values Output Parameters: void Revision History: Revision 01.00 2000/05/08 17:22:00 K. TURPIE () AUTHOR Team-unique Header: References and Credits: Design Notes: Uses the SeaDAS routines bin_init and bin2ll in the source bin_csub.c or library libsdsbin.a To compile use: cc -32 -shared -KPIC -o binloc.so binloc.c libsdsbin.a or cc -32 -KPIC -c binloc.c cc -32 -KPIC -c bin_csub.c cc -32 -shared -o binloc.so binloc.o bin_csub.o or cc -32 -KPIC -c binloc.c cc -32 -shared -o binloc.so binloc.o libsdsbin.a!END-*/-C Code follows:-binloc.c (C handling for IDL routine binloc; calls bin2ll to do the job.)void binloc( int argc, void *argv ) long i, j, k, a, *b, nbins, nrows, *nbin, *bbin, tbin; float *lbin, *lat, *lon; b = (long *) argv0; nbins = *(long *) argv1; nrows = *(long *) argv2; lat = (float *) argv3; lon = (float *) argv4; bin_init( nrows, &nbin, &bbin, &lbin, &tbin ); for ( k=0; k nbins; k+ ) bin2ll( bk, &(latk), &(lonk) );bin_csub.c is taken from SeaDAS. bin2ll below does the conversion./* This file contains all subroutines for the seawifs level 3 binning*/#include #include #include #definePI3.141592653589793/* global variables only used inside this file, the call to bin_init may pass these variables to the calling routineNUMROWS : total number of rows for binning *NUMBIN: no. of bin in each row *BASEBIN: 1st bin of each row TOTBINS: total bin numbers *LATBIN: center latitude of each row */static long *NUMBIN, *BASEBIN;static float *LATBIN;static long TOTBINS;static long NUMROWS;static float SEAM_LON;/* bin_init: given the total row number, this subroutine returns above variables.*/int bin_init(nrow, nbin, bbin, lbin, tbin)long nrow;long *nbin, *bbin, *tbin;float *lbin; int i; double radfac; static int first=1; if (first = 0) /* free the memory if not the first call */ free(NUMBIN); free(BASEBIN); free(LATBIN); SEAM_LON = -180.0; /* this value should be passed in */ NUMROWS = nrow; NUMBIN = (long *) calloc(NUMROWS, sizeof(long); BASEBIN = (long *) calloc(NUMROWS, sizeof(long); LATBIN = (float *) calloc(NUMROWS, sizeof(float); radfac = PI / 180.0; for (i=0; iNUMROWS; i+) *(LATBIN+i) = (i + 0.5) * (180.0 / NUMROWS) - 90.0; *(NUMBIN+i) = (long) (cos(*(LATBIN+i)*radfac) * (2.0 * NUMROWS) + 0.5); *BASEBIN = 1; for (i=1; i 0 & *(BASEBIN+old_row-1) bin) row = old_row; else if (bin TOTBINS) bin = TOTBINS; /* north pole */ /* binary search for row in range 1.NUMROWS */ rlow = 1; /* 1-relative */ rhi = NUMROWS; /* 1-relative */ while (1) rmid = (rlow + rhi - 1) / 2; /* 0-relative */ if (*(BASEBIN+rmid) bin) rhi = rmid; else rlow = rmid + 1; if (rlow = rhi) row = rlow; break; old_row = row; *lat = *(LATBIN+row-1); *lon = 360.0 * (bin - *(BASEBIN+row-1) + 0.5) / *(NUMBIN+row-1); *lon = *lon + SEAM_LON; /* note, *lon returned here may be in 0 to 360 */* given the lat/lon, return the bin number lon has to be in the range of -180.0 to 180.0*/int ll2bin(lat, lon, bin)float lat, lon;long *bin; /* 1-relative */ long row, col; /* 0-relative */ row = (long) (90.0 + lat) * (float) NUMROWS / 180.0); col = (long) (float) (*(NUMBIN+row) * (lon - SEAM_LON) / 360.0); *bin = *(BASEBIN+row) + col;/* given the lat/lon, return the row, column lon has to be in the range of -180.0 to 180.0 6/96. Cast double in order to get same bin #s for SeaWiFS smap9 and SeaDAS smigen.*/int ll2rc(lat, lon, row, col)float lat, lon;long *row, *col; /* 1-relative */ *row = (long) (90.0 + (double) lat) * (double) NUMROWS / 180.0); *col = (long) (double) (*(NUMBIN + (*row) * (double) lon - SEAM_LON) / 360.0); *row = *row + 1; *col = *col + 1; /* given row/column, return lat/lon*/int rc2ll(row, col, lat, lon)long row, col; /* 1-relative */float *lat, *lon; *lat = *(LATBIN+row-1); *lon = SEAM_LON + (360.0 * (col - 0.5) / *(NUMBIN+row-1);/* given a row/column number, return the bin number (1-relative) */int rc2bin(row, col, bin)long row, col, *bin; *bin = *(BASEBIN+row-1) + col - 1;/* given a bin number, return the row and column (both are 1-relative) heuristic and binary search algorithm is used*/int bin2rc(bin, row, col)long bin; /* 1-relative */long *row, *col; /* 1-relative */ long rlow, rhi, rmid; static long old_row=0; /* 1-relative */ if (old_row 0 & *(BASEBIN+old_row-1) bin) *row = old_row; else if (bin TOTBINS) bin = TOTBINS; /* north pole */ /* binary search for row in range 1.NUMROWS */ rlow = 1; /* 1-relative */ rhi = NUMROWS; /* 1-relative */ while (1) rmid = (rlow + rhi - 1) / 2; /* 0-relative */ if (*(BASEBIN+rmid) bin) rhi = rmid; else rlow = rmid + 1; if (rlow = rhi) *row = rlow; break; old_row = *row; *col = bin - *(BASEBIN + (*row) - 1) + 1;/* given a bin number, return the center lat/lon of that bin number no heuristic or binary search algorithm is used this routine is very slow due to the array reference (I think)*/*int old_bin2ll(bin, lat, lon)long bin;float *lat, *lon; long row; row = NUMROWS-1; while (bin BASEBINrow) row-; *lat = LATBINrow; *lon = 360.0 * (bin - BASEBINrow + 0.5) / NUMBINrow;*/* update version, much faster than using array reference */int old_bin2ll(bin, lat, lon)long bin;float *lat, *lon; long row; long *tmpptr; row = NUMROWS - 1; tmpptr = BASEBIN + NUMROWS - 1; while (bin 0)+2 ) qkey = StrMid( qkey, 0, Mpos )+Q+StrMid( qkey, Mpos+1, lq-Mpos )+* If ( (StrMid( qkey, (pathpos0)+3, 2 ) eq AP) or $ (StrMid( qkey, (pathpos0)+3, 2 ) eq SP) ) then begin; Its an