Hi all, I am seriously hoping someone out there can help me with
this . I am trying to georeference MODIS L1B hdf files using envi
batch mode. I found a bit about this on the net (particularly the
Ocean Color forums) and have tried to emulate their method, although
as they wern't using L1B, I think I am encountering some sort of
additional issue.
What I have so far is pasted below. This actually does output a
georeferenced-looking file, however the values are somewhat different
to the result if I georeference it with ENVI normally. I think the
issue is either in making the GLT, or in the georeferencing itself.
(nb. I also tried the MODIS Conversion Toolkit, which georeferences
the file just fine, but as it converts it to a .img is not suitable
for what I need)
Any help greatly appreciated,
Cheers,
Helen
open lon file
ENVI_OPEN_DATA_FILE, filename, r_fid=x_fid, /hdf_sd,
hdfsd_dataset=1, hdfsd_interleave=0
if (x_fid eq -1) then begin
envi_batch_exit
return
endif
ENVI_FILE_QUERY, x_fid, ns=xns, nl=xnl, nb=xnb
x_pos=lindgen(xnb)
xdims=[-1L, 0, xns-1, 0, xnl-1]
;open lat file
ENVI_OPEN_DATA_FILE, filename, r_fid=y_fid, /hdf_sd,
hdfsd_dataset=0, hdfsd_interleave=0
if (y_fid eq -1) then begin
envi_batch_exit
return
endif
ENVI_FILE_QUERY, y_fid, ns=yns, nl=ynl, nb=ynb
y_pos=lindgen(ynb)
ydims=[-1L, 0, yns-1, 0, ynl-1]
;open file
ENVI_OPEN_DATA_FILE, filename, r_fid=therm_fid, /hdf_sd,
hdfsd_dataset=4, hdfsd_interleave=0
if (therm_fid eq -1) then begin
envi_batch_exit
return
endif
ENVI_FILE_QUERY, therm_fid, ns=ns, nl=nl, nb=nb,
data_type=data_type, bnames=bnames
pos=lindgen(nb)
dims=[-1L,0, ns-1, 0, nl-1]
;CD TO THE OUTPUT DIRECTORY
cd, Outputdir
; Figure out what UTM zone we're in.
lat_data=envi_get_data(fid=y_fid, dims=ydims, pos=0)
lon_data=envi_get_data(fid=x_fid, dims=xdims, pos=0)
lat_lon=[lat_data, lon_data]
zone = fix(31.0 + lat_lon[1]/6.0)
south = (lat_lon[0] lt 0)
nlat= n_elements(lat_data)
nlon = n_elements(lon_data)
lat=mean(lat_data) ;find the average (so approx center of image) to
define the UTM zone.
lon=mean(lon_data)
LongTemp = (Lon+180)-fix((Lon+180)/360)*360-180; // -180.00 .. 179.9
ZoneNumber = FIX((LongTemp + 180.0)/6.0) + 1
if Lat GE 56.0 AND Lat LE 64.0 AND LongTemp GE 3.0 AND LongTemp LE
12.0 then ZoneNumber = 32
;special zones for Svalbard
IF lat GE 72.0 and lat LE 84.0 then begin
if longtemp ge 0.0 and longtemp lt 9.0 then zonenumber =32
if longtemp GE 9.0 and longtemp lt 21.0 then zonenumber = 33
if longtemp ge 21.0 and longtemp lt 33.0 then zonenumber = 35
if longtemp ge 33.0 and longtemp lt 42.0 then zonenumber = 37
ENDIF
south = (lat lt 0)
; Make the GLT.
zone=zonenumber
envi_file_query, therm_fid, sname=sname
out_name1='GLT_file_'+fname
;pixel_size=[1000.0, 1000.0]
rotation=0.0
i_proj = envi_proj_create(/geographic)
o_proj = envi_proj_create(/utm, zone=zone, south=south)
envi_glt_doit, i_proj=i_proj, o_proj=o_proj, out_name=out_name1,
r_fid=glt_fid,$
x_fid=x_fid, y_fid=y_fid, x_pos=0, y_pos=0,pixel_size=pixel_size,
rotation=rotation
;therm_fid contains all the emissive bands.We only need bands 29-32
for the SO2 retrieval.
t_fid=
[therm_fid,therm_fid,therm_fid,therm_fid,therm_fid,therm_fid ,therm_fid,therm_fid,therm_fid,therm_fid,
$
therm_fid,therm_fid,therm_fid,therm_fid,therm_fid,therm_fid, therm_fid]
t_pos=[11,10,9,8]
; Georeference the image from the GLT.
out_name2=fname+'georefferenced'
envi_doit, 'envi_georef_from_glt_doit', fid=t_fid, $
glt_fid=glt_fid, out_name=out_name2, pos=t_pos, $
subset=dims, r_fid=r_fid
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
