如何用IDL处理Shapefile数据

如何用IDL处理Shapefile数据

使用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

  本文介绍基于 ENVI 软件，对 不含有任何地理参考信息 的栅格遥感影像添加 地理坐标系 或 投影坐标系 等 地理参考信息 的方法。

  我们先来看一下本文需要实现的需求。现有以下两景遥感影像，其位于不同的空间位置；但由于二者均不含任何地理参考信息，导致其在 ENVI 软件中打开后会自动重叠在一起；如下图所示。

  那么我们就以其中一景遥感影像为例，对其添加地理参考信息。

  明确了具体需求，接下来就可以开始 *** 作。首先，我们在 ENVI 软件中打开对应的两景遥感影像；其次，在需要添加地理参考信息的图像名称处右键，选择“ View Metadata ”。

  d出如下所示的元数据浏览窗口。

  这里我们需要注意：如果大家打开的元数据浏览窗口的左侧列表中含有“ Map Info ”这个选项，那么我们直接单击，将其打开，并选择“ Edit Metadata ”进行本文后续的 *** 作即可；而如果是像本文中一样，遥感影像元数据窗口没有“ Map Info ”这个选项，那么我们就需要点击上图中“ Edit Metadata ”，随后在d出的“ Set Raster Metadata ”窗口中点击左上角的“ Add... ”选项，将d出另一个“ Add Metadata Items ”窗口。

  随后，在“ Add Metadata Items ”窗口中选择“ Spatial Reference ”选项，并点击“ OK ”。

  稍等片刻（这段时间中， ENVI 软件可能会出现如同卡死一般的闪烁，大家不用管，多等待一会即可），可以看到在“ Set Raster Metadata ”窗口中，已经出现如下所示的“ Spatial Reference ”选项。

  我们对“ Spatial Reference ”选项进行编辑即可。其中，首先需要选择地理坐标系或投影坐标系的种类；其次配置遥感图像的空间分辨率，也就是每一个像元的 X 大小和 Y 大小；再次，“ Tie Point ”中，前两个选项（“ Pixel X ”与“ Pixel Y ”）为我们 参考点 （这个参考点具体是什么，我们稍后会介绍）在图像中的位置，后两个选项（“ Map X ”与“ Map Y ”）则是该 参考点 实际的空间位置——如果我们选择的是 地理坐标系 ，那么这里就是实际的 经纬度 ；如果我们选择的是 投影坐标系 ，那么这里就是实际的 投影数值 。最后，配置坐标系的旋转角度，一般填 0 就可以。我在这里只是做一个示范，因此下图中的各参数也都是乱填的，大家依据实际情况来配置各参数即可。

  关于这个“ 参考点 ”，这里有必要再多提几句。 参考点 其实就是该图像中，某一个已知 实际空间坐标信息 、已知 其在图像中位置 的点；我们需要将这个点在 图像中的位置 （以行列号的形式表示，行数与列数均从 0 开始算起，遥感影像左上角的 像元 的左上角 点 为第 0 行第 0 列）与该点在 实际中的位置 输入进去，然后软件再依据我们所选择的坐标系与图像空间分辨率，对图像中每一个像元的空间位置进行计算，从而最终生成一个带有地理参考信息的栅格图像。

  随后点击“ OK ”，即可完成对该图像的地理参考信息的配置。我们再一次查看该图像的元数据，可以发现此时其已经含有“ Map Info ”这个选项，且其中的参数都已经是刚刚我们设定的参数了。

  这里我们再依据结果图像，来再解释一下参考点的意义。通过上图我们可以知道，我在本文中是将“ Pixel X ”与“ Pixel Y ”均为 0 的这个点作为参考点，并将其空间位置（“ Map X ”与“ Map Y ”）均设置为 1 ；那么在结果图中，我们通过 Crosshairs 功能、 Cursor Value 功能确定一下该点的位置，如下图所示；可以看到“ Pixel X ”与“ Pixel Y ”均为 0 的这个点（图中黄色圈内），其经、纬度就近似为1°与1°（之所以是近似，是因为我也不是完全选中了这个参考点，而是近似选中）。

---