|
|
| Re: problem with subset one image by another [message #59428 is a reply to message #59335] |
Thu, 27 March 2008 08:56   |
Jean H.
Messages: 472
Registered: July 2006
|
Senior Member |
|
|
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
|
|
|
|
| Re: problem with subset one image by another [message #59430 is a reply to message #59335] |
Thu, 27 March 2008 08:40 |
negra
Messages: 9
Registered: March 2008
|
Junior Member |
|
|
>
> 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
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
|
|
|
|
|
|
|
|
|
|
| Re: problem with subset one image by another [message #59438 is a reply to message #59335] |
Wed, 26 March 2008 19:10 |
negra
Messages: 9
Registered: March 2008
|
Junior Member |
|
|
On 27 мар, 01:52, James Kuyper <jameskuy...@verizon.net> wrote:
> negra wrote:
>> On 22 ÍÁÒ, 22:12, David Fanning <n...@dfanning.com> wrote:
>>> negra writes:
>>> Well, there probably is nothing easier than this in IDL,
>>> either, so I think folks are either waiting for you to
>>> realize this or waiting for you to explain in a little
>>> more detail exactly what it is you are trying to do. :-)
>
>>> Cheers,
>
>>> David
>>> --
>> I have a large number of georeferenced MODIS images in ENVI. And it's
>> just attempt to get rid from leaving for area of interest of
>> territory.
>
> Are these Level 2 images, or Level 3? Level 2 products are organized
> around 5-minute granules of data, and are generally in HDF-EOS Swath
> format. For best results, a Geolocation (MOD03 or MYD03) file should be
> used for georeferencing such data.
>
> Level 3 images are rebinned to a rectangular grid in som particular map
> projection (often sinusoidal), and are stored in HDF-EOS Grid format.
>
> The techniques that should be used are somewhat different in the two cases.
>
> And, as David said, we need a whole lot more details about what it is
> you are trying to do.
It's images processed from Level1 or Level 0 to Level 2 products. Some
I take them from our receiving station, others, when we have problems
with receiving, I take as 5-minute granules of data as Level 0 (pds)
from ftp. Then I process them till level 2 product. And georefernce in
IDL (because I have IDL routine for processing a large number of data)
to geografical projection. And then I calculate indexes what I need.
For example it's Vegetation, Snow, Dust and other calculations.
For different arithmetic procedures in ENVI&IDL I must have same size
images. May be it would be easier if I will subset by coordinates of
four corner?
But some times images smaller than the area of interest(AOI). Because
AOI it's the territory of whole our country. In this case, I have
procedure of filling by null data till needed size. But in future I
plan to use this procedure (If I wrote it :-)) certainly) to smaller
regions.
Gulshat,
P.S. Excuse me, for my English If something wrong.
|
|
|
|
| Re: problem with subset one image by another [message #59441 is a reply to message #59335] |
Wed, 26 March 2008 12:52 |
jameskuyper
Messages: 79
Registered: October 2007
|
Member |
|
|
negra wrote:
> On 22 О©╫О©╫О©╫, 22:12, David Fanning <n...@dfanning.com> wrote:
>> negra writes:
>> Well, there probably is nothing easier than this in IDL,
>> either, so I think folks are either waiting for you to
>> realize this or waiting for you to explain in a little
>> more detail exactly what it is you are trying to do. :-)
>>
>> Cheers,
>>
>> David
>> --
> I have a large number of georeferenced MODIS images in ENVI. And it's
> just attempt to get rid from leaving for area of interest of
> territory.
Are these Level 2 images, or Level 3? Level 2 products are organized
around 5-minute granules of data, and are generally in HDF-EOS Swath
format. For best results, a Geolocation (MOD03 or MYD03) file should be
used for georeferencing such data.
Level 3 images are rebinned to a rectangular grid in som particular map
projection (often sinusoidal), and are stored in HDF-EOS Grid format.
The techniques that should be used are somewhat different in the two cases.
And, as David said, we need a whole lot more details about what it is
you are trying to do.
|
|
|
|
| Re: problem with subset one image by another [message #59486 is a reply to message #59335] |
Fri, 28 March 2008 21:18 |
Wasit.Weather
Messages: 62
Registered: February 2008
|
Member |
|
|
On Mar 28, 8:54 pm, negra <chem...@mail.ru> wrote:
> On 28 мар, 22:22, Jean H <jghas...@DELTHIS.ucalgary.ANDTHIS.ca> wrote:
>
>
>
>
>
>> negra wrote:
>>> On 27 ÍÁÒ, 21:56, Jean H <jghas...@DELTHIS.ucalgary.ANDTHIS.ca> 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 = 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
>
>>> It's write that something wrong but I don't know how will be right.
>>> Where is mistake?
>
>>> pos = lindgen(file_nb)
>>> input = ENVI_GET_DATA(fid=file_fid, dims=file_dims, POS=pos)
>>> ;
>>> subset = input[file1_xf[0]:file1_xf[1], file1_yf[0]:file1_yf[1]]
>
>> This will not work, as file1_xf contains the Cartesian coordinate of the
>> mask.
>> first, do
>
>> envi_convert_file_coordinates, subset_fid, subset_xf, subset_yf,
>> Mask_xmap, Mask_ymap
>
>> so with the above, you have, with respect to the subset image, the
>> Cartesian coordinates of the Mask corners.
>
>> Then do
>> subset = input[subset_xf[0]:subset_xf[1],subset_yf[0]:subset_yf[1]]
>
>>> out_names = 'subset.img'
>>> ENVI_WRITE_ENVI_FILE, subset, DATA_TYPE=4 , OUT_NAME=out_names,
>>> R_FID=s_fid
>
>> ENVI_WRITE_ENVI_FILE requires the other following keywords:
>> NB=integer | NL=integer | NS=integer
>> OFFSET=value
>
>> Jean- Скрыть цитируемый текст -
>
>> - Показать цитируемый текст -
>
> I tried to do as you suggested, but it was unsuccessful.
> May be it must look like this?
> min_x = max(file1_xmap[0], file_xmap[0]); Attempt to store into an
> expression: <DOUBLE ( 63.295452)>.
>
> max_x = min(file1_xmap[1], file_xmap[1])
> min_y = max(file1_ymap[0], file_ymap[0])
> max_y = min(file1_ymap[1], file_ymap[1])
> subset = input[min_x:max_x, min_y:max_y]
> But it give me mistake on the first line.
> Do you know where is mistake?
>
> Gulshat- Hide quoted text -
>
> - Show quoted text -
Hi Gulshat!
Your name is interesting. Something like Uyghurian or Hungarian name.
|
|
|
|
| Re: problem with subset one image by another [message #59487 is a reply to message #59335] |
Fri, 28 March 2008 18:54 |
negra
Messages: 9
Registered: March 2008
|
Junior Member |
|
|
On 28 мар, 22:22, Jean H <jghas...@DELTHIS.ucalgary.ANDTHIS.ca> wrote:
> negra wrote:
>> On 27 ÍÁÒ, 21:56, Jean H <jghas...@DELTHIS.ucalgary.ANDTHIS.ca> 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 = 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
>
>> It's write that something wrong but I don't know how will be right.
>> Where is mistake?
>
>> pos = lindgen(file_nb)
>> input = ENVI_GET_DATA(fid=file_fid, dims=file_dims, POS=pos)
>> ;
>> subset = input[file1_xf[0]:file1_xf[1], file1_yf[0]:file1_yf[1]]
>
> This will not work, as file1_xf contains the Cartesian coordinate of the
> mask.
> first, do
>
> envi_convert_file_coordinates, subset_fid, subset_xf, subset_yf,
> Mask_xmap, Mask_ymap
>
> so with the above, you have, with respect to the subset image, the
> Cartesian coordinates of the Mask corners.
>
> Then do
> subset = input[subset_xf[0]:subset_xf[1],subset_yf[0]:subset_yf[1]]
>
>> out_names = 'subset.img'
>> ENVI_WRITE_ENVI_FILE, subset, DATA_TYPE=4 , OUT_NAME=out_names,
>> R_FID=s_fid
>
> ENVI_WRITE_ENVI_FILE requires the other following keywords:
> NB=integer | NL=integer | NS=integer
> OFFSET=value
>
> Jean- Скрыть цитируемый текст -
>
> - Показать цитируемый текст -
I tried to do as you suggested, but it was unsuccessful.
May be it must look like this?
min_x = max(file1_xmap[0], file_xmap[0]); Attempt to store into an
expression: <DOUBLE ( 63.295452)>.
max_x = min(file1_xmap[1], file_xmap[1])
min_y = max(file1_ymap[0], file_ymap[0])
max_y = min(file1_ymap[1], file_ymap[1])
subset = input[min_x:max_x, min_y:max_y]
But it give me mistake on the first line.
Do you know where is mistake?
Gulshat
|
|
|
|
| Re: problem with subset one image by another [message #59506 is a reply to message #59335] |
Fri, 28 March 2008 09:22 |
Jean H.
Messages: 472
Registered: July 2006
|
Senior Member |
|
|
negra wrote:
> On 27 О©╫О©╫О©╫, 21:56, Jean H <jghas...@DELTHIS.ucalgary.ANDTHIS.ca> 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 = 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
>
> It's write that something wrong but I don't know how will be right.
> Where is mistake?
>
> pos = lindgen(file_nb)
> input = ENVI_GET_DATA(fid=file_fid, dims=file_dims, POS=pos)
> ;
> subset = input[file1_xf[0]:file1_xf[1], file1_yf[0]:file1_yf[1]]
This will not work, as file1_xf contains the Cartesian coordinate of the
mask.
first, do
envi_convert_file_coordinates, subset_fid, subset_xf, subset_yf,
Mask_xmap, Mask_ymap
so with the above, you have, with respect to the subset image, the
Cartesian coordinates of the Mask corners.
Then do
subset = input[subset_xf[0]:subset_xf[1],subset_yf[0]:subset_yf[1]]
> out_names = 'subset.img'
> ENVI_WRITE_ENVI_FILE, subset, DATA_TYPE=4 , OUT_NAME=out_names,
> R_FID=s_fid
ENVI_WRITE_ENVI_FILE requires the other following keywords:
NB=integer | NL=integer | NS=integer
OFFSET=value
Jean
|
|
|
|
| Re: problem with subset one image by another [message #59508 is a reply to message #59428] |
Fri, 28 March 2008 08:30 |
negra
Messages: 9
Registered: March 2008
|
Junior Member |
|
|
On 27 мар, 21:56, Jean H <jghas...@DELTHIS.ucalgary.ANDTHIS.ca> 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 = 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
>
It's write that something wrong but I don't know how will be right.
Where is mistake?
pos = lindgen(file_nb)
input = ENVI_GET_DATA(fid=file_fid, dims=file_dims, POS=pos)
;
subset = input[file1_xf[0]:file1_xf[1], file1_yf[0]:file1_yf[1]]
out_names = 'subset.img'
ENVI_WRITE_ENVI_FILE, subset, DATA_TYPE=4 , OUT_NAME=out_names,
R_FID=s_fid
|
|
|
|
|
|
|
|
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
