Thank you James and Devin for your excellent suggestions and
comments!!! I also made it without any errors. I think the problem was
that I used too may envi_file_query statements.
There was no difference between using integer and float for albedo.
However, as devin suggested, MOD43B3 albedo is stored as integer. I
have no reason to use floating point for this... To avoid other
unexpected errors, I need to use integer for albedo from now on.
Here is my final code.
-----------------------------------
PRO MOD43B3_MakeImage_1km_final
WorkDir = 'd:\MODIS_ALL\'
WorkDirOut = 'd:\MODIS43\Processed\'
batch_st = strcompress(WorkDir + 'batch_MOD43.txt', /remove_all)
WorkDirSat = 'D:\MODIS43\MOD43\'
StrMOD = 'MOD43'
albedo = intarr(2, 10, 1200, 1200)
albedo_black = intarr(1200,1200)
albedo_white = intarr(1200,1200)
qc = ulonarr(2, 1200, 1200)
qc1 = ulonarr(1200, 1200)
qc2 = ulonarr(1200, 1200)
OpenR, lun, batch_st, /Get_Lun
numdates = file_lines(batch_st)
Dates = StrArr(numdates)
; Read input dates from batch file
OpenR, lun, batch_st, /Get_lun
ReadF, lun, Dates
Free_Lun, lun, /force
close, /all
FOR j = 0L, numDates-1 DO BEGIN
;- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
- - - -
print, "Now processing MOD43 data from date: ", Dates[j], '
File ', j+1, ' out of ', numDates, $
' overpassing dates. '
Filename = WorkdirSat+Dates[j]
FileID = HDF_OPEN(filename, /read)
sdFileID = HDF_SD_Start(filename, /read) ;The returned
value of this function is the SD ID of the HDF file
sdsID_albedo = HDF_SD_Select(sdFileID, 0) ; Albedo
sdsID_qc = HDF_SD_Select(sdFileID, 1) ; QC
hdf_sd_getdata, sdsID_albedo, albedo
hdf_sd_getdata, sdsID_qc, qc
albedo_black[*,*] = albedo[0, 9, *, *]
albedo_white[*,*] = albedo[1, 9, *, *]
qc1[*,*] = qc[0,*,*]
qc2[*,*] = qc[1,*,*]
;VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV VVVVV
;V V V V V V
;V V V 1. Input map Information(SIN projection) V V V
;V V V 2. Convert SIN to TM Korea V V V
;V V V V V V
;VVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVVV VVVVV
; 1. Input map information (SIN)
mc=[0.5, 0.5, 11119968.509D, 4447338.766D]
ps=[926.6254331D, 926.6254331D]
units =envi_translate_projection_units('Meters')
params1=[6371007.181D, 0, 0, 0]
Projection_Name1 = 'SIN_MODIS'
map_info=ENVI_MAP_INFO_CREATE(type=16, name=Projection_Name1,
params=params1, $
UNITS = units, MC = mc, PS = ps)
;----------------------------------------------------------- ---
; albedo black and white
ENVI_WRITE_ENVI_FILE, albedo_black,
out_name='albedo_black_map.img', map_info=map_info, out_dt=4,$
r_fid=albedo_black_map, sensor_type=32, interleave = 0
ENVI_WRITE_ENVI_FILE, albedo_white,
out_name='albedo_white_map.img', map_info=map_info, out_dt=4,$
r_fid=albedo_white_map, sensor_type=32, interleave = 0
envi_file_query, albedo_black_map, ns=ns, nl=nl, nb=nb
dims = [-1, 0, ns-1, 0, nl-1]
pos_albedo_black = lindgen(nb)
;----------------------------------------------------------- ---
; qc1 and qc2
ENVI_WRITE_ENVI_FILE, qc1 , out_name='qc1_map.img',
map_info=map_info, out_dt=2, $
r_fid=qc1_map, sensor_type=32
ENVI_WRITE_ENVI_FILE, qc2 , out_name='qc2_map.img',
map_info=map_info, out_dt=2, $
r_fid=qc2_map, sensor_type=32
envi_file_query, qc1_map, ns=ns, nl=nl, nb=nb
dims = [-1, 0, ns-1, 0, nl-1]
pos_qc1 = lindgen(nb)
; 2. Convert SIN to TM Korea
; Names of MODIS products
; MOD11A1.A2002001.h28v05.004.2003182143054.hdf
; MOD07_L2.A2003148.1150.004.2003173192335.hdf
; MOD43B3.A2002177.h28v05.004.2003246195929.hdf
StrDate = STRMID(Dates[j], 9, 7)
StrSat = STRMID(Dates[j], 0, 7)
out_name_Albedo_black = WorkDirOut+StrSat+'\'+StrDate+'\'+StrDate
+'_albedo_black_1km'+'.img'
out_name_Albedo_white = WorkDirOut+StrSat+'\'+StrDate+'\'+StrDate
+'_albedo_white_1km'+'.img'
out_name_QC1 = WorkDirOut+StrSat+'\'+StrDate+'\'+StrDate
+'_QC1_1km'+'.img'
out_name_QC2 = WorkDirOut+StrSat+'\'+StrDate+'\'+StrDate
+'_QC2_1km'+'.img'
;----------------------------------------------------------- ---
; Converting to Korea TM
DATUM = 'Tokyo mean'
Projection_Name2= 'Korea - TM (Middle)'
Params2 = [6377397.2, 6356079.0, 38.000000D, 127.002890D,
200000.0, 500000.0, 1.000000]
OUT_Proj = ENVI_PROJ_CREATE(type=3, name=Projection_Name2,
datum=Datum, params=Params2)
;----------------------------------------------------------- ---
; albedo black and white
envi_convert_file_map_projection, fid=albedo_black_map,
pos=pos_albedo_black, dims=dims, o_proj=OUT_Proj, $
o_pixel_size=[1000, 1000],out_name=out_name_albedo_black,
warp_method=2, r_fid=albedo_black_TM, $
resampling=0, background=0
envi_convert_file_map_projection, fid=albedo_white_map,
pos=pos_albedo_black, dims=dims, o_proj=OUT_Proj, $
o_pixel_size=[1000, 1000],out_name=out_name_albedo_white,
warp_method=2, r_fid=albedo_white_TM, $
resampling=0, background=0
;----------------------------------------------------------- ---
; qc1 and qc2
envi_convert_file_map_projection, fid=qc1_map, pos=pos_qc1,
dims=dims, o_proj=OUT_Proj, $
o_pixel_size=[1000, 1000],out_name=out_name_QC1, warp_method=2,
r_fid=qc1_TM, $
resampling=0, background=0 ; NOTICE fid=qc1_map!
envi_convert_file_map_projection, fid=qc2_map, pos=pos_qc1,
dims=dims, o_proj=OUT_Proj, $
o_pixel_size=[1000, 1000],out_name=out_name_qc2, warp_method=2,
r_fid=qc2_TM, $
resampling=0, background=0 ; NOTICE fid=qc2_map!
Close, /all, /force
;------------------------------------------------------
; Cleaning memory
;envi_file_mng, id= albedo_black_map, /remove
;envi_file_mng, id= albedo_white_map, /remove
;envi_file_mng, id= albedo_qc1_map, /remove
;envi_file_mng, id= albedo_qc2_map, /remove
;------------------------------------------------------
; Done with SDS, close the interface
HDF_SD_ENDACCESS, SDSID_albedo
HDF_SD_ENDACCESS, SDSID_qc
HDF_SD_END, sdFileID
HDF_Close, FileID
ENDFOR
print, "Look's O.K!!!"
Free_Lun, lun, /force
END
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
