地学分析与三维可视化实习一.docx
《地学分析与三维可视化实习一.docx》由会员分享，可在线阅读，更多相关《地学分析与三维可视化实习一.docx（12页珍藏版）》请在冰豆网上搜索。
地学分析与三维可视化实习一
地学分析与三维可视化实习报告一
XX：
方志策
班级：
064131
学号：
20211003574
作业一、
任务一、利用对象图形法创立一个三维立方体，并将各顶点设置为不同的颜色。
IDL代码如下：
oWindow=OBJ_NEW（'IDLgrWindow',dimension=[400,400]）
;建立一个图像窗口，初始设置图像大小为400*400
oView=OBJ_NEW（'IDLgrView',viewPlane_Rect=[-1,-1,3,3],zClip=[3,-3],eye=2）
oModel=OBJ_NEW（'IDLgrModel'）
;创立多边形
oPoly=OBJ_NEW（'IDLgrPolygon'）
oView.add,oModel&oModel.add,oPoly
verts=[[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=1
oModel.rotate,[1,0,1],60
;设置立体图像显示出来的角度为60度
oWindow.draw,oView
vertscolor=[[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=1
oWindow.draw,oView
end
运行出来的结果如下列图：
任务二、读取head.dat中的数据，进展体数据的显示，并实现切片提取操作。
IDL代码如下：
provolume
device,decomposed=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
;显示垂直切片
window,1,xsize=100,ysize=57
tv,volume[23,*,*]
window,2,xsize=80,ysize=57
tv,volume[*,20,*]
;显示水平切片
window,3,xsize=80,ysize=100
tv,volume[*,*,20]
end
显示结果如下
任务三、法显示DEM曲面
利用对象图形法创立曲面对象和纹理对象并进展叠加显示
代码如下：
imageFile=FILEPATH（'elev_t.jpg',SUBDIRECTORY=['examples','data']）
;读取图像文件
READ_JPEG,imageFile,image
demFile=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
;TITLE
WINDOW,0,TITLE='ElevationData'&SHADE_SURF,dem
oModel=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->GetProperty,XRANGE=xr,YRANGE=yr,ZRANGE=zr
;计算归一化显示比例，并在各个方向平移-0.5，从而使图像居中
xs=NORM_COORD（xr）&xs[0]=xs[0]-0.5
ys=NORM_COORD（yr）&ys[0]=ys[0]-0.5
zs=NORM_COORD（zr）&zs[0]=zs[0]-0.5
oSurface->SetProperty,XCOORD_CONV=xs,YCOORD_CONV=ys,ZCOORD=zs;TEXTURE_MAP
oSurface->SetProperty,TEXTURE_MAP=oImage,color=[255,255,255]
oModel->Add,oSurface&oView->Add,oModel
oModel->ROTATE,[1,0,0],-90&oModel->ROTATE,[0,1,0],30
oModel->ROTATE,[1,0,0],30&oWindow->Draw,oView;OBJ_DESTROY
XOBJVIEW,oModel,/BLOCK,SCALE=1&OBJ_DESTROY,[oView,oImage]
end
作业二、
任务一、绘制函数
，其中
的网格曲面图。
IDL代码如下：
protest_surface
;对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（-xx^2-yy^2）
;绘制曲面
surface,z
end
得到结果如下图：
任务二、读取head.dat中的三维动画数组，播放该动画，并存储该动画的像素映射图
IDL代码如下：
protest_animate
file=filepath（'head.dat',subdirectory=['examples','data']）
;读取三维动画数据
head=read_binary（file,data_dims=[80,100,57]）;
xinteranimate,set=[80,100,57],/showload
;初始化动画工具
fori=0,56doxinteranimate,frame=i,image=head[*,*,i]
;将数组加载到动画工具的缓冲区
xinteranimate,50,/keep_pixmaps
;播放该动画，并存储像素映射图
End
得到结果如下列图：
任务三、载入avhrr.png中的数据，并加以显示
将其转换到“InterruptedGoode〞投影坐标系下，并用iimage命令显示。
将第一问中得到的地图转换到“Mollweide〞投影坐标系下，并用iimage命令显示。
IDL代码如下：
protest_projection
file=filepath（'avhrr.png',subdirectory=['examples','data']）
;读取数据
data=read_png（file,r,g,b）
red0=rebin（r[data],360,180）
green0=rebin（g[data],360,180）
blue0=rebin（b[data],360,180）
;对原始数据进展重采样
iimage,red=red0,green=green0,blue=blue0,dimensions=[500,600],view_grid=[1,3]
;数据显示
smap=map_proj_init（'InterruptedGoode'）
;创立InterruptedGoode投影
red1=map_proj_image（red0,map_structure=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）
;创立Mollweide投影
red2=map_proj_image（red1,uvrange,image_structure=smap,map_structure=mapstruct,mask=mask,xindex=xindex2,yindex=yindex2）
green2=map_proj_image（green1,xindex=xindex2,yindex=yindex2）
blue2=map_proj_image（blue1,xindex=xindex2,yindex=yindex2）
;投影转换
iimage,red=red2,green=green2,blue=blue2,alpha=mask*255b,/view_next
;显示投影后的地图
End
得到结果如下：
升级会员 