利用ENVI_SUBSET_VIA_ROI_DOIT函数实现矢量文件裁剪栅格数据。
测试数据：栅格文件是Mercator投影，矢量文件是Geographic坐标系
1 infile='F:Temp_DataFY3_BEIJING_HFII_DAY_200808031211_Mercator.tif' 2 outfile='F:Temp_DataFY3_BEIJING_HFII_DAY_200808031211_Mercator_Sub.dat' 3 shpfile='F:Shapefileeijing.shp' 4 5 envi_open_file,infile,r_fid=tfid 6 envi_file_query,tfid,dims=dims,ns=ns,nl=nl,nb=nb 7 oproj=envi_get_projection(fid=tfid) 8 iproj=envi_proj_create(/geographic) 9 10 oshp=obj_new('IDLffShape',shpfile) 11 oshp->GetProperty,n_Entities=n_ent,attribute_info=attr_info,n_attributes=n_attr,Entity_type=ent_type 12 roi_ids=lonarr(n_ent) 13 for i=0,n_ent-1 do begin 14 entity=oshp->GetEntity(i) 15 if ptr_valid(entity.parts) ne 0 then begin 16 tempLon=reform((*entity.vertices)[0,*]) 17 tempLat=reform((*entity.vertices)[1,*]) 18 envi_convert_projection_coordinates, $ 19 tempLon,tempLat,iproj, $ 20 xmap,ymap,oproj 21 envi_convert_file_coordinates, $ 22 tfid,xf,yf,xmap,ymap 23 roi_ids[i]=envi_create_roi(ns=ns,nl=nl) 24 envi_define_roi,roi_ids[i],xpts=xf,ypts=yf,/polygon 25 endif 26 27 if i eq 0 then begin 28 xmin=round(min(xf,max=xmax)) 29 ymin=round(min(yf,max=ymax)) 30 endif else begin 31 xmin=xmin<round(min(xf)) 32 xmax=xmax>round(max(xf)) 33 ymin=ymin<round(min(yf)) 34 ymax=ymax>round(max(yf)) 35 endelse 36 oshp->DestroyEntity,entity 37 endfor 38 obj_destroy,oshp 39 xmin=xmin>0 40 xmax=xmax<ns 41 ymin=ymin>0 42 ymax=ymax<nl 43 44 dims=[-1,xmin,xmax,ymin,ymax] 45 envi_doit,'envi_subset_via_roi_doit', $ 46 fid=tfid, $ 47 dims=dims, $ 48 ns=ns,nl=nl, $ 49 pos=indgen(nb), $ 50 background=0, $ 51 roi_ids=roi_ids, $ 52 proj=oproj, $ 53 out_name=outfile 54 envi_file_mng,id=tfid,/remove 55 envi_file_mng,id=rfid,/remove 56 envi_delete_rois,/all
结果如下：
