如何用IDL处理Shapefile数据使用IDL处理shapefile格式，需要了解IDLffShape对象，IDL帮助中有一些说明和代码，但过于简单，不熟悉的人很难上手，现对几个关键点进行说明：感和GIS不分家，IDL擅长处理遥感数据，但偶尔也需要用来处理一些GIS数据，不过还好IDL能处理Shapefile数一、读取shapefile文件1.首先要打开文件我们用Arcview带的数据做例子吧，就用那个国界数据吧。创建和销毁idlffshape分别使用的是IDL处理对象的通用命令OBJ_new和Obj_Destroy,每建立一个对象都要记着要销毁，否则会出现内存不足问题。pro readshapefileshapefile="C:ESRIESRIDATAWORLDcountry.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 readshapefileshapefile="C:ESRIESRIDATAWORLDcountry.shp" ;定义shape文件位置oshp=Obj_New("IDLffShape",filename)oshp->getproperty,n_entities=n_ent,Attribute_info=attr_info,n_attributes=n_attr,Entity_type=ent_typeprint,"实体个数:"n_entprint,"属性表字段数:",n_attrprint,"实体类型代码:",ent_typeObj_destroy,oshp ；销毁一个shape对象end3.再了解几个概念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 readshapefileshapefile="C:ESRIESRIDATAWORLDcountry.shp" ;定义shape文件位置oshp=Obj_New("IDLffShape",filename)oshp->getproperty,n_entities=n_ent,Attribute_info=attr_info,n_attributes=n_attr,Entity_type=ent_typeFOR 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 beginprint vert[index][0],vert[index][1]endforendforObj_destroy,oshp ；销毁一个shape对象end4.读属性属性表的结构属性表结构存储在Attribute_info中，前面代码已经获得了这个结构体（attr_info）,下面的代码是打印每一个字段的结构pro readshapefileshapefile="C:ESRIESRIDATAWORLDcountry.shp" ;定义shape文件位置oshp=Obj_New("IDLffShape",filename)oshp->getproperty,n_entities=n_ent,Attribute_info=attr_info,n_attributes=n_attr,Entity_type=ent_typeFOR i=0,n_attr-1 do begin ;循环PRINT, "字段序号: ",iPRINT, "字段名: ", attr_info[i].namePRINT, "字段类型代码: ", attr_info[i].typePRINT, "字段宽度: ", attr_info[i].widthPRINT, "精度: ", attr_info[i].precisionendforObj_destroy,oshp ；销毁一个shape对象end读属性表中的值读属性表，跟读取实体有些类似，用GetAttributes方法pro readshapefileshapefile="C:ESRIESRIDATAWORLDcountry.shp" ;定义shape文件位置oshp=Obj_New("IDLffShape",filename)oshp->getproperty,n_entities=n_ent,Attribute_info=attr_info,n_attributes=n_attr,Entity_type=ent_typeFOR i=0,n_ent-1 do begin ;循环，n_ent跟记录数是一样的attr=oshp->GetAttributes(i) ;读取第i个记录for index in n_attr-1 do beginprint attr.(index)endforendforObj_destroy,oshp ；销毁一个shape对象end二、写入shapefile数据写入shapefile的一半过程是，首先初始化idlffshape对象，定义属性表结构，定义实体类型，写入坐标值，写入属性值，最后销毁对象初始化写入数据也用Obj_new初始化，不过需要设置输出实体的类型，并设置该数据可写，这里面重要的就是需要知道实体的类型代码，我们常用的就是1，3和5Point 1PolyLine 3Polygon 5MultiPoint 8PointZ 11PolyLineZ 13PolygonZ 15MultiPointZ 18PointM 21PolyLineM 23PolygonM 25MultiPointM 28MultiPatch 31对象初始化的代码如下:pro writeshapefileshapefile="d:datacitys.shp"oshp=obj_new("IDLffshape",new_shapefile,Entity_type=3,/update);;其他代码obj_destroy,oshp创建属性表结构和实体类型对所有类型的实体，创建属性表的方法都已一样的。用AddAttribute方法，一般用法为：oshp->AddAttribute, 字段名称,字段类型,字段宽度[, PRECISION=integer] 精度只有浮点和双精度等情况下采用，字符和整形可以缺省，也可以设置为0关键的还是需要知道常用的几种字段类型3 Longword integer5 Double-precision floating-point7 String没错只有三种，这不是idl的错，shapefile只定义了这三种定义实体类型的方法比较简单：entNew = {IDL_SHAPE_ENTITY}entNew.SHAPE_TYPE = 1 ;1为实体类型，表示点写入实体和属性这两个过程一般同时进行的,用代码表示吧：Pro writepointshapefile="d: estcitys.shp"oshp=OBJ_NEW("IDLffshape",shapefile,Entity_type=1,/update);定义实体类型entNew = {IDL_SHAPE_ENTITY}entNew.SHAPE_TYPE = 1 ;1为实体类型，表示点;添加坐标,加那个地方呢，我爱北京天安门吧entNew.ISHAPE=0entNew.BOUNDS[0] = 116.391188entNew.BOUNDS[1] = 39.904546entNew.BOUNDS[2] = 0.00000000entNew.BOUNDS[3] = 0.00000000entNew.BOUNDS[4] = 116.391188entNew.BOUNDS[5] = 39.904546entNew.BOUNDS[6] = 0.00000000entNew.BOUNDS[7] = 0.00000000entNew.N_VERTICES = 1;加属性了;先定义属性表结构oshp->AddAttribute,"id",3,8,PRECISION=0oshp->AddAttribute,"name",7,20,PRECISION=0oshp->AddAttribute,"longitude",5,8,PRECISION=4oshp->AddAttribute,"latitude",5,8,PRECISION=4;还要把实体写入到shp对象中oshp -> PutEntity, entNew;获得属性表结构对象new_attr = oshp ->GetAttributes(/ATTRIBUTE_STRUCTURE)new_attr.ATTRIBUTE_0 = 1new_attr.ATTRIBUTE_1 = "北京天安门"new_attr.ATTRIBUTE_2 = 116.3911new_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)=2new_attr.(1)="兰州"new_attr.(2)=103.8676new_attr.(3)=36.0480oshp -> PutEntity, entNewoshp -> SetAttributes,1,new_attr;OBJ_DESTROY,oshpprint,"end"End以上是加入点类型的数据,比较简单，来个复杂点的，加两个多边形吧pro writepolygonshapefile="d: estForbidden_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=0entNew.BOUNDS[0] = min(coor[0,*])entNew.BOUNDS[1] = min(coor[1,*])entNew.BOUNDS[2] = 0.00000000entNew.BOUNDS[3] = 0.00000000entNew.BOUNDS[4] = max(coor[0,*])entNew.BOUNDS[5] = max(coor[1,*])entNew.BOUNDS[6] = 0.00000000entNew.BOUNDS[7] = 0.00000000pvertice=coorentNew.VERTICES=PTR_NEW(pvertice,/no_copy)entNew.N_VERTICES = 5;还要把实体写入到shp对象中oshp -> PutEntity, entNew;加属性;先定义属性表结构oshp->AddAttribute,"id",3,8,PRECISION=0oshp->AddAttribute,"name",7,20,PRECISION=0;获得属性表结构对象new_attr = oshp ->GetAttributes(/ATTRIBUTE_STRUCTURE)new_attr.ATTRIBUTE_0 = 1new_attr.ATTRIBUTE_1 = "river";把属性写入到shp对象中oshp -> SetAttributes,0,new_attr;这里面的0是指实体的索引值，等于entNew.ISHAPEcoor=[[116.3858622895445,39.92099455865304],$[116.3863498312803,39.91211319734286],$[116.3952884054441,39.91246510632352],$[116.3948307781919,39.92118603918453],$[116.3858622895445,39.92099455865304]]entNew.ISHAPE=1entNew.BOUNDS = [min(coor[0,*]),min(coor[1,*]),0,0,max(coor[0,*]),max(coor[1,*]),0,0]pvertice=coorentNew.VERTICES=PTR_NEW(pvertice,/no_copy)entNew.N_VERTICES = (size(coor))[2]entNew.N_Parts=2P_parts=[0,5,9]entNew.Parts=Ptr_new(P_parts,/no_copy);还要把实体写入到shp对象中oshp -> PutEntity, entNew;加属性new_attr.ATTRIBUTE_0 = 1new_attr.ATTRIBUTE_1 = "Forbidden_City";把属性写入到shp对象中oshp -> SetAttributes,1,new_attr;这里面的0是指实体的索引值，等于entNew.ISHAPEobj_destroy,oshpprint,"end"end三、获得完整示例代码行文仓促，文中的代码可能有误，撰写了3个较为完整的示例代码,放在了我的代码库中，感兴趣的话可以到googleCode上下载。下载地址为:http://code.google.com/p/datatools/source/browse/#svn/trunk/IDL/Shapefile3个示例代码的名称分别为:readshapefile.prowritepoint.prowritepolygon.pro另外,大家还从该站点获得利用IDL创建aster图像索引图的程序.四、后记用IDL处理矢量数据，始终比较复杂，能够处理的格式也有限。如果跟矢量数据大交道比较多的话，建议尝试Python
