地學分析與三維可視化實習一

1、地學分析與三維可視化實習報告一姓名：方志策班級：064131學號：20131003574提交日期：2015.11.5作業(yè)一、任務一、利用對象圖形法創(chuàng)建一個三維立方體，并將各頂點設置為不同的顏色。IDL代碼如下：oWindow=OBJ_NEW('IDLgrWindow',dimension=400,400);建立一個圖像窗口，初始設置圖像大小為400*400oView=OBJ_NEW('IDLgrView',viewPlane_Rect=-1,-1,3,3,zClip=3,-3,eye=2) oModel=OBJ_NEW('IDLgrModel')

2、;創(chuàng)建多邊形oPoly=OBJ_NEW('IDLgrPolygon')oView.add,oModel&oModel.add,oPolyverts=0,0,0,1,0,0,1,1,0,0,1,0,0,0,1,1,0,1,1,1,1,0,1,1 ;設置立方體頂點connect=4,0,1,2,3,4,0,1,5,4,4,1,2,6,5,4,2,3,7,6,4,3,0,4,7,4,4,5,6,7;將一個立方體各個定點依次標明，依次連接八個面構成一個閉合的立方體oPoly.setproperty,data=verts,polygons=connect,style=1oMode

3、l.rotate,1,0,1,60;設置立體圖像顯示出來的角度為60度oWindow.draw,oViewvertscolor=0,0,0,1,0,0,1,1,0,0,1,0,0,0,1,1,0,1,1,1,0,0,1,1*255 ;用RGB成像依次構造立方體各個部位的不同顏色oPoly.setproperty,vert_color=vertscolor,style=2,shading=1oWindow.draw,oViewend運行出來的結果如下圖：任務二、讀取head.dat中的數據，進行體數據的顯示，并實現切片提取操作。IDL代碼如下：pro volume device,decompos

4、ed = 1 file = filepath('head.dat',subdirectory = 'examples','data') ;讀取體數據到volume數組 volume = read_binary(file,data_dims = 80,100,57) ;利用XVolume交互顯示數據 xvolume,volume,/interpolate slice = extract_slice(volume,40,40,40,50,28,30,0,0) window,0,xsize = 40,ysize = 40 tv,slice ;顯示垂直切

5、片 window,1,xsize = 100,ysize = 57 tv,volume23,*,* window,2,xsize = 80,ysize = 57 tv,volume*,20,* ;顯示水平切片 window,3,xsize = 80,ysize = 100 tv,volume*,*,20end顯示結果如下 任務三、讀取elev_t.jpg圖像以及elevbin.dat中的DEM數據：利用直接圖形法顯示DEM曲面利用對象圖形法創(chuàng)建曲面對象和紋理對象并進行疊加顯示代碼如下：imageFile=FILEPATH('elev_t.jpg',SUBDIRECTORY=&#

6、39;examples','data') ;讀取圖像文件 READ_JPEG,imageFile,imagedemFile=FILEPATH('elevbin.dat',SUBDIRECTORY='examples','data') ;讀取DEM數據dem=READ_BINARY(demFile,DATA_DIMS=64,64)dem=CONGRID(dem,128,128,/INTERP)DEVICE,DECOMPOSED=0,RETAIN=2;TITLEWINDOW,0,TITLE='Elevation Dat

7、a'&SHADE_SURF,demoModel=OBJ_NEW('IDLgrModel')oView=OBJ_NEW('IDLgrView')oWindow=OBJ_NEW('IDLgrWindow',RETAIN=2,COLOR_MODEL=0) oSurface=OBJ_NEW('IDLgrSurface',dem,STYLE=2)oImage=OBJ_NEW('IDLgrImage',image,INTERLEAVE=0,/INTERPOLATE) oSurface->GetProper

8、ty,XRANGE=xr,YRANGE=yr,ZRANGE=zr;計算歸一化顯示比例，并在各個方向平移-0.5，從而使圖像居中xs=NORM_COORD(xr)&xs0=xs0-0.5ys=NORM_COORD(yr)&ys0=ys0-0.5zs=NORM_COORD(zr)&zs0=zs0-0.5oSurface->SetProperty,XCOORD_CONV=xs,YCOORD_CONV=ys,ZCOORD=zs ;TEXTURE_MAPoSurface->SetProperty,TEXTURE_MAP=oImage,color=255,255,255

9、 oModel->Add,oSurface&oView->Add,oModeloModel->ROTATE,1,0,0,-90&oModel->ROTATE,0,1,0,30oModel->ROTATE,1,0,0,30&oWindow->Draw,oView ;OBJ_DESTROYXOBJVIEW,oModel,/BLOCK,SCALE=1&OBJ_DESTROY,oView,oImage end 作業(yè)二、任務一、繪制函數z=x*exp-x2+ y2，其中-2x,y2的網格曲面圖。IDL代碼如下：pro test_surf

10、ace ;對x,y采樣 x = (findgen(41) - 20)/10 y = (findgen(41) - 20)/10 ;對x,y網格化 temp_x = make_array(n_elements(y),value = 1) temp_y = make_array(n_elements(x),value = 1) xx = x # temp_x yy = temp_y # y ;計算函數的值 z = xx * exp(-xx2 - yy2) ;繪制曲面 surface,zend 得到結果如圖所示：任務二、讀取head.dat中的三維動畫數組，播放該動畫，并存儲該動畫的像素映射圖IDL

11、代碼如下：pro test_animate file = filepath('head.dat',subdirectory = 'examples','data') ;讀取三維動畫數據head = read_binary(file,data_dims = 80,100,57);xinteranimate,set = 80,100,57,/showload ;初始化動畫工具 for i = 0,56 do xinteranimate,frame = i,image = head*,*,i ;將數組加載到動畫工具的緩沖區(qū) xinteranimate,

12、50,/keep_pixmaps ;播放該動畫，并存儲像素映射圖End得到結果如下圖：任務三、載入avhrr.png中的數據，并加以顯示將其轉換到“Interrupted Goode”投影坐標系下，并用iimage命令顯示。將第一問中得到的地圖轉換到“Mollweide”投影坐標系下，并用iimage命令顯示。IDL代碼如下：pro test_projection file = filepath('avhrr.png',subdirectory = 'examples','data') ;讀取數據 data = read_png(file,r,g

13、,b) red0 = rebin(rdata,360,180) green0 = rebin(gdata,360,180) blue0 = rebin(bdata,360,180) ;對原始數據進行重采樣 iimage,red = red0,green = green0,blue = blue0,dimensions = 500,600,view_grid = 1,3 ;數據顯示 smap = map_proj_init('Interrupted Goode') ;創(chuàng)建Interrupted Goode投影 red1 = map_proj_image(red0,map_stru

14、cture = smap,mask = mask,uvrange = uvrange,xindex = xindex,yindex = yindex) green1 = map_proj_image(green0,xindex = xindex,yindex = yindex) blue1 = map_proj_image(blue0,xindex = xindex,yindex = yindex) ;投影轉換 iimage,red = red1,green = green1,blue = blue1,alpha = mask*255b,/view_next ;顯示轉換后的地圖 mapstruct = map_proj_init('Mollweide',/gctp) ;創(chuàng)建Mollweide投影 red2 = map_proj_image(red1,uvrange,image_structure = smap,map_structure = mapstruct,mask = mask,xindex = xindex2,yindex = yindex2)