如何用IDL处理Shapefile数据

转载：http://blog.sina.com.cn/s/blog_6e51df7f010104e9.html

遥感和GIS不分家，IDL擅长处理遥感数据，但偶尔也需要用来处理一些GIS数据，不过还好IDL能处理Shapefile数据。使用IDL处理shapefile格式，需要了解IDLffShape对象，IDL帮助中有一些说明和代码，但过于简单，不熟悉的人很难上手，现对几个关键点进行说明：
一、读取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+GDAL