negra wrote:
>> Without more details about what you want to do it is difficult,
>> bordering upon impossible, for us to give you any useful advice about
>> how to do it.
>
> here is idl routine what I written. The beginning of it, is working.
> pro spat_subset
> cd, 'C:\Scandata\L1a'
> HKMfiles = FILE_SEARCH('MOD02HKM.*.img',count=numfiles)
> PRINT, '# NDSI files:',N_ELEMENTS(FILE_SEARCH('MOD02HKM.*.img'))
> print, FILE_SEARCH(HKMfiles)
> for j=0,numfiles-1 DO BEGIN
> HKM_name = HKMfiles[j]
> print,HKM_name
> ;first restore all base save files
> ;
> envi, /restore_base_save_files
> ;
> ;Initialize ENVI and send all errors and warnings to the file
> batch.txt
> ;
> envi_batch_init, log_file='batch.txt'
> ;
> ;Open the input files
> ;kz_hkm mask file
> envi_open_file, 'C:\Scandata\L1A\mask\kz_hkm', r_fid=file1_fid
> if (file1_fid eq -1) then begin
> envi_batch_exit
> return
> endif
> envi_open_file,'C:\Scandata\L1A\'+ HKMfiles[j], r_fid=file_fid
> if (file_fid eq -1) then begin
> envi_batch_exit
> return
> endif
> envi_file_query, file1_fid, dims=file1_dims, ns=file1_ns, nl=file1_nl,
> nb=file1_nb
> file1_dims = [-1L,0,file1_ns-1,0,file1_nl-1]
> file1_mapinfo = envi_get_map_info(fid=file1_fid)
> print, file1_mapinfo
> file1_xf = [0,file1_ns-1]
> file1_yf = [0,file1_nl-1]
> envi_convert_file_coordinates, file1_fid, file1_xf, file1_yf,
> file1_xmap, file1_ymap, /to_map
> print, 'UL corner:',file1_xmap[0],file1_ymap[0]
> print, 'LR corner:',file1_xmap[1],file1_ymap[1]
> ;Longitude 44.39011111 - 88.38219444
> ;Latitude 36.20264722 - 56.35832500
> ;subset #1
>
> envi_file_query, file_fid, dims=file_dims, ns=file_ns, nl=file_nl,
> nb=file_nb
> file_dims = [-1L,0,file_ns-1,0,file_nl-1]
> file_mapinfo = envi_get_map_info(fid=file_fid)
> print, file_mapinfo
> pos = lindgen(file_nb)
> out_namea = HKMfiles[j]+'subset.img'
> file_mapinfo = envi_get_map_info(fid=file_fid)
> _xfactor = file1_mapinfo.ps[0]/file_mapinfo.ps[0]
> _yfactor = file1_mapinfo.ps[1]/file_mapinfo.ps[1]
> print, [_xfactor, _yfactor]
>
> file_xf = [0,file_ns-1]
> file_yf = [0,file_nl-1]
> envi_convert_file_coordinates, file_fid, file_xf, file_yf, file_xmap,
> file_ymap, /to_map
you have to do the opposite: take the projected coordinates of the mask
image, and convert it to the Cartesian coordinate of the image to subset.
then, don't do the "resize", as you are just changing the pixel size and
not the covered area. do subset =
image[maskXmin:maskXmax,maskYmin:maskYmax]
Jean
> print, 'UL corner:',file_xmap[0],file_ymap[0]
> print, 'LR corner:',file_xmap[1],file_ymap[1]
> ;Longitude 63.29925556 - 106.42788056
> ;Latitude 42.01251944 - 64.14574167
> ; I have coordinates of image corners
>
> ;--------------------------------------
> ;here I think must have somthing like this, but I'm not sure
> if file_xmap[0]< file1_xmap[0]then file1_xmap[0] else file_xmap[0]
> if file_xmap[1]> file1_xmap[1]then file1_xmap[1] else file_xmap[1]
> if file_ymap[0]< file1_ymap[0]then file1_ymap[0] else file_ymap[0]
> if file_ymap[1]> file1_ymap[1]then file1_ymap[1] else file_ymap[1]
> then I think must been calculation number of pixels
>
> ;-----------------------------------------------------
> envi_doit, 'resize_doit', $
> fid=file_fid, pos=pos, dims=file_dims, interp=0, rfact=[_xfactor,
> _yfactor], $
> out_name=out_namea, r_fid=file2_fid
> endfor
> end
|
