|
ocssw
1.0
|
00001 #include <stdio.h> 00002 #include <stdlib.h> 00003 #include <time.h> 00004 #include <string.h> 00005 00006 #include <clo.h> 00007 #include <hdf_utils.h> 00008 #include <genutils.h> 00009 00010 #define PROGRAM "ice2hdf" 00011 #define VERSION "1.2" 00012 #define FMT_L2HDF 2 00013 00014 #define DATA_CITATION "Meier, W., F. Fetterer, K. Knowles, M. Savoie, M. J. Brodzik. 2006, updated quarterly. Sea ice concentrations from Nimbus-7 SMMR and DMSP SSM/I passive microwave data. Boulder, Colorado USA: National Snow and Ice Data Center. Digital media." 00015 #define WEB_SITE "http://nsidc.org/data/docs/daac/nsidc0051_gsfc_seaice.gd.html" 00016 #define FTP_SITE "ftp://sidads.colorado.edu/pub/DATASETS/seaice/polar-stereo/nasateam/near-real-time" 00017 00018 #define HEADER_SIZE 300 00019 00020 #define NORTH_ROWS 448 00021 #define NORTH_COLS 304 00022 #define SOUTH_ROWS 332 00023 #define SOUTH_COLS 316 00024 00025 #define MAX_DATA 250 00026 #define POLAR_HOLE 251 00027 #define MISSING_DATA 255 00028 00029 // data file structure 00030 typedef struct data_file_t { 00031 char* filename; 00032 uint8 *dataPtr; 00033 int32 numRows; 00034 int32 numCols; 00035 int hasHeader; 00036 int32 year; 00037 int32 doy; 00038 } data_file_t; 00039 00040 /* -------------------------------------------------------------------- */ 00041 /* allocateDataFile */ 00042 /* -------------------------------------------------------------------- */ 00043 data_file_t* allocateDataFile(char* filename, int32 numRows, int32 numCols) { 00044 data_file_t* dataFile = (data_file_t*) malloc(sizeof(data_file_t)); 00045 00046 dataFile->filename = strdup(filename); 00047 dataFile->dataPtr = (uint8*) malloc(numRows * numCols); 00048 dataFile->numRows = numRows; 00049 dataFile->numCols = numCols; 00050 dataFile->hasHeader = 0; 00051 dataFile->year = 0; 00052 dataFile->doy = 0; 00053 00054 return dataFile; 00055 } 00056 00057 /* -------------------------------------------------------------------- */ 00058 /* readFile */ 00059 /* -------------------------------------------------------------------- */ 00060 void readFile(data_file_t* dataFile) { 00061 char header[HEADER_SIZE]; 00062 int result; 00063 FILE *fin; 00064 00065 fin = fopen(dataFile->filename, "rb"); 00066 if (fin == NULL) { 00067 fprintf(stderr, "-E- %s line %d: Could not open file, \"%s\" .\n", 00068 __FILE__, __LINE__, dataFile->filename); 00069 exit(EXIT_FAILURE); 00070 } 00071 00072 if (dataFile->hasHeader) { 00073 result = fread(header, 1, HEADER_SIZE, fin); 00074 if (result != HEADER_SIZE) { 00075 fprintf(stderr, 00076 "-E- %s line %d: Error reading file header, \"%s\" .\n", 00077 __FILE__, __LINE__, dataFile->filename); 00078 exit(EXIT_FAILURE); 00079 } 00080 00081 // 00082 // parse the header 00083 // 00084 int32 numRows = atoi(&header[12]); 00085 int32 numCols = atoi(&header[6]); 00086 dataFile->year = atoi(&header[102]); 00087 dataFile->doy = atoi(&header[108]); 00088 00089 if (numRows != dataFile->numRows) { 00090 fprintf(stderr, 00091 "-E- %s line %d: Bad number of rows(%d) in \"%s\" should be %d\n", 00092 __FILE__, __LINE__, numRows, dataFile->filename, 00093 dataFile->numRows); 00094 exit(EXIT_FAILURE); 00095 } 00096 00097 if (numCols != dataFile->numCols) { 00098 fprintf(stderr, 00099 "-E- %s line %d: Bad number of columns(%d) in \"%s\" should be %d\n", 00100 __FILE__, __LINE__, numCols, dataFile->filename, 00101 dataFile->numCols); 00102 exit(EXIT_FAILURE); 00103 } 00104 00105 } 00106 00107 // 00108 // read the data 00109 // 00110 result = fread(dataFile->dataPtr, 1, dataFile->numRows * dataFile->numCols, 00111 fin); 00112 if (result != dataFile->numRows * dataFile->numCols) { 00113 fprintf(stderr, "-E- %s line %d: Error reading file data, \"%s\" \n", 00114 __FILE__, __LINE__, dataFile->filename); 00115 exit(EXIT_FAILURE); 00116 } 00117 00118 } 00119 00120 /* -------------------------------------------------------------------- */ 00121 /* maskData */ 00122 /* Set the data to 255 anywhere the mask dataset is set to the */ 00123 /* given values. */ 00124 /* -------------------------------------------------------------------- */ 00125 void maskData(data_file_t* dataFile, data_file_t* maskFile, int numVals, 00126 int* vals) { 00127 if (dataFile->numRows != maskFile->numRows) { 00128 fprintf(stderr, 00129 "-E- %s line %d: data file \"%s\" and mask file \"%s\" have different number of Rows\n", 00130 __FILE__, __LINE__, dataFile->filename, maskFile->filename); 00131 exit(EXIT_FAILURE); 00132 } 00133 if (dataFile->numCols != maskFile->numCols) { 00134 fprintf(stderr, 00135 "-E- %s line %d: data file \"%s\" and mask file \"%s\" have different number of Columns\n", 00136 __FILE__, __LINE__, dataFile->filename, maskFile->filename); 00137 exit(EXIT_FAILURE); 00138 } 00139 00140 int i, j; 00141 for (i = 0; i < dataFile->numRows * dataFile->numCols; i++) { 00142 for (j = 0; j < numVals; j++) { 00143 if (maskFile->dataPtr[i] == vals[j]) { 00144 dataFile->dataPtr[i] = MISSING_DATA; 00145 } 00146 } 00147 } 00148 } 00149 00150 /* -------------------------------------------------------------------- */ 00151 /* fillHoles */ 00152 /* Fill the holes in the data by iterativly averaging the surrounding */ 00153 /* pixels. Assume anything over 250 is a hole. */ 00154 /* */ 00155 /* -------------------------------------------------------------------- */ 00156 void fillHoles(data_file_t* dataFile) { 00157 int delta = 1; // size of the box around the pixel 00158 00159 uint8 * data2; 00160 int holesFound = 1; 00161 int row, col; 00162 int r, c; 00163 int r1, r2; 00164 int c1, c2; 00165 int count; 00166 int sum; 00167 00168 data2 = (uint8*) malloc(dataFile->numRows * dataFile->numCols); 00169 if (data2 == NULL) { 00170 fprintf(stderr, "-E- %s line %d: Could not malloc.\n", __FILE__, 00171 __LINE__); 00172 exit(EXIT_FAILURE); 00173 } 00174 00175 while (holesFound) { 00176 holesFound = 0; 00177 00178 // copy the data so we don't key off of newly filled values 00179 memcpy(data2, dataFile->dataPtr, dataFile->numRows * dataFile->numCols); 00180 00181 for (row = 0; row < dataFile->numRows; row++) { 00182 for (col = 0; col < dataFile->numCols; col++) { 00183 00184 if (data2[row * dataFile->numCols + col] > MAX_DATA) { 00185 holesFound = 1; 00186 00187 // calc limits of the box around the pixel 00188 r1 = row - delta; 00189 if (r1 < 0) 00190 r1 = 0; 00191 r2 = row + delta; 00192 if (r2 >= dataFile->numRows) 00193 r2 = dataFile->numRows - 1; 00194 00195 c1 = col - delta; 00196 if (c1 < 0) 00197 c1 = 0; 00198 c2 = col + delta; 00199 if (c2 >= dataFile->numCols) 00200 c2 = dataFile->numCols - 1; 00201 00202 // intialize the average variables 00203 count = 0; 00204 sum = 0; 00205 00206 // sum surrounding pixels 00207 for (r = r1; r <= r2; r++) { 00208 for (c = c1; c <= c2; c++) { 00209 if (data2[r * dataFile->numCols + c] <= MAX_DATA) { 00210 count++; 00211 sum += data2[r * dataFile->numCols + c]; 00212 } 00213 } // for c1 to c2 00214 } // for r1 to r2 00215 00216 if (count > 1) { 00217 dataFile->dataPtr[row * dataFile->numCols + col] = sum 00218 / count; 00219 } 00220 00221 } // pixel is a hole 00222 } // for cols 00223 } // for rows 00224 00225 } // while holes 00226 00227 free(data2); 00228 } 00229 00230 /* -------------------------------------------------------------------- */ 00231 /* main */ 00232 /* -------------------------------------------------------------------- */ 00233 int main(int argc, char* argv[]) { 00234 clo_optionList_t * list; 00235 clo_option_t* option; 00236 00237 char tmpStr[2048]; 00238 int initialNumOptions; 00239 int numExtraOptions; 00240 00241 data_file_t *northData; 00242 data_file_t *southData; 00243 int northNumRegions; 00244 int* northRegionList; 00245 data_file_t *northRegion; 00246 int southNumRegions; 00247 int* southRegionList; 00248 data_file_t *southRegion; 00249 00250 int32 sd_id; 00251 int32 hdfStat; 00252 idDS ds_id; 00253 00254 char northFilename[1024]; 00255 char southFilename[1024]; 00256 char outFilename[1024]; 00257 char northRegionFilename[1024]; 00258 char southRegionFilename[1024]; 00259 00260 sprintf(tmpStr, "ice2hdf %s (%s %s)", VERSION, __DATE__, __TIME__); 00261 clo_setVersion(tmpStr); 00262 00263 list = clo_createList(); 00264 00265 strcpy(tmpStr, "Usage: ice2hdf [options] northFile southFile outFile\n\n"); 00266 strcat(tmpStr, 00267 "Program to convert the NSIDC ice data into an HDF file that l2gen can use.\n"); 00268 strcat(tmpStr, "Documentation - "); 00269 strcat(tmpStr, WEB_SITE); 00270 strcat(tmpStr, "\nData - "); 00271 strcat(tmpStr, FTP_SITE); 00272 strcat(tmpStr, "\n"); 00273 clo_setHelpStr(tmpStr); 00274 00275 clo_addOption(list, "fill_pole", CLO_TYPE_BOOL, "true", "Fill the Arctic polar hole with 100% ice"); 00276 strcpy(tmpStr, "Mask the values in the given\n northern regions\n"); 00277 strcat(tmpStr, " 0: Lakes, extended coast\n"); 00278 strcat(tmpStr, " 1: Non-regional ocean\n"); 00279 strcat(tmpStr, " 2: Sea of Okhotsk and Japan\n"); 00280 strcat(tmpStr, " 3: Bering Sea\n"); 00281 strcat(tmpStr, " 4: Hudson Bay\n"); 00282 strcat(tmpStr, " 5: Baffin Bay/Davis Strait/Labrador Sea\n"); 00283 strcat(tmpStr, " 6: Greenland Sea\n"); 00284 strcat(tmpStr, " 7: Kara and Barents Seas\n"); 00285 strcat(tmpStr, " 8: Arctic Ocean\n"); 00286 strcat(tmpStr, " 9: Canadian Archipelago\n"); 00287 strcat(tmpStr, " 10: Gulf of St. Lawrence\n"); 00288 strcat(tmpStr, " 11: Land\n"); 00289 strcat(tmpStr, " 12: Coast"); 00290 clo_addOption(list, "region_mask_north", CLO_TYPE_INT, "[0,11,12]", tmpStr); 00291 clo_addOption(list, "region_file_north", CLO_TYPE_IFILE, 00292 "$OCDATAROOT/common/region_n.msk", 00293 "\n Region mask file for the northern hemisphere"); 00294 00295 strcpy(tmpStr, "Mask the values in the given\n southern regions\n"); 00296 strcat(tmpStr, " 2: Weddell Sea\n"); 00297 strcat(tmpStr, " 3: Indian Ocean\n"); 00298 strcat(tmpStr, " 4: Pacific Ocean\n"); 00299 strcat(tmpStr, " 5: Ross Sea\n"); 00300 strcat(tmpStr, " 6: Bellingshausen Amundsen Sea\n"); 00301 strcat(tmpStr, " 11: Land\n"); 00302 strcat(tmpStr, " 12: Coast"); 00303 clo_addOption(list, "region_mask_south", CLO_TYPE_INT, "[11,12]", tmpStr); 00304 clo_addOption(list, "region_file_south", CLO_TYPE_IFILE, 00305 "$OCDATAROOT/common/region_s.msk", 00306 "\n Region mask file for the southern hemisphere"); 00307 00308 initialNumOptions = clo_getNumOptions(list); 00309 clo_readArgs(list, argc, argv); 00310 numExtraOptions = clo_getNumOptions(list) - initialNumOptions; 00311 00312 if (numExtraOptions != 3) { 00313 printf("ERROR - wrong number of parameters given\n\n"); 00314 clo_printUsage(list); 00315 return 1; 00316 } 00317 00318 option = clo_getOption(list, initialNumOptions); 00319 parse_file_name(option->key, northFilename); 00320 00321 option = clo_getOption(list, initialNumOptions + 1); 00322 parse_file_name(option->key, southFilename); 00323 00324 option = clo_getOption(list, initialNumOptions + 2); 00325 parse_file_name(option->key, outFilename); 00326 00327 northRegionList = clo_getInts(list, "region_mask_north", &northNumRegions); 00328 parse_file_name(clo_getString(list, "region_file_north"), 00329 northRegionFilename); 00330 00331 southRegionList = clo_getInts(list, "region_mask_south", &southNumRegions); 00332 parse_file_name(clo_getString(list, "region_file_south"), 00333 southRegionFilename); 00334 00335 time_t currentTime = time(NULL); 00336 char* currentTimeStr = strdup(asctime(gmtime(¤tTime))); 00337 currentTimeStr[24] = '\0'; 00338 00339 // 00340 // read the data files 00341 // 00342 northData = allocateDataFile(northFilename, NORTH_ROWS, NORTH_COLS); 00343 northData->hasHeader = 1; 00344 readFile(northData); 00345 00346 if(clo_getBool(list, "fill_pole")) { 00347 int i; 00348 for(i=0; i<NORTH_ROWS * NORTH_COLS; i++) { 00349 if(northData->dataPtr[i] == POLAR_HOLE) { 00350 northData->dataPtr[i] = MAX_DATA; 00351 } 00352 } 00353 } 00354 00355 southData = allocateDataFile(southFilename, SOUTH_ROWS, SOUTH_COLS); 00356 southData->hasHeader = 1; 00357 readFile(southData); 00358 00359 if (northNumRegions > 0) { 00360 northRegion = allocateDataFile(northRegionFilename, NORTH_ROWS, 00361 NORTH_COLS); 00362 northRegion->hasHeader = 1; 00363 readFile(northRegion); 00364 maskData(northData, northRegion, northNumRegions, northRegionList); 00365 } 00366 00367 if (southNumRegions > 0) { 00368 southRegion = allocateDataFile(southRegionFilename, SOUTH_ROWS, 00369 SOUTH_COLS); 00370 southRegion->hasHeader = 1; 00371 readFile(southRegion); 00372 maskData(southData, southRegion, southNumRegions, southRegionList); 00373 } 00374 00375 fillHoles(northData); 00376 fillHoles(southData); 00377 00378 // 00379 // open the HDF file 00380 // 00381 sd_id = SDstart(outFilename, DFACC_CREATE); 00382 if (sd_id == FAIL) { 00383 fprintf(stderr, "-E- %s line %d: Could not create HDF file, %s .\n", 00384 __FILE__, __LINE__, outFilename); 00385 return (HDF_FUNCTION_ERROR); 00386 } 00387 00388 /* */ 00389 /* Create the north and south data array SDSes */ 00390 /* ---------------------------------------------------------------- */ 00391 /* */ 00392 PTB( CreateSDS( sd_id, /* file id */ 00393 "north", /* short name */ 00394 "North Pole Ice Fraction", /* long name */ 00395 NULL, /* standard name */ 00396 NULL, /* units */ 00397 0,0, /* range */ 00398 0,0, /* slope */ 00399 DFNT_INT8, /* HDF number type */ 00400 2, /* rank */ 00401 northData->numRows, northData->numCols, 0, /* dimension sizes */ 00402 "north_rows", "north_cols", NULL /* dimension names */ 00403 )); 00404 00405 PTB( CreateSDS( sd_id, /* file id */ 00406 "south", /* short name */ 00407 "South Pole Ice Fraction", /* long name */ 00408 NULL, /* standard name */ 00409 NULL, /* units */ 00410 0,0, /* range */ 00411 0,0, /* slope */ 00412 DFNT_INT8, /* HDF number type */ 00413 2, /* rank */ 00414 southData->numRows, southData->numCols, 0, /* dimension sizes */ 00415 "south_rows", "south_cols", NULL /* dimension names */ 00416 )); 00417 00418 ds_id.fid = sd_id; 00419 ds_id.fftype = FMT_L2HDF; 00420 00421 // Modify to use new HDF4/NCDF4 libhdfutils library JMG 09/28/12 00422 PTB( SetChrGA (ds_id, "Title", "Sea Ice Concentration Daily")); 00423 PTB( SetChrGA (ds_id, "Creation Date", currentTimeStr )); 00424 PTB( SetChrGA (ds_id, "Created By", "Don Shea, SAIC, from NSIDC NRT data (http://nsidc.org)")); 00425 PTB( SetChrGA (ds_id, "Reference1", DATA_CITATION )); 00426 PTB( SetChrGA (ds_id, "Software Name", PROGRAM)); 00427 PTB( SetChrGA (ds_id, "Software Version", VERSION)); 00428 PTB( SetI32GA (ds_id, "year", northData->year)); 00429 PTB( SetI32GA (ds_id, "day_of_year", northData->doy)); 00430 00431 PTB( sd_writedata(sd_id, "north", northData->dataPtr, 0, 0, 0, northData->numRows, northData->numCols,0)); 00432 PTB( sd_writedata(sd_id, "south", southData->dataPtr, 0, 0, 0, southData->numRows, southData->numCols,0)); 00433 00434 // 00435 // close the HDF file 00436 // 00437 hdfStat = SDend(sd_id); 00438 if (hdfStat == FAIL) { 00439 fprintf(stderr, "-E- %s line %d: Could not close HDF file, %s .\n", 00440 __FILE__, __LINE__, outFilename); 00441 return (HDF_FUNCTION_ERROR); 00442 } 00443 00444 return (LIFE_IS_GOOD); 00445 }
1.7.6.1