@read_standard_hdf pro nb_cloud_mask_modis_prior input_path = "/apollo/cloud/archive/Validation_Data/MYD08_M3/" dlon = 1.0 dlat = 1.0 lat_min = -90.0 lat_max = 90.0 lon_min = -180.0 lon_max = 180.0 nlon = (lon_max-lon_min)/dlon nlat = (lat_max-lat_min)/dlat nmonth = 12 ntimes = 3 missing_input = -999.0 missing_output = -999.0 count_min = 10.0 count_table = make_array(nlon, nlat, nmonth, ntimes, /FLOAT, VALUE = 0.0) cld_frac_sum_table = make_array(nlon, nlat, nmonth, ntimes, /FLOAT, VALUE = 0.0) cld_frac_table = make_array(nlon, nlat, nmonth, ntimes, /FLOAT, VALUE=missing_output) multi_sum_table = make_array(nlon, nlat, nmonth, ntimes, /FLOAT, VALUE = 0.0) multi_frac_table = make_array(nlon, nlat, nmonth, ntimes, /FLOAT, VALUE=missing_output) count_2d = make_array(nlon, nlat, /FLOAT, VALUE=0) files = file_search(input_path,'*.hdf') n_files = n_elements(files) input_length = strlen(input_path) for i_file = 0, n_files - 1 do begin root_name = strmid(files[i_file],input_length) temp = strsplit(root_name, '.', /extract) date_string = temp[1] year = fix(strmid(date_string, 1,4)) doy = fix(strmid(date_string, 5,3)) compute_month_day, doy, leapyr(year), month, day print, 'processing ', i_file, ' of ', n_files, root_name, year, doy, month ;--- read data cfrac_mean_mean_temp = read_standard_hdf( files[i_file],'Cloud_Fraction_Mean_Mean') cfrac_day_mean_temp = read_standard_hdf( files[i_file],'Cloud_Fraction_Day_Mean_Mean') cfrac_night_mean_temp = read_standard_hdf( files[i_file],'Cloud_Fraction_Night_Mean_Mean') k = month - 1 idx = where(cfrac_mean_mean_temp ne missing_input, complement = idx_missing, ncomplement=cc_missing,cc) d = 0 count_2d[*,*] = 1 if (cc_missing gt 0) then cfrac_mean_mean_temp[idx_missing] = 0.00 if (cc gt 0) then begin count_2d[idx] = 1 count_table[*,*,k,d] = count_table[*,*,k,d] + count_2d cld_frac_sum_table[*,*,k,d] = cld_frac_sum_table[*,*,k,d] + cfrac_mean_mean_temp endif idx = where(cfrac_day_mean_temp ne missing_input, complement = idx_missing, ncomplement=cc_missing,cc) d = 1 count_2d[*,*] = 1 if (cc_missing gt 0) then cfrac_day_mean_temp[idx_missing] = 0.00 if (cc gt 0) then begin count_2d[idx] = 1 count_table[*,*,k,d] = count_table[*,*,k,d] + count_2d cld_frac_sum_table[*,*,k,d] = cld_frac_sum_table[*,*,k,d] + cfrac_day_mean_temp endif idx = where(cfrac_night_mean_temp ne missing_input, complement = idx_missing, ncomplement=cc_missing,cc) d = 2 count_2d[*,*] = 1 if (cc_missing gt 0) then cfrac_night_mean_temp[idx_missing] = 0.00 if (cc gt 0) then begin count_2d[idx] = 1 count_table[*,*,k,d] = count_table[*,*,k,d] + count_2d cld_frac_sum_table[*,*,k,d] = cld_frac_sum_table[*,*,k,d] + cfrac_night_mean_temp endif cycle: endfor idx = where(count_table gt count_min,cc) if (cc gt 0) then begin cld_frac_table[idx] = cld_frac_sum_table[idx] / count_table[idx] endif cld_frac_table_smoothed = smooth(cld_frac_table,[3,3,2,1],/Edge_Mirror,Missing=missing_output) ;multi_frac_table_smoothed = smooth(multi_frac_table,[3,3,1,1],/Edge_Mirror,Missing=missing_output) ;window, 0, xsize = 800, ysize = 600 ;load_color_table, 33, 0, 1.0, 249 ;erase, color = 3 ;map_set, 0.0, 0.0, limit=[-90,-180,90,180], position = [0.05,0.15,0.95,0.95], /noerase ;view, reform(cld_frac_table[*,*,0,2]), coltab=13, sds_min=-101, sds_max=1.0, title = 'Night Total Cloud Fraction', /current ;map_continents, color = 3, thick = 2, /lores ;view, reform(cld_frac_table[*,*,0,0]), coltab=13, sds_min=-101, sds_max=1.0, title = 'Day+Night Total Cloud Fraction' ;view, reform(cld_frac_table[*,*,0,1]), coltab=13, sds_min=-101, sds_max=1.0, title = 'Day Total Cloud Fraction' ;view, reform(cld_frac_table[*,*,0,2]), coltab=13, sds_min=-101, sds_max=1.0, title = 'Night Total Cloud Fraction' ;view, reform(multi_frac_table[*,*,0,0]), coltab=13, sds_min=-101, sds_max=1.0, title = 'Day+Night Multilayer Fraction' ;view, reform(multi_frac_table[*,*,0,1]), coltab=13, sds_min=-101, sds_max=1.0, title = 'Day Multilayer Fraction' ;view, reform(multi_frac_table[*,*,0,2]), coltab=13, sds_min=-101, sds_max=1.0, title = 'Night Multilayer Fraction' save, file = 'modis_prior.sav' ;--------------------------------- nc_filename = 'nb_cloud_mask_modis_prior.nc' header = timestamp() sd_id = ncdf_create(nc_filename,/CLOBBER) ncdf_attput,sd_id, "description", header, /CHAR, /global ncdf_attput,sd_id, "number_longitudes", float(nlon), /FLOAT, /global ncdf_attput,sd_id, "number_latitudes", float(nlat), /FLOAT, /global ncdf_attput,sd_id, "number_months", float(nmonth), /FLOAT, /global ncdf_attput,sd_id, "number_times", float(ntimes), /FLOAT, /global ncdf_attput,sd_id, "latitude_min", float(lat_min), /FLOAT, /global ncdf_attput,sd_id, "latitude_max", float(lat_max), /FLOAT, /global ncdf_attput,sd_id, "longitude_min", float(lon_min), /FLOAT, /global ncdf_attput,sd_id, "longitude_max", float(lon_max), /FLOAT, /global ncdf_attput,sd_id, "longitude_spacing", float(dlon), /FLOAT, /global ncdf_attput,sd_id, "latitude_spacing", float(dlat), /FLOAT, /global ;--- get dimensions nx_id = NCDF_DIMDEF(sd_id, 'longitude', nlon) ny_id = NCDF_DIMDEF(sd_id, 'latitude', nlat) nm_id = NCDF_DIMDEF(sd_id, 'month', nmonth) nd_id = NCDF_DIMDEF(sd_id, 'diurnal', ntimes) ;--- create variable sds_id1 = NCDF_VARDEF(sd_id, 'cloud_fraction_table', [nx_id,ny_id,nm_id,nd_id], /FLOAT) NCDF_CONTROL, sd_id, /ENDEF ; Put the file into data mode ;---------------------------------------------- ; write to netcdf file ;---------------------------------------------- NCDF_VARPUT, sd_id, sds_id1, cld_frac_table ;--- close netcdf files ncdf_close, sd_id stop end ;=================================================================== ; routine to read hdf sds's ;=================================================================== pro read_hdf, filename, sds_name, sds_data, Error_Status catch, Error_Status if (Error_Status ne 0) then begin print, "error in read_hdf, no ",sds_name sds_data = -1 return endif sd_id = hdf_sd_start(filename, /READ) sds_id = hdf_sd_select(sd_id, hdf_sd_nametoindex(sd_id,sds_name)) hdf_sd_getdata, sds_id, sds_data hdf_sd_endaccess, sds_id hdf_sd_end, sd_id catch, /CANCEL end