|
Hi Angelica
I'm sending my script
PRO mod13a3_reproject
COMPILE_OPT STRICTARR
; ****************************************************************************************
; DESCRIPTION: Program to Reproject MODIS STACK subset based on ROI
; *********************************************************************
; AUTHOR: David Lopez A & Stef Lhermitte
; *********************************************************************
; DATE: 2010/12/06
; *********************************************************************
; FORMATS
; *********************************************************************
F80='(80("-"))'
F80E='(80("="))'
F120='(120("-"))'
F120E='(120("="))'
F120S='(120("*"))'
; *********************************************************************
; PARAMETERS
; *********************************************************************
ifolder = '/media/ION300/rawdata/modis13a3/stack/'
tfolder = '/media/ION300/rawdata/modis13a3/stack/'
ofolder = '/media/ION300/rawdata/modis13a3/stack/subset/'
vfolder = '/media/ION300/rawdata/modis13a3/stack/'
img_fname = ifolder+'MOD13A3.B0.stack'
evf_fname = vfolder+'borde_.evf'
; *********************************************************************
PRINT, FORMAT=F120S
if (FILE_TEST(ifolder, /DIRECTORY) ne 1) then begin
PRINT, 'Input folder does not exist: ', ifolder
PRINT, FORMAT=F120S
return
endif
if (FILE_TEST(tfolder, /DIRECTORY) ne 1) then begin
PRINT, 'Working folder does not exist: ', tfolder
PRINT, FORMAT=F120S
return
endif
if (FILE_TEST(ofolder, /DIRECTORY) ne 1) then begin
PRINT, 'Output folder does not exist: ', ofolder
PRINT, FORMAT=F120S
return
endif
if (FILE_TEST(vfolder, /DIRECTORY) ne 1) then begin
PRINT, 'Vector folder does not exist: ', vfolder
PRINT, FORMAT=F120S
return
endif
; Open image file to reproject
ENVI_OPEN_FILE, img_fname, r_fid=fid
if (fid[0] eq -1) then begin
PRINT, 'Input image does not exist: '+ im_fname
PRINT, FORMAT=F120S
return
endif
PRINT, 'Opening input image :'+ img_fname
ENVI_FILE_QUERY, fid, bnames=bnames, ns=ns, nl=nl, nb=nb
i_proj = ENVI_GET_PROJECTION(fid=fid)
map_info = ENVI_GET_MAP_INFO(fid=fid)
; Open vector file with ROI
PRINT, 'Opening input vector :'+ evf_fname
evf_exist = FILE_TEST(evf_fname)
if (evf_exist ne 1) then begin
PRINT, 'Input vector does not exist: '+ evf_fname
PRINT, FORMAT=F120S
return
endif
evf_id = ENVI_EVF_OPEN(evf_fname)
; Get map_info from vector file
dims = [-1,0,ns-1,0,nl-1]
ENVI_EVF_INFO, evf_id, projection=o_proj, num_recs=num_recs
PRINT, FORMAT=F80
PRINT, 'Number of Records: ',num_recs
for i=0,num_recs-1 do begin
xy = envi_evf_read_record(evf_id, i, parts_ptr=parts_ptr)
PRINT, 'Number of nodes in Record ' + strtrim(i+1,2) + ': ', n_elements(xy[0,*])
PRINT, 'Multiparts polygon ptr: ' & PRINT, STRING(parts_ptr)
ixmap=TRANSPOSE(xy[0,parts_ptr[0]:(parts_ptr[1]-1)])
iymap=TRANSPOSE(xy[1,parts_ptr[0]:(parts_ptr[1]-1)])
ENVI_CONVERT_PROJECTION_COORDINATES, ixmap, iymap, o_proj, oxmap, oymap, i_proj
ENVI_CONVERT_FILE_COORDINATES, fid, xf, yf, oxmap, oymap
dims=[-1,max([min(FLOOR(xf)),dims[1]]),min([max(CEIL(xf)),dims[2]]),max([min(FLOOR(yf)),dims[3]]),min([max(CEIL(yf)),dims[4]])]
endfor
PRINT,' Dims for reprojecting : ' + STRJOIN(STRTRIM(STRING(dims),2), ' - ')
pos = INDGEN(nb)
; Reproject
ENVI_CONVERT_FILE_MAP_PROJECTION, fid=fid, pos=pos, dims=dims, o_proj=o_proj, $
out_name=img_fname+'.reproj' , out_bname=bnames, warp_method=2, resampling=0, background=0, $
r_fid=fid_repr
; Get info for reprojected file
ENVI_FILE_QUERY, fid_repr, bnames=bnames, ns=ns, nl=nl, nb=nb, dims=dims
; Define ROI
roi_id = ENVI_CREATE_ROI(name='Study area', ns=ns, nl=nl)
; Get vector info from vector file
PRINT, FORMAT=F80
PRINT, 'Number of Records: ',num_recs
for i=0,num_recs-1 do begin
xy = envi_evf_read_record(evf_id, i)
PRINT, 'Number of nodes in Record ' + strtrim(i+1,2) + ': ', n_elements(xy[0,*])
PRINT, 'Multiparts polygon ptr: ' & PRINT, STRING(parts_ptr)
xmap=TRANSPOSE(xy[0,parts_ptr[0]:(parts_ptr[1]-1)])
ymap=TRANSPOSE(xy[1,parts_ptr[0]:(parts_ptr[1]-1)])
ENVI_CONVERT_FILE_COORDINATES, fid_repr, xf, yf, xmap, ymap
ENVI_DEFINE_ROI, roi_id, /polygon, xpts=xf, ypts=yf
endfor
; Create mask file based on ROI
ENVI_DOIT, 'envi_roi_to_image_doit', fid=fid_repr, roi_ids=roi_id, $
out_name='StudyAreaMask', class_values=[1], r_fid=temp_mask
ENVI_FILE_MNG, id=temp_mask, /REMOVE
ENVI_SETUP_HEAD, fname='StudyAreaMask', ns=ns, nl=nl, nb=1, $
interleave=0, data_type=1, $
offset=0, map_info=map_info, $
bnames='Mask Band', r_fid=mask_id, $
/write, /open
print, 'Dims before :' +dims
; Set output dims
dims=[-1,max([min(FLOOR(xf)),0]),min([max(CEIL(xf)),ns-1]),max([min(FLOOR(yf)),0]),min([max(CEIL(yf)),nl-1])]
print, 'Dims after :' +dims
; Apply mask and subset
ENVI_MASK_APPLY_DOIT, FID = fid_repr, POS = pos, DIMS = dims, M_FID = mask_id, M_POS = 0, VALUE = 0, $
out_name=img_fname+'.masked' , out_bname=bnames, r_fid=masked_id
PRINT, FORMAT=F120S
END
Regards
|