IDL实现矢量文件裁剪栅格数据

利用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

结果如下：