actual_albedo_hdf.zip_7KC_actual_albedo_albedo_hdf_反照率

/* !C*********************************************************************** !DESCRIPTION: This program reads MODIS BRDF/albedo product MOD43A1, computes the black and white sky albedos (with Wolfgang's simple polynomial equation) and then assuming optical depth is an input, computes the actual albedo (linear combination of white and black-sky albedo) and prints it out in HDF format. !INPUTS: MOD43A1 Optical depth Solar zenith angle !OUTPUTS: Actual albedo (blue sky albedo) !DEVELOPERS: Feng Gao, Crystal Schaaf (Boston University) [email protected], [email protected] !REVISIONS: 1/2000, 5/2000, 6/2015 !END ******************************************************************* */ #include <stdio.h> #include <math.h> #include <string.h> #include <assert.h> #include "mfhdf.h" #define SUCCESS 0 #define FAILURE -1 #define Max_Str_Len 100 #ifndef PI #define PI 3.1415926535 #endif #define D2R PI/180.0 #define ALBEDO_fillv 32767 #define ALBEDO_factor 0.001 #define ALBEDO_offset 0 char m_band_names[10][20] = {"Band1","Band2","Band3","Band4", "Band5","Band6","Band7","vis","nir","shortwave"}; typedef struct { float iso; float vol; float geo; } PARAMETERS; typedef struct { int aerosol_type; int bandno; float solar_zenith; float optical_depth; float skyl; } SKYL; typedef struct { int nrows; int ncols; int nbands; int16 ***par_data; float64 scale; float64 offset; int16 fillv; int32 hdf_id; int32 sds_id[10]; } BRDF_SDS; typedef struct { int nrows; int ncols; int nbands; int16 ***albedo_data; float64 scale; float64 offset; int16 fillv; int32 hdf_id; int32 sds_id[10]; } ALBEDO_SDS; void usage(void); void read_skyl_table (float skyl_lut[2][10][90][50]); float get_skyl (SKYL s,float skyl_lut[2][10][90][50]); float cal_bsa(PARAMETERS p, SKYL s); float cal_wsa(PARAMETERS p); float cal_actual_albedo(PARAMETERS p, SKYL s); int readSDSattribute(BRDF_SDS *b); void writeSDSattribute(ALBEDO_SDS *a); void alloc_1dim_contig (void **ptr,int d1,int elsize); void alloc_3dim_contig (void ****ptr,int d1,int d2,int d3,int elsize); void CloseAlbedoFile(ALBEDO_SDS *a); void CloseInput(BRDF_SDS *b); char parname[Max_Str_Len]; char outname[Max_Str_Len]; void main(int argc, char **argv) { int i,j,band; float albedo; float skyl_lut[2][10][90][50]; intn status; int32 sd_id,sds_id,sds_index,sdout_id,sdsout_index,sdsout_id; int32 n_datasets,n_file_attrs,index; int32 rank,data_type,n_attrs; char name[MAX_NC_NAME]; ALBEDO_SDS a; BRDF_SDS b; PARAMETERS p; SKYL s; // MODFILE *par,*out; /* see if commmand line are correct or not */ if(argc!=9) { usage(); exit(1); } /* parse command line */ for(i=1;i<argc;i++){ if(strcmp(argv[i],"-par")==0) strcpy(parname,argv[++i]); else if(strcmp(argv[i],"-out")==0) strcpy(outname,argv[++i]); else if(strcmp(argv[i],"-od")==0) s.optical_depth=atof(argv[++i]); else if(strcmp(argv[i],"-szn")==0) s.solar_zenith=atof(argv[++i]); else{ printf("\nWrong option:%s\n",argv[i]); usage(); exit(1); } } printf("\nRead SKYL table \n"); read_skyl_table(skyl_lut); /* open BRDF parameters HDF-EOS file and get factor scale and offset*/ printf("Read BRDF SDS attribute\n"); /* read neccessary BRDF parameters SDS attributes */ if(0 != readSDSattribute(&b)){ printf("Read hdf info failed: %s\n", parname); exit (1); } /* initialize albedo SDS attribute */ a.fillv=ALBEDO_fillv; a.scale=ALBEDO_factor; a.offset=ALBEDO_offset; a.nrows=b.nrows; a.ncols=b.ncols; a.nbands=b.nbands; printf("Write ALbedo SDS attribute\n"); writeSDSattribute(&a); printf("Allocate memory\n"); /* allocate all memory */ //alloc_3dim_contig((void ****)&a.albedo_data,1,a.ncols,a.nbands,sizeof(int16)); //alloc_3dim_contig((void ****)&b.par_data,b.ncols,b.nbands,3,sizeof(int16)); int16 (*buf_par)[3] = (int16 (*)[3])malloc(sizeof(int16)*b.ncols*3); int16 *buf_alb = (int16 *)malloc(sizeof(int16)*b.ncols); printf("Computing actual albedo\n"); /* to process each line */ int32 start[3] = {0, 0, 0}; int32 edge[3] = {1, b.ncols, 3}; for(band=0; band<b.nbands; band++){ for(i=0; i<b.nrows; i++){ start[0]=i; status =SDreaddata(b.sds_id[band], start, NULL ,edge,buf_par); assert(status == SUCCEED); /* process each pixel */ for(j=0; j<b.ncols; j++){ /* check parameters and see if they are valid values */ if(buf_par[j][0]==b.fillv ) { buf_alb[j] = a.fillv; continue; } /* calculate real parameter values */ p.iso=(float)(buf_par[j][0]-b.offset)*b.scale; p.vol=(float)(buf_par[j][1]-b.offset)*b.scale; p.geo=(float)(buf_par[j][2]-b.offset)*b.scale; /* NOTE: need to make sure you get the right LUT band order: # bandwidth: 0: 0.620-0.670 (red) # 1: 0.841-0.876 (nir) # 2: 0.459-0.479 # 3: 0.545-0.565 # 4: 1.230-1.250 # 5: 1.628-1.652 # 6: 2.105-2.155 # BB 7: 0.400-0.700 (vis) # BB 8: 0.700-4.000 (nir) # BB 9: 0.250-4.000 (sw) you may need index your band order to LUT band order, such as: s.bandno=index[band]; */ s.bandno=band; /* assume continental type, but finally need to get it from BRDF_type 0: continental 1: maritime */ s.aerosol_type=0; /* get SKYL from SKYL lookup table (it depends on optical depth, solar zenith angle, aerosol type, and bands) */ s.skyl=get_skyl(s,skyl_lut); albedo=cal_actual_albedo(p,s); buf_alb[j]=albedo/a.scale+a.offset; } // for col /* write albedo data */ status=SDwritedata(a.sds_id[band], start, NULL, edge, buf_alb); HEprint(stdout, 0); assert(status == SUCCEED); } // for row } // for band CloseInput(&b); CloseAlbedoFile(&a); printf("Finished\n\n"); } void usage(void) { printf("\nUsage: actual_albedo.exe [-par][-od][-szn][-out] \n\n"); printf(" -par <par_file> input MCD43A1 (V005 or V006, in HDF format)\n"); printf(" -od <optical_depth> input optical depth (float, range: 0.0-1.0)\n"); printf(" -szn <solar_zenith> input solar zenith angle you want to compute\n"); printf(" (float, range: 0.0-89.0 degrees)\n"); printf(" -out <albedo_file> output file to save actual albedo (in HDF format)\n\n"); } /******************************************************* read SKYL LUT from pre-generated LUT file (skyl_lut.dat) 2 aerosol types 10 bands (7 MODIS bands + 3 broad bands (VIS, NIR, SW) 90 degrees (0-89 degrees with 1 degree step 50 optical depth ( 0-1 with 0.02 step) ********************************************************/ void read_skyl_table (float skyl_lut[2][10][90][50]) { char str[Max_Str_Len]; FILE *in; int aerosol,band,szn,od,NBANDS; NBANDS=10; /* MODIS defaults */ if((in=fopen("skyl_lut.dat","r"))==NULL) { printf("Can't open SKYL LUT file (skyl_lut.dat) \n"); exit(1); } for(aerosol=0;aerosol<2;aerosol++) { fscanf(in,"%s %s\n",str,str); for(band=0;band<NBANDS;band++) { fscanf(in,"%s %s\n",str,str); for(od=0;od<51;od++) fscanf(in,"%s ",str); for(szn=0;szn<90;szn++){ fscanf(in,"%s ",str); for(od=0;od<50;od++) { fscanf(in,"%f ",&(skyl_lut[aerosol][band][szn][od])); /*printf("