On 27 =CD=C1=D2, 21:56, Jean H <jghas...@[EMAIL PROTECTED]
>
wrote:=
> 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 =3D FILE_SEARCH('MOD02HKM.*.img',count=3Dnumfiles)
> > PRINT, '# NDSI files:',N_ELEMENTS(FILE_SEARCH('MOD02HKM.*.img'))
> > print, FILE_SEARCH(HKMfiles)
> > for j=3D0,numfiles-1 DO BEGIN
> > HKM_name =3D 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=3D'batch.txt'
> > ;
> > ;Open the input files
> > ;kz_hkm mask file
> > envi_open_file, 'C:\Scandata\L1A\mask\kz_hkm', r_fid=3Dfile1_fid
> > if (file1_fid eq -1) then begin
> > =9A =9A =9Aenvi_batch_exit
> > =9A =9A =9Areturn
> > endif
> > envi_open_file,'C:\Scandata\L1A\'+ HKMfiles[j], r_fid=3Dfile_fid
> > if (file_fid eq -1) then begin
> > =9A =9A =9Aenvi_batch_exit
> > =9A =9A =9Areturn
> > endif
> > envi_file_query, file1_fid, dims=3Dfile1_dims, ns=3Dfile1_ns,
nl=3Dfile1=
_nl,
> > nb=3Dfile1_nb
> > file1_dims =3D [-1L,0,file1_ns-1,0,file1_nl-1]
> > file1_mapinfo =3D envi_get_map_info(fid=3Dfile1_fid)
> > print, file1_mapinfo
> > file1_xf =3D [0,file1_ns-1]
> > file1_yf =3D [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=3Dfile_dims, ns=3Dfile_ns,
nl=3Dfile_nl,=
> > nb=3Dfile_nb
> > file_dims =3D [-1L,0,file_ns-1,0,file_nl-1]
> > file_mapinfo =3D envi_get_map_info(fid=3Dfile_fid)
> > print, file_mapinfo
> > pos =9A=3D lindgen(file_nb)
> > out_namea =3D HKMfiles[j]+'subset.img'
> > file_mapinfo =3D envi_get_map_info(fid=3Dfile_fid)
> > _xfactor =3D file1_mapinfo.ps[0]/file_mapinfo.ps[0]
> > _yfactor =3D file1_mapinfo.ps[1]/file_mapinfo.ps[1]
> > print, [_xfactor, _yfactor]
>
> > file_xf =3D [0,file_ns-1]
> > file_yf =3D [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. =9Ado subset =3D
> 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
>
It's write that something wrong but I don't know how will be right.
Where is mistake?
pos =3D lindgen(file_nb)
input =3D ENVI_GET_DATA(fid=3Dfile_fid, dims=3Dfile_dims, POS=3Dpos)
;
subset =3D input[file1_xf[0]:file1_xf[1], file1_yf[0]:file1_yf[1]]
out_names =3D 'subset.img'
ENVI_WRITE_ENVI_FILE, subset, DATA_TYPE=3D4 , OUT_NAME=3Dout_names,
R_FID=3Ds_fid
|
