Messages from Usenet group comp.lang.idl-pvwave, compiled by Paulo Penteado
problem with subset one image by another
Hello all! Could anybody help me with spatial subsetting one image by another in idl. It's easy in Envi, but I must process a large number of files. Thanks! Gulshat *** Sent From/Enviado desde: http://groups.expo.st ***
Re: problem with subset one image by another
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.
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.
negra wrote: ... > 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. OK, that makes it a little bit clearer what you are doing. What we need is more details about what it is that you want to do, and haven't figured out how to do. What have you tried? What's wrong with the approaches you have tried?
(If I wrote it :-)) certainly) to smaller >> regions. > > OK, that makes it a little bit clearer what you are doing. What we need > is more details about what it is that you want to do, and haven't > figured out how to do. What have you tried? What's wrong with the > approaches you have tried? In helps for ENVI&IDL, ittvis and other web sites I don't find something what I could use. I don't know even from what to begin. And without experience in programming it's still difficultly, write this procedure independently.
negra wrote: > (If I wrote it :-)) certainly) to smaller >>> regions. >> OK, that makes it a little bit clearer what you are doing. What we need >> is more details about what it is that you want to do, and haven't >> figured out how to do. What have you tried? What's wrong with the >> approaches you have tried? > > > In helps for ENVI&IDL, ittvis and other web sites I don't find > something what I could use. I don't know even from what to begin. And > without experience in programming it's still difficultly, write this > procedure independently. Don't concentrate on how to do it; that's never the right place to begin. First, concentrate on clearly describing what it is you want to do. What do you start with, and what do you want to have when it's finished? Don't worry about the details of how to program it and what tools to use until you've taken care of the first part. 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.
> > 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
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
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
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
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
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.
> 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)>. look at the help file for "max". You must give it an array of values... the 2nd argument is a variable name that will contain the index of the max value... so try max([var1,var2]) You will want to round the coordinates properly as well! Jean > > Gulshat
negra wrote: ... > 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]) As Jean has already pointed out, you'd have to organize things differently to use MAX() and MIN() here. However, it's far simpler to avoid MAX() and MIN() altogether. max_x = file1_xmap[1] < file_xmap[1] min_y = file1_ymap[0] > file_ymap[0] max_y = file1_ymap[1] < file_ymap[1] It all looks pretty confusing if you're used to a language (there are quite a few) where '<' and '>' are comparison operators, but in IDL the corresponding comparison operators are called 'lt' and 'gt'. The '<' and '>' operators return the minimum and maximum, respectively, of their operands.
I want to say thanks to all who tried to help me. All has turned out.
On Apr 4, 10:32 am, negra <chem...@mail.ru> wrote: > I want to say thanks to all who tried to help me. All has turned out. Can you share the whole PRO...END program here. Or send me by email. Actually, that is what I am trying to work it out now. email: wasit.weather@gmail.com