I'd recommend that you use the MODIS Conversion Toolkit (MCTK). It is
an ENVI plugin that was specifically designed for problems like this.
You can download it from the ITT VIS website:
http://www.ittvis.com/UserCommunity/CodeLibrary.aspx
Search for "MODIS" and it will be one of the first entries in the
returned list. I received your email about this, too, and replied
with the above recommendation.
On Sep 10, 8:05 am, hethomas <het...@googlemail.com> wrote:
> 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
|