使用IDL处理shapefile格式，需要了解IDLffShape对象，IDL帮助中有一些说明和代码，但过于简单，不熟悉的人很难上手，现对几个关键点进行说明：感和GIS不分家，IDL擅长处理遥感数据，但偶尔也需要用来处理一些GIS数据，不过还好IDL能处理Shapefile数一、读取shapefile文件1.首先要打开文件
我们用Arcview带的数据做例子吧，就用那个国界数据吧。
创建和销毁idlffshape分别使用的是IDL处理对象的通用命令OBJ_new和Obj_Destroy,每建立一个对象都要记着要销毁，否则会出现内存不足问滑备题。
pro readshapefile
shapefile='C:\ESRI\ESRIDATA\WORLD\country.shp' 定义shape文件位置
oshp=Obj_New('IDLffShape',filename)
print oshp
中间处理代码
Obj_destroy,oshp ；销毁一个shape对象
end
如果读取错误，oshp会返回-1，否则得话返回的就是一个结构体。
2.获取整体描述信息
读完shape对象后，就需要读几何数据和属性表数据了。Shapefile数据由几何体(或郑简实体Entity)和属性表两部分组成，而几何体一般又包括点（point）、线（polyline）和多边形（polygon）（当然也有其它类型，但不常用）。属性表包括属性表结构、字段个数和记录个数，属性表记录信丛毁数与实地必须一一对应，属性表的结构又包含字段名，字段类型，字段长度和精确度。
在IDL读取数据前，需要了解一些全局属性，知道有多少个几何体和记录，属性表中有多少个字段，就需要用GetProperty方法，它查询shape文件的属性，包括实体类型，实体个数，属性表结构，属性表字段个数，记录数等，代码如下:
pro readshapefile
shapefile='C:\ESRI\ESRIDATA\WORLD\country.shp' 定义shape文件位置
oshp=Obj_New('IDLffShape',filename)
oshp->getproperty,n_entities=n_ent,Attribute_info=attr_info,n_attributes=n_attr,Entity_type=ent_type
print,'实体个数:'n_ent
print,'属性表字段数:',n_attr
print,'实体类型代码:',ent_type
Obj_destroy,oshp ；销毁一个shape对象
end
3.再了解几个概念
bouns
存储的是每个实体的范围，是一个有8个元素的数组（[x0,x1,x2,x3,x4,x5,x6,x7]）,其中
x0 —最小x值，x1 —最小y值，x2最小z值（高度），x3 — 最小M值（测量值，一般不用）x4-x7就是最大值了。注意，这里是每个实体的范围，而不是整个地图的范围，所以如果要求整个地图的范围，还需要整个再求一遍最大值和最小值。
bounds还有一个作用，就是存储点的坐标，点类型数据没有安排另外的对象来存储，直接用bounds来管理。
VERTICES
线和面实体的所有坐标都存在这里面，是一个指针型数组，存储的是实体的所有拐点.读的时候比较容易，写入的时候，需要先建一个指针变量将坐标赋值到指针变量,然后将指针变量赋值给vertices.数组结构如下：
[[x,y],[x,y],[x,y],[x,y],[x,y],[x,y]](如果是2维的话)
[[x,y,z],[x,y,z],[x,y,z],[x,y,z],[x,y,z],[x,y,z]](如果是3维的话)
N_VERTICES
拐点的个数，不需要解释了.
N_PARTS和Parts
处理复杂对象的需要注意了，如有内环的多边形。所有的拐点坐标都存在vertices中，Parts也是一个指针数组，存储的是每个弧段的起始索引值。N_PARTS表示有几个弧段.
ISHAPE
表示实体的序号，是一个整形变量，读取的时候一般不需要注意，写的时候需要定义，序号不能重复。
4.开始读坐标了
如果要一次性读取全部实体，可以用ent=oshp->getentity(/all),但大部分时间都需要一个个的处理，就需要用循环
pro readshapefile
shapefile='C:\ESRI\ESRIDATA\WORLD\country.shp' 定义shape文件位置
oshp=Obj_New('IDLffShape',filename)
oshp->getproperty,n_entities=n_ent,Attribute_info=attr_info,n_attributes=n_attr,Entity_type=ent_type
FOR i=0,n_ent-1 do begin 循环
ent=oshp->getentity(i) 读取第i个实体
bounds=ent.bounds 读取实体的边界
n_vert=ent.n_vertices 实体中包括拐点或顶点的个数，只有polyline和polygon具有该属性
vert=*(ent.vertices) 实体的顶点，只有polyline和polygon具有该属性
n_parts=ent.parts 只有polygon具有该属性
part=*(ent.parts) part坐标
输出几何体范围
print,'min x=',bounds[0]
print,'min y=',bounds[1]
print, 'max x=',bound[3]
print, 'max y=',bound[4]
如果是点的话，输出点坐标
print,bounds[0],bounds[1]
如果是线或面的话，输出点坐标
for index in n_vert-1 do begin
print vert[index][0],vert[index][1]
endfor
endfor
Obj_destroy,oshp ；销毁一个shape对象
end
4.读属性属性表的结构
属性表结构存储在Attribute_info中，前面代码已经获得了这个结构体（attr_info）,下面的代码是打印每一个字段的结构
pro readshapefile
shapefile='C:\ESRI\ESRIDATA\WORLD\country.shp' 定义shape文件位置
oshp=Obj_New('IDLffShape',filename)
oshp->getproperty,n_entities=n_ent,Attribute_info=attr_info,n_attributes=n_attr,Entity_type=ent_type
FOR i=0,n_attr-1 do begin 循环
PRINT, '字段序号: ',i
PRINT, '字段名: ', attr_info[i].name
PRINT, '字段类型代码: ', attr_info[i].type
PRINT, '字段宽度: ', attr_info[i].width
PRINT, '精度: ', attr_info[i].precision
endfor
Obj_destroy,oshp ；销毁一个shape对象
end
读属性表中的值
读属性表，跟读取实体有些类似，用GetAttributes方法
pro readshapefile
shapefile='C:\ESRI\ESRIDATA\WORLD\country.shp' 定义shape文件位置
oshp=Obj_New('IDLffShape',filename)
oshp->getproperty,n_entities=n_ent,Attribute_info=attr_info,n_attributes=n_attr,Entity_type=ent_type
FOR i=0,n_ent-1 do begin 循环，n_ent跟记录数是一样的
attr=oshp->GetAttributes(i) 读取第i个记录
for index in n_attr-1 do begin
print attr.(index)
endfor
endfor
Obj_destroy,oshp ；销毁一个shape对象
end
二、写入shapefile数据
写入shapefile的一半过程是，首先初始化idlffshape对象，定义属性表结构，定义实体类型，写入坐标值，写入属性值，最后销毁对象
初始化
写入数据也用Obj_new初始化，不过需要设置输出实体的类型，并设置该数据可写，这里面重要的就是需要知道实体的类型代码，我们常用的就是1，3和5
Point 1
PolyLine 3
Polygon 5
MultiPoint 8
PointZ 11
PolyLineZ 13
PolygonZ 15
MultiPointZ 18
PointM 21
PolyLineM 23
PolygonM 25
MultiPointM 28
MultiPatch 31
对象初始化的代码如下:
pro writeshapefile
shapefile='d:\data\citys.shp'
oshp=obj_new('IDLffshape',new_shapefile,Entity_type=3,/update)
其他代码
obj_destroy,oshp
创建属性表结构和实体类型
对所有类型的实体，创建属性表的方法都已一样的。用AddAttribute方法，一般用法为：
oshp->AddAttribute, 字段名称,字段类型,字段宽度[, PRECISION=integer]
精度只有浮点和双精度等情况下采用，字符和整形可以缺省，也可以设置为0
关键的还是需要知道常用的几种字段类型
3 Longword integer
5 Double-precision floating-point
7 String
没错只有三种，这不是idl的错，shapefile只定义了这三种
定义实体类型的方法比较简单：
entNew = {IDL_SHAPE_ENTITY}
entNew.SHAPE_TYPE = 1 1为实体类型，表示点
写入实体和属性
这两个过程一般同时进行的,用代码表示吧：
Pro writepoint
shapefile='d:\test\citys.shp'
oshp=OBJ_NEW('IDLffshape',shapefile,Entity_type=1,/update)
定义实体类型
entNew = {IDL_SHAPE_ENTITY}
entNew.SHAPE_TYPE = 1 1为实体类型，表示点
添加坐标,加那个地方呢，我爱北京天安门吧
entNew.ISHAPE=0
entNew.BOUNDS[0] = 116.391188
entNew.BOUNDS[1] = 39.904546
entNew.BOUNDS[2] = 0.00000000
entNew.BOUNDS[3] = 0.00000000
entNew.BOUNDS[4] = 116.391188
entNew.BOUNDS[5] = 39.904546
entNew.BOUNDS[6] = 0.00000000
entNew.BOUNDS[7] = 0.00000000
entNew.N_VERTICES = 1
加属性了
先定义属性表结构
oshp->AddAttribute,'id',3,8,PRECISION=0
oshp->AddAttribute,'name',7,20,PRECISION=0
oshp->AddAttribute,'longitude',5,8,PRECISION=4
oshp->AddAttribute,'latitude',5,8,PRECISION=4
还要把实体写入到shp对象中
oshp ->PutEntity, entNew
获得属性表结构对象
new_attr = oshp ->GetAttributes(/ATTRIBUTE_STRUCTURE)
new_attr.ATTRIBUTE_0 = 1
new_attr.ATTRIBUTE_1 = '北京天安门'
new_attr.ATTRIBUTE_2 = 116.3911
new_attr.ATTRIBUTE_3 = 39.904546
把属性写入到shp对象中
oshp ->SetAttributes,0,new_attr这里面的0是指实体的索引值，等于entNew.ISHAPE
再加一个吧，就兰州了
entNew.BOUNDS = [103.867694,36.048088,0,0,103.867694,36.048088,0,0]
new_attr.(0)=2
new_attr.(1)='兰州'
new_attr.(2)=103.8676
new_attr.(3)=36.0480
oshp ->PutEntity, entNew
oshp ->SetAttributes,1,new_attr
OBJ_DESTROY,oshp
print,'end'
End
以上是加入点类型的数据,比较简单，来个复杂点的，加两个多边形吧
pro writepolygon
shapefile='d:\test\Forbidden_City.shp'
oshp=obj_new('IDLffshape',shapefile,Entity_type=5,/update)
定义实体类型
entNew = {IDL_SHAPE_ENTITY}
entNew.SHAPE_TYPE = 5
添加坐标
coor=[[116.3852041484393,39.9214192520002],$
[116.3856922399481,39.91151453640624],$
[116.3960721525212,39.9118040463524],$
[116.3955102491546,39.92183809311693],$
[116.3852041484393,39.9214192520002]]
entNew.ISHAPE=0
entNew.BOUNDS[0] = min(coor[0,*])
entNew.BOUNDS[1] = min(coor[1,*])
entNew.BOUNDS[2] = 0.00000000
entNew.BOUNDS[3] = 0.00000000
entNew.BOUNDS[4] = max(coor[0,*])
entNew.BOUNDS[5] = max(coor[1,*])
entNew.BOUNDS[6] = 0.00000000
entNew.BOUNDS[7] = 0.00000000
pvertice=coor
entNew.VERTICES=PTR_NEW(pvertice,/no_copy)
entNew.N_VERTICES = 5
还要把实体写入到shp对象中
oshp ->PutEntity, entNew
加属性
先定义属性表结构
oshp->AddAttribute,'id',3,8,PRECISION=0
oshp->AddAttribute,'name',7,20,PRECISION=0
获得属性表结构对象
new_attr = oshp ->GetAttributes(/ATTRIBUTE_STRUCTURE)
new_attr.ATTRIBUTE_0 = 1
new_attr.ATTRIBUTE_1 = 'river'
把属性写入到shp对象中
oshp ->SetAttributes,0,new_attr这里面的0是指实体的索引值，等于entNew.ISHAPE
coor=[[116.3858622895445,39.92099455865304],$
[116.3863498312803,39.91211319734286],$
[116.3952884054441,39.91246510632352],$
[116.3948307781919,39.92118603918453],$
[116.3858622895445,39.92099455865304]]
entNew.ISHAPE=1
entNew.BOUNDS = [min(coor[0,*]),min(coor[1,*]),0,0,max(coor[0,*]),max(coor[1,*]),0,0]
pvertice=coor
entNew.VERTICES=PTR_NEW(pvertice,/no_copy)
entNew.N_VERTICES = (size(coor))[2]
entNew.N_Parts=2
P_parts=[0,5,9]
entNew.Parts=Ptr_new(P_parts,/no_copy)
还要把实体写入到shp对象中
oshp ->PutEntity, entNew
加属性
new_attr.ATTRIBUTE_0 = 1
new_attr.ATTRIBUTE_1 = 'Forbidden_City'
把属性写入到shp对象中
oshp ->SetAttributes,1,new_attr这里面的0是指实体的索引值，等于entNew.ISHAPE
obj_destroy,oshp
print,'end'
end
三、获得完整示例代码
行文仓促，文中的代码可能有误，撰写了3个较为完整的示例代码,放在了我的代码库中，感兴趣的话可以到googleCode上下载。下载地址为:
http://code.google.com/p/datatools/source/browse/#svn/trunk/IDL/Shapefile
3个示例代码的名称分别为:
readshapefile.pro
writepoint.pro
writepolygon.pro
另外,大家还从该站点获得利用IDL创建aster图像索引图的程序.
四、后记
用IDL处理矢量数据，始终比较复杂，能够处理的格式也有限。如果跟矢量数据大交道比较多的话，建议尝试Python
能拆此兄_对氐牡谝徊?,是和安装流程一样的旅袭界面。 然后,会有对话框询问你是否真的“want to”卸载。 在回答之后,会进入一个红色的界扒唯面,提示的是“removing” 这是进度
从开码笑始做这个课题到现在就没少用IDL读FITS文件。这个方面用mrdfits比较容易，基本就是一行搞定数据，几行搞定文件头，用了不知多少次。其实在读FITS的时候就在想，把写FITS也搞明白吧，不过惰性太大，一直都回避这个问题。 今天合作者建议我罩慎把数据平滑一下重新计算。我用的那个程序的输入就是一个FITS文件，这就意味着我需要重新写一个平滑后的FITS文件，于是今天不得不去看看怎么写FITS文件了。原来知道和mrdfits对应的有mwrfits，专门写FITS文件的。我有一个数组迟闷含a和文件头head，于是按照说明里写 IDL>mwrfits,a,'out.fits',head 这样倒是可以生成一个FITS文件，查看了也正常，可是我用来处理的那个程序就是不认。不得已，参考了一下别人的程序，用writefits IDL>writefits, 'out.fits', a, head 这样生成的FITS文件就能被识别了。原因为何，有待研究。 在文件头某些信息改变的情况下还需要改一下文件头里的参数，可以用fxaddpar，例如改变参数'NAXIS'的值 IDL>fxaddpar,'NAXIS',2欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）