Hi, all. I would like to add the coastline on my solar zenith angle
(SOLZA) image layers on Korean Peninsula retrieved from MODIS 07
atmospheric profile product. These images were georeferenced on the
basis of the information in the geolocation fields (data set # 0 and
1, latitude and longitude, respectively) of MODIS 07 HDF file. In
order to add the high resolution coastline, I downloaded
Map_GSHHS_shoreline (David Fanning, 2008).
However, at least in my understanding so far, Map_GSHHS_shoreline
seems to draw a coastline on a blank window, and I have no idea how to
link this procedure to the classical MODIS image-processing method.
Please take a look at my sourcecode and give me any suggestions,
especially section 3.2.
I know that someone in this forum may not like the idea to post a full
source code here. I am sorry if I am bothering him, but I hope this
humble source code can help others to begin MODIS data processing in
near future.
Thank you very much in dvance.
----------------------------------------
PRO MOD07_MakeImage_solza_062310 ; Ta, Td, Solza only
;*********************************************************** **************
;******** 1. Preprocessing.
*********
;*********************************************************** **************
; 1.1. Designate paths of workspace and output image folders.
WorkSp07 = 'E:\MODIS07\MOD07_Workspace\'
WorkSpOut07 = 'E:\MODIS07\Image_output\'
batch_st07 = strcompress(WorkSp07+'List_MOD07.txt',/remove_all)
WorkSpHDF = 'E:\MODIS07\MOD07\'
numDates = file_lines(batch_st07)
Dates = StrArr(numDates)
; 1.2. Read input dates from batch file
OpenR, lun, batch_st07, /Get_lun
ReadF, lun, Dates
Free_Lun, lun, /force
FOR j = 0, numDates-1 DO BEGIN
;FOR j = 0, 50 DO BEGIN
s_time = systime(1)
print, "Now processing MOD07 data from date: ", Dates[j], '
File ', j+1, ' out of ', numDates
; 1.3. Open HDF files.
HDFname = WorkSpHDF+Dates[j]
FileOpen = HDF_OPEN(HDFname, /Read)
sdFileID = HDF_SD_Start(HDFname, /Read)
sdsID_lat = HDF_SD_Select(sdFileID, 0)
sdsID_lon = HDF_SD_Select(sdFileID, 1)
sdsID_SOLZA = HDF_SD_Select(sdFileID, 3)
hdf_sd_getdata, sdsID_lat, lat
hdf_sd_getdata, sdsID_lon, lon
hdf_sd_getdata, sdsID_lon, SOLZA_temp
ENVI_WRITE_ENVI_FILE, lat, out_name='lat_temp.img'
ENVI_WRITE_ENVI_FILE, lon, out_name='lon_temp.img'
envi_open_file, 'lat_temp.img', r_fid=lat_fid
envi_open_file, 'lon_temp.img', r_fid=lon_fid
envi_file_query, lat_fid, ns=ns, nl=nl, nb=nb
dims = [-1, 0, ns-1, 0, nl-1]
pos = lindgen(nb)
;*********************************************************** **************
;** 2. CONSTRUCTING GLT FILE FOR GEOREFERENCING **
;*********************************************************** **************
; 2.1. Designate datum and projection name.
DATUM = 'Tokyo mean'
Projection_Name = 'Korea - TM (Middle)'
Params = [6377397.2, 6356079.0, 38.000000, 127.002890, 200000.0,
500000.0, 1.000000]
IN_proj = ENVI_PROJ_CREATE(/Geographic)
OUT_Proj = ENVI_PROJ_CREATE(type=3, name=Projection_Name,
datum=Datum, params=Params)
;2.2. Building GLT file for georeferencing
ENVI_DOIT, 'ENVI_GLT_DOIT', DIMS=dims, I_PROJ=IN_proj,
O_PROJ=OUT_proj, R_FID=GLT_fid,$
ROTATION=0, X_FID=lon_fid,
X_POS=pos[0], Y_FID=lat_fid, Y_POS=pos[0], out_name='GLT.img'
;2.3. Scale HDF data to physical values
SOLZA = CONV_PHYSICAL(sdsID_SOLZA, SOLZA_temp) ; SOLAR zenith
angle [degrees]
;*********************************************************** **************
;** 3. Add Korea TM coordinates on georeferenced files *********
;*********************************************************** **************
;3.1. Printout these arrays on image files.
StrDate = STRMID(Dates[j], 10,7)
StrFile = STRMID(Dates[j], 10, 12)
StrTime = STRMID(Dates[j], 18, 4)
StrSat = STRMID(Dates[j], 0, 5)
ENVI_WRITE_ENVI_FILE, SOLZA, out_name='SOLZA_temp.img'
envi_open_file, 'SOLZA_temp.img', r_fid=SOLZA_fid
envi_file_query, SOLZA_fid, ns=ns, nl=nl, nb=nb
pos = lindgen(nb)
out_name_SOLZA = WorkSpOut07+StrSat+'/'+StrDate+'/'+StrSat+StrFile
+'_SOLZA0_5km'+'.img'
ENVI_DOIT, 'ENVI_GEOREF_FROM_GLT_DOIT', FID=SOLZA_fid,
GLT_FID=GLT_fid, POS=pos, R_fid=SOLZA_G_fid, out_name=out_name_SOLZA
; 3.2. Add coast lines for Korean penninsular. - not yet
completed...!!!
Set_plot, 'z'
datafile = 'C:\Program Files\ITT\IDL708\lib
\GSHHS_2.0\gshhs_h.b'
Device, set_resolution = [500, 350]
pos = [0.1,0.1, 0.9, 0.8]
Map_Set, 37.0, 127.0, Position=pos, Scale=15e6, /Mercator, /
NoBorder
Polyfill, [pos[0], pos[0], pos[2], pos[2], pos[0]], $
[pos[1], pos[3], pos[3], pos[1], pos[1]], $
/Normal, Color=FSC_Color('Almond')
Map_GSHHS_Shoreline, datafile, Level=3, /Outline
b = tvrd()
write_jpeg, 'test09.jpg', b
;*********************************************************** **************
;****** 4. Finish up the image processing
*************
;*********************************************************** **************
; 4.1. Clean memories.
NextHdf:
ENVI_FILE_MNG, ID=lat_fid, /remove ;, /Delete
ENVI_FILE_MNG, ID=lon_fid, /remove ;, /Delete
ENVI_FILE_MNG, ID=SOLZA_G_fid, /remove
NextHdf_2:
; 4.2. Done with SDS, close the interface
HDF_SD_ENDACCESS,SdsID_LAT
HDF_SD_ENDACCESS,SdsID_LON
HDF_SD_ENDACCESS,SdsID_SOLZA
HDF_SD_END,sdFileID
HDF_CLOSE,FileOpen
Close, /all, /force
ENDFOR
e_time = systime(1)
print, "It's done!!! Elapsed time for this procedure: ", e_time -
s_time
close, /all
END
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
