pro IDL_VIIRS_WaterMap_Combined, Inpath, Outpath,keywords ;-------Check IDL Environment for NC file---------- if(NCDF_EXISTS() eq 0) then begin print, 'NC is not supported by the current IDL environment' return endif ;-------Check IDL Environment for HDF file------------ if(HDF_EXISTS() eq 0) then begin print, 'HDF is not supported by the current IDL environment' return endif Result = FILE_SEARCH(Inpath + keywords, /FULLY_QUALIFY_PATH) RFlag = Keyword_set(Result) IF RFlag eq 0 then begin print, 'No hdf file exists in your file directory' return ENDIF NFile = N_elements(Result) IF NFile gt 0 then begin FileNames = Result indx = Sort(FileNames) FileNames = FileNames[indx] print, 'Input hdf file: ', FileNames EndIF Else Begin return EndElse for ifile = 0, NFile - 1 do begin filename = FileNames[ifile] print, filename publicname = strmid(filename,0, strlen(filename) - 4) publicname_1 = strmid(filename,strlen(Inpath), strlen(filename) - strlen(Inpath) - 4) ;---check whether the png file for the current data has been created---- png_filename = Outpath + publicname_1 + '.png' kml_filename = Outpath + publicname_1 + '.kml' kmz_filename = Outpath + publicname_1 + '.kmz' png_filename1 = publicname_1 + '.png' tiff_filename = Outpath + publicname_1 + '.tif' pngFlag = FILE_TEST(png_filename) if pngFlag eq 1 then begin print,'png file already exists' continue EndIF Else Begin ;-----------The current data does not have its bmp file------------ lst_hdf_read, filename, 'WaterDetection', watermask lst_hdf_attr_read, filename, 'WaterDetection', 'ProjectionMinLongitude', ul_lon lst_hdf_attr_read, filename, 'WaterDetection', 'ProjectionMaxLatitude', ul_lat lst_hdf_attr_read, filename, 'WaterDetection', 'ProjectionMaxLongitude', lr_lon lst_hdf_attr_read, filename, 'WaterDetection', 'ProjectionMinLatitude', lr_lat lst_hdf_attr_read, filename, 'WaterDetection', 'ProjectionResolution', resolution usr_ul_lon = ul_lon usr_ul_lat = ul_lat usr_lr_lon = lr_lon usr_lr_lat = lr_lat dimsize1 = size(watermask) col1 = dimsize1[1] row1 = dimsize1[2] subset_col_begin = fix((usr_ul_lon - ul_lon)/resolution + 0.5) subset_row_begin = fix((ul_lat - usr_ul_lat)/resolution + 0.5) subset_col_end = fix((usr_lr_lon - ul_lon)/resolution) subset_row_end = fix((ul_lat - usr_lr_lat)/resolution) if subset_col_begin lt 0 then begin subset_col_begin = 0 endif if subset_row_begin lt 0 then begin subset_row_begin = 0 endif if subset_col_end ge col1 then begin subset_col_end = col1 - 1 endif if subset_row_end ge row1 then begin subset_row_end = row1 - 1 endif subset_watermask = watermask(subset_col_begin:subset_col_end, subset_row_begin: subset_row_end) ;------------------Display data--------------------------------- ;-----------define user definded color table------------------- R = [0, 0, 196, 180, 255, 0, 200, 100, 0, 50, 0, 200,255,255,255,255,255,255] G = [0, 0, 162, 0, 255, 255, 200, 100, 0, 255,255, 255,255,255,200,150,100, 0] B = [0,100, 114, 230, 255, 255, 200, 100, 255, 100, 0, 0,150, 0, 0, 50, 0, 0] num_color = n_elements(R) dimsize = size(subset_watermask) width = dimsize[1] height = dimsize[2] indx = where(subset_watermask eq 88 or subset_watermask eq 38, count) IF count gt 0 then begin subset_watermask[indx] = 38 EndIF indx = where(subset_watermask gt 100, count) IF count gt 0 then begin subset_watermask[indx] = subset_watermask[indx] - 100 EndIF ;-------------density slice------------------ count = 0L FullImg1 = bytarr(width,height) index = bytarr(width,height) indx = where(subset_watermask le 1,count) IF count gt 0 then begin FullImg1[indx] = 0 EndIF index = bytarr(width,height) indx = where(subset_watermask eq 15 ,count) IF count gt 0 then begin FullImg1[indx] = 1 EndIF indx = where(subset_watermask ge 16 and subset_watermask le 17,count) IF count gt 0 then begin FullImg1[indx] = 2 EndIF indx = where(subset_watermask eq 38, count) IF count gt 0 then begin FullImg1[indx] = 3 EndIF indx = where(subset_watermask eq 20,count) IF count gt 0 then begin FullImg1[indx] = 4 EndIF indx = where(subset_watermask eq 27,count) IF count gt 0 then begin FullImg1[indx] = 5 EndIF indx = where(subset_watermask eq 30,count) IF count gt 0 then begin FullImg1[indx] = 6 EndIF indx = where(subset_watermask eq 50 or subset_watermask eq 150,count) IF count gt 0 then begin FullImg1[indx] = 7 EndIF indx = where(subset_watermask eq 99 or subset_watermask eq 199,count) IF count gt 0 then begin FullImg1[indx] = 8 EndIF indx = where( (subset_watermask gt 100 and subset_watermask le 120) or (subset_watermask gt 200 and subset_watermask le 220) ,count) IF count gt 0 then begin FullImg1[indx] = 9 EndIF for i = 10, num_color - 1 do begin indx = where( (subset_watermask le (120 + (i - 9) * 10) $ and subset_watermask gt (120 + (i - 10) * 10)) or (subset_watermask le (220 + (i - 9) * 10) and subset_watermask gt (220 + (i - 10) * 10)) ,count) IF count gt 0 then begin FullImg1[indx] = i EndIF endfor blackindx = where((R eq 0 and G eq 0 and B eq 0), count) write_png, png_filename, FullImg1, R, G, B, $ TRANSPARENT=FullImg1[blackindx],/order ; FullImg1[*,*] = (FullImg1[*,*] + 1)* 14 g_tags = { ModelPixelScaleTag: [ resolution, resolution, 0 ], $ ModelTiepointTag: [0, 0, 0,usr_ul_lon,usr_ul_lat, 0 ], $ GTModelTypeGeoKey: 2s, $ ; (ModelTypeProjected) GTRasterTypeGeoKey: 1s, $ ; (RasterPixelIsArea) GeographicTypeGeoKey: 4326, $ ; (user-defined) GeogGeodeticDatumGeoKey: 6326, $ ; User-Defined GeogLinearUnitsGeoKey: 9001s, $ ; Linear_Meter GeogAngularUnitsGeoKey: 9102s, $ ; Angular_Degree GeogSemiMajorAxisGeoKey: 6378273.0d, $ ; GeogSemiMinorAxisGeoKey:6356889.449d} write_tiff, tiff_filename, subset_watermask, BITS_PER_SAMPLE = 8, GEOTIFF = g_tags NAME = png_filename1 DESC = 'VIIRS/ABI/AHI flood detection map' WRITE_KML, png_filename, NAME, DESC, usr_ul_lat, usr_ul_lon, usr_lr_lat, usr_lr_lon, kml_filename, kmz_filename EndElse endfor end pro lst_hdf_read, filename, SD_name, data sd_id=HDF_SD_START(filename,/READ) if sd_id eq -1 then begin print, 'Fail to open hdf file' return endif index = HDF_SD_NAMETOINDEX(sd_id, SD_name) print, 'Index:' print, index SDS_id = HDF_SD_SELECT(sd_id, HDF_SD_NAMETOINDEX(sd_id, SD_name)) if sd_id eq -1 then begin print, 'Fail to open hdf file' return endif HDF_SD_GETDATA, SDS_id, data ,COUNT=count, START=START, STRIDE= stride HDF_SD_ENDACCESS,sds_id HDF_SD_END,sd_id end pro lst_hdf_attr_read, filename, SD_name, ATRR_name, data sd_id=HDF_SD_START(filename,/READ) if sd_id eq -1 then begin print, 'Fail to open hdf file' return endif SDS_id = HDF_SD_SELECT(sd_id, HDF_SD_NAMETOINDEX(sd_id, SD_name)) if sd_id eq -1 then begin print, 'Fail to open hdf file' return endif attr_index = HDF_SD_ATTRFIND(sds_id, ATRR_name) HDF_SD_ATTRINFO,SDS_id,attr_index, DATA = data, $ NAME=name, TYPE=type, COUNT=count HDF_SD_ENDACCESS,sds_id HDF_SD_END,sd_id end pro lst_hdf_dimsize_read, filename, SD_name, dataType, dims, ndims sd_id=HDF_SD_START(filename,/READ) if sd_id eq -1 then begin print, 'Fail to open hdf file' return endif SDS_id = HDF_SD_SELECT(sd_id, HDF_SD_NAMETOINDEX(sd_id, SD_name)) if sd_id eq -1 then begin print, 'Fail to open hdf file' return endif HDF_SD_GETINFO, SDS_id,DIMS = dims, HDF_TYPE = dataType, NDIMS = ndims HDF_SD_ENDACCESS,sds_id HDF_SD_END,sd_id end PRO WRITE_KML, IMAGE_FILE, NAME, DESC, UL_LAT, UL_LON, LR_LAT, LR_LON, KML_FILE, KMZ_FILE ;- Write a KML file for Google Earth ;- Create KML file openw, lun, kml_file, /get_lun printf, lun, '' printf, lun, '' printf, lun, '' printf, lun, name, format='("", a, "")' printf, lun, desc, format='("", a, "")' printf, lun, '' printf, lun, name, format='("", a, "")' printf, lun, '' printf, lun, '' printf, lun, ul_lat, format='("", f10.4, "")' printf, lun, ul_lon, format='("", f10.4, "")' printf, lun, lr_lat, format='("", f10.4, "")' printf, lun, lr_lon, format='("", f10.4, "")' printf, lun, '0.0' printf, lun, '' printf, lun, '10' printf, lun, '' printf, lun, '' free_lun, lun ;- Create KMZ file ;command = 'zip -r -D -j' + ' ' + kmz_file + ' ' + kml_file + ' ' + image_file ;spawn, command END