shp文件的读取

来自：http://blog.sina.com.cn/s/blog_6e51df7f0100ya11.html

首先了解一下shp文件的一些简单知识

Shapefile文件是美国环境系统研究所（ESRI）所研制的GIS文件系统格式文件，是工业标准的矢量数据文件。Shapefile将空间特征表中的非拓扑几何对象和属性信息存储在数据集中，特征表中的几何对象存为以坐标点集表示的图形文件—SHP文件，Shapefile文件并不含拓扑（Topological）数据结构。一个Shape文件包括三个文件：一个主文件(*.shp)，一个索引文件(*.shx)，和一个dBASE(*.dbf)表。主文件是一个直接存取，变长度记录的文件，其中每个记录描述构成一个地理特征（Feature）的所有vertices坐标值。在索引文件中，每条记录包含对应主文件记录距离主文件头开始的偏移量，dBASE表包含SHP文件中每一个Feature的特征属性，表中几何记录和属性数据之间的一一对应关系是基于记录数目的ID。在dBASE文件中的属性记录必须和主文件中的记录顺序是相同的。图形数据和属性数据通过索引号建立一一对应的关系。

Shapefile中坐标文件（.shp）由固定长度的文件头和接着的变长度空间数据记录组成。文件头由100字节的说明信息组成的，主要说明文件的长度、Shape类型、整个Shape图层的范围等等，这些信息构成了空间数据的元数据。在导入空间数据时首先要读入文件头获取Shape文件的基本信息，并以此信息为基础建立相应的元数据表。而变长度空间数据记录是由固定长度的记录头和变长度记录内容组成，其记录结构基本类似，每条记录都有记录头和记录内容组成（空间坐标对）。记录头的内容包括记录号（Record Number）和坐标记录长度（ContentLength）两个记录项，Shapefile文件中的记录号都是从1开始的，坐标记录长度是按16位字来衡量的。记录内容包括目标的几何类型（ShapeType）和具体的坐标记录（X，Y），记录内容因要素几何类型的不同，其具体的内容和格式都有所不同。对于具体的记录主要包括空Shape记录，点记录，线记录和多边形记录。

属性文件（.dbf）用于记录属性信息。它是一个标准的DBF文件，也是由头文件和实体信息两部分构成。其中文件头部分的长度是不定长的，它主要对DBF文件作了一些总体说明，其中最主要的是对这个DBF文件的记录项的信息进行了详细的描述，比如对每个记录项的名称，数据类型，长度等信息都有具体的说明。属性文件的实体信息部分就是一条条属性记录，每条记录都是由若干个记录项构成，因此只要依次循环读取每条记录就可以了。

索引文件（.shx）主要包含坐标文件的索引信息，文件中每个记录包含对应的坐标文件记录距离坐标文件的文件头的偏移量。通过索引文件可以很方便地在坐标文件中定位到指定目标地坐标信息。索引文件也是由文件头和实体信息两部分构成的，其中文件头部分是一个长度固定（100bytes）的记录段，其内容与坐标文件的文件头基本一致。它的实体信息以记录为基本单位，每一条记录包括偏移量（Offset）和记录段长度（Content Length）两个记录项。

接下来我们再看一下代码（c++）：

SHP文件的读取

FeatureClass* CGISMapDoc::ImportShapeFileData( FILE* fpShp, FILE* fpDbf )
{
//读Shp文件头开始
int fileCode = -1;
int fileLength = -1;
int version = -1;
int shapeType = -1;
fread(&fileCode , sizeof(int) , 1 , fpShp) ;
fileCode = ReverseBytes(fileCode) ;

if (fileCode != 9994)
{
CString strTemp ;
strTemp.Format(" WARNING filecode %d ", fileCode );
AfxMessageBox(strTemp);
}

for( int i = 0 ; i < 5 ; i ++ )
fread(&fileCode , sizeof(int) , 1 , fpShp) ;

fread(&fileLength , sizeof(int) , 1 , fpShp) ;
fileLength = ReverseBytes(fileLength) ;

fread(&version , sizeof(int) , 1 , fpShp) ;
fread(&shapeType , sizeof(int) , 1 , fpShp) ;

double tempOriginX , tempOriginY ;
fread( &tempOriginX , sizeof(double) , 1 , fpShp ) ;
fread( &tempOriginY , sizeof(double) , 1 , fpShp ) ;

double xMaxLayer , yMaxLayer ;
fread( &xMaxLayer , sizeof(double) , 1 , fpShp ) ;
fread( &yMaxLayer , sizeof(double) , 1 , fpShp ) ;

double* skip = new double[4] ;
fread( skip , sizeof(double) , 4 , fpShp ) ;
delete []skip ;
skip = 0 ;
//读Shp文件头结束

int uniqueID = this->m_pDataSource->GetUniqueID() ;
FeatureClass* pShpDataSet = 0 ;
//根据目标类型创建相应的图层DataSet。
switch( shapeType )
{
case 1 :
pShpDataSet = (FeatureClass*)&(m_pDataSource->CreateDataSet(uniqueID , POINTDATASET , layerName)) ;
break ;
case 3 :
case 23 :
pShpDataSet = (FeatureClass*)&(m_pDataSource->CreateDataSet(uniqueID , LINEDATASET , layerName)) ;
break ;
case 5 :
pShpDataSet = (FeatureClass*)&(m_pDataSource->CreateDataSet(uniqueID , POLYGONDATASET , layerName)) ;
break ;
}

if ( pShpDataSet == 0 ) return 0;

// 读DBF文件头---------begin------------
struct DBFHeader
{
char m_nValid;
char m_aDate[3];
char m_nNumRecords[4];
char m_nHeaderBytes[2];
char m_nRecordBytes[2];
char m_nReserved1[3];
char m_nReserved2[13];
char m_nReserved3[4];
}dbfheader;

struct DBFFIELDDescriptor
{
char m_sName[10];//应该为char m_sName[11]
char m_nType;
char m_nAddress[4];
char m_nFieldLength;
char m_nFieldDecimal;
char m_nReserved1[2];
char m_nWorkArea;
char m_nReserved2[2];
char m_nSetFieldsFlag;
char m_nReserved3[8];
};

fread(&dbfheader,sizeof(DBFHeader),1,fpDbf);
/*int recordsNum = *((int*)dbfheader.m_nNumRecords);
int headLen = *((short*)dbfheader.m_nHeaderBytes);
int everyRecordLen = *((short*)dbfheader.m_nRecordBytes);

if ( recordsNum == 0 ||  headLen == 0 || everyRecordLen == 0 )
return 0 ;

int fieldCount = (headLen - 1 - sizeof(DBFHeader))/sizeof(DBFFIELDDescriptor);

DBFFIELDDescriptor *pFields = new DBFFIELDDescriptor[fieldCount];
for ( i = 0; i < fieldCount; i ++ )
fread(&pFields[i],sizeof(DBFFIELDDescriptor),1,fpDbf);

char endByte;
fread(&endByte,sizeof(char),1,fpDbf);

if ( endByte != 0x0D)
{
delete []pFields;
pFields = 0;
return 0;
}*/

Fields& fields = pShpDataSet->GetFields();
DBFFIELDDescriptor field ;
BYTE endByte = ' ';
char fieldName[12];
int fieldDecimal, fieldLen, everyRecordLen = 0 ;
while ( !feof(fpDbf) )
{
fread(&endByte,sizeof(BYTE),1,fpDbf);
if ( endByte == 0x0D)   break ;
fread(&field,sizeof(DBFFIELDDescriptor),1,fpDbf);

fieldName[0] = endByte;
for ( i = 0; i < 10; i ++ )
fieldName[i+1] = field.m_sName[i];
fieldName[11] = '/0';

fieldDecimal = field.m_nFieldDecimal;
fieldLen = field.m_nFieldLength;
switch( field.m_nType )
{
case 'C':
fields.AddField(fieldName,fieldName,FIELD_STRING,fieldLen);
break;
case 'F':
fields.AddField(fieldName,fieldName,FIELD_DOUBLE,fieldLen);
break;
case 'N':
{
if ( fieldDecimal == 0 )
fields.AddField(fieldName,fieldName,FIELD_INT,fieldLen);
else fields.AddField(fieldName,fieldName,FIELD_DOUBLE,fieldLen);
}
break;
}
everyRecordLen += fieldLen ;
}
// 读DBF文件头---------end------------

while( !feof(fpShp) )
{
//读记录头开始
int recordNumber = -1 ;
int contentLength = -1 ;
fread(&recordNumber , sizeof(int) , 1 , fpShp) ;
fread(&contentLength , sizeof(int) , 1 , fpShp) ;
recordNumber = ReverseBytes(recordNumber) ;
contentLength = ReverseBytes(contentLength) ;
//读记录头结束

switch( shapeType )
{
case 1: // '/001'
//读取点目标开始
{
Fields &featureFields = pShpDataSet->GetFields();
Feature *pFeature = new Feature(recordNumber , 1 , &featureFields) ;

int pointShapeType ;
fread(&pointShapeType , sizeof(int) , 1 , fpShp) ;
double xValue , yValue ;
fread(&xValue , sizeof(double) , 1 , fpShp) ;
fread(&yValue , sizeof(double) , 1 , fpShp) ;

GeoPoint *pNewGeoPoint = new GeoPoint( xValue , yValue ) ;
pFeature->SetBound(xValue , yValue , 0 , 0 ) ;
pFeature->SetGeometry(pNewGeoPoint) ;
this->LoadAttributeData(pFeature,fpDbf,everyRecordLen);
pShpDataSet->AddRow(pFeature) ;
}
break ;
//读取点目标结束

case 3: // '/003'
//读取线目标开始
{
Fields &featureFields = pShpDataSet->GetFields();
Feature *pFeature = new Feature(recordNumber , 1 , &featureFields) ;

int arcShapeType ;
fread(&arcShapeType , sizeof(int) , 1 , fpShp) ;

double objMinX , objMinY , objMaxX , objMaxY ;
fread(&objMinX , sizeof(double) , 1 , fpShp) ;
fread(&objMinY , sizeof(double) , 1 , fpShp) ;
fread(&objMaxX , sizeof(double) , 1 , fpShp) ;
fread(&objMaxY , sizeof(double) , 1 , fpShp) ;

GeoPolyline *pNewGeoLine = new GeoPolyline();
double width = objMaxX - objMinX ;
double height = objMaxY - objMinY ;
pFeature->SetBound(objMinX , objMinY , width , height) ;

int numParts , numPoints ;
fread(&numParts , sizeof(int) , 1 , fpShp) ;
fread(&numPoints , sizeof(int) , 1 , fpShp) ;
//存储各段线的起点索引
int* startOfPart = new int[numParts] ;
for( int i = 0 ; i < numParts ; i++ )
{
int indexFirstPoint ;
fread(&indexFirstPoint , sizeof(int) , 1 , fpShp) ;
startOfPart[i] = indexFirstPoint ;
}

//处理单个目标有多条线的问题
pNewGeoLine->SetPointsCount( numParts ) ;

for( i = 0 ; i < numParts ; i++ )
{
GeoPoints& points = pNewGeoLine->GetPoints(i) ;
int curPosIndex = startOfPart[i] ;
int nextPosIndex = 0 ;
int curPointCount = 0 ;
if( i == numParts - 1 )
curPointCount = numPoints - curPosIndex ;
else
{
nextPosIndex = startOfPart[i + 1] ;
curPointCount = nextPosIndex - curPosIndex ;
}
points.SetPointCount( curPointCount ) ;
//加载一条线段的坐标
for( int iteratorPoint = 0 ; iteratorPoint < curPointCount ; iteratorPoint ++ )
{
double x , y ;
fread(&x , sizeof(double) , 1 , fpShp) ;
fread(&y , sizeof(double) , 1 , fpShp) ;
GeoPoint newVertex(x, y);
points.SetPoint(iteratorPoint, newVertex);
}
}
delete []startOfPart ;
startOfPart = 0 ;
pFeature->SetGeometry(pNewGeoLine) ;
this->LoadAttributeData(pFeature,fpDbf,everyRecordLen);
pShpDataSet->AddRow(pFeature) ;
}
break ;
//读取线目标结束

case 5: // '/005'
//读取面目标开始
{
Fields &featureFields = pShpDataSet->GetFields();
Feature *pFeature = new Feature(recordNumber , 1 , &featureFields) ;
int polygonShapeType ;
fread(&polygonShapeType  , sizeof(int) , 1 ,fpShp) ;
if( polygonShapeType != 5 )
AfxMessageBox( "Error: Attempt to load non polygon shape as polygon." ) ;

double objMinX , objMinY , objMaxX , objMaxY ;
fread(&objMinX , sizeof(double) , 1 , fpShp) ;
fread(&objMinY , sizeof(double) , 1 , fpShp) ;
fread(&objMaxX , sizeof(double) , 1 , fpShp) ;
fread(&objMaxY , sizeof(double) , 1 , fpShp) ;

GeoPolygon *pNewGeoPolygon = new GeoPolygon();
double width = objMaxX - objMinX ;
double height = objMaxY - objMinY ;
pFeature->SetBound(objMinX , objMinY , width , height) ;

int numParts , numPoints ;
fread(&numParts , sizeof(int) , 1 , fpShp) ;
fread(&numPoints , sizeof(int) , 1 , fpShp) ;
//存储各个面的起点索引
int* startOfPart = new int[numParts] ;
for( int i = 0 ; i < numParts ; i++ )
{
int indexFirstPoint ;
fread(&indexFirstPoint , sizeof(int) , 1 , fpShp) ;
startOfPart[i] = indexFirstPoint ;
}

//处理单个目标有多面问题
pNewGeoPolygon->SetPointsCount( numParts ) ;

for( i = 0 ; i < numParts ; i++ )
{
GeoPoints& points = pNewGeoPolygon->GetPoints(i) ;
int curPosIndex = startOfPart[i] ;
int nextPosIndex = 0 ;
int curPointCount = 0 ;
if( i == numParts - 1 )
curPointCount = numPoints - curPosIndex ;
else
{
nextPosIndex = startOfPart[i + 1];
curPointCount = nextPosIndex - curPosIndex ;
}
points.SetPointCount( curPointCount ) ;
//加载一个面(多边形)的坐标
for( int iteratorPoint = 0 ; iteratorPoint < curPointCount ; iteratorPoint ++ )
{
double x , y ;
fread(&x , sizeof(double) , 1 , fpShp) ;
fread(&y , sizeof(double) , 1 , fpShp) ;
GeoPoint newVertex(x, y);
points.SetPoint(iteratorPoint, newVertex);
}
}
delete []startOfPart ;
startOfPart = 0 ;
pFeature->SetGeometry(pNewGeoPolygon) ;
this->LoadAttributeData(pFeature,fpDbf,everyRecordLen);
pShpDataSet->AddRow(pFeature) ;
}
break ;
//读取面目标结束

case 23: // '/027'
//读取Measure形线目标开始
{
Fields &featureFields = pShpDataSet->GetFields();
Feature *pFeature = new Feature(recordNumber , 1 , &featureFields) ;
int arcMShapeType ;
fread(&arcMShapeType , sizeof(int) , 1 , fpShp) ;

double objMinX , objMinY , objMaxX , objMaxY ;
fread(&objMinX , sizeof(double) , 1 , fpShp) ;
fread(&objMinY , sizeof(double) , 1 , fpShp) ;
fread(&objMaxX , sizeof(double) , 1 , fpShp) ;
fread(&objMaxY , sizeof(double) , 1 , fpShp) ;

GeoPolyline *pNewGeoLine = new GeoPolyline();
double width = objMaxX - objMinX ;
double height = objMaxY - objMinY ;
pFeature->SetBound(objMinX , objMinY , width , height) ;

int numParts , numPoints ;
fread(&numParts , sizeof(int) , 1 , fpShp) ;
fread(&numPoints , sizeof(int) , 1 , fpShp) ;
//存储各段线的起点索引
int* startOfPart = new int[numParts] ;
for( int i = 0 ; i < numParts ; i++ )
{
int indexFirstPoint ;
fread(&indexFirstPoint , sizeof(int) , 1 , fpShp) ;
startOfPart[i] = indexFirstPoint ;
}

//处理单个目标有多条线的问题
pNewGeoLine->SetPointsCount( numParts ) ;

for( i = 0 ; i < numParts ; i++ )
{
GeoPoints& points = pNewGeoLine->GetPoints(i) ;
int curPosIndex = startOfPart[i] ;
int nextPosIndex = 0 ;
int curPointCount = 0 ;
if( i == numParts - 1 )
curPointCount = numPoints - curPosIndex ;
else
{
nextPosIndex = startOfPart[i + 1] ;
curPointCount = nextPosIndex - curPosIndex ;
}
points.SetPointCount( curPointCount ) ;
//加载一条线段的坐标
for( int iteratorPoint = 0 ; iteratorPoint < curPointCount ; iteratorPoint ++ )
{
double x , y ;
fread(&x , sizeof(double) , 1 , fpShp) ;
fread(&y , sizeof(double) , 1 , fpShp) ;
GeoPoint newVertex(x, y);
points.SetPoint(iteratorPoint, newVertex);
}
}
delete []startOfPart ;
startOfPart = 0 ;

double* value = new double[2 + numPoints] ;
fread( value , sizeof(double) , 2+numPoints, fpShp) ;
delete []value ;
value = 0 ;

pFeature->SetGeometry(pNewGeoLine) ;
this->LoadAttributeData(pFeature,fpDbf,everyRecordLen);
pShpDataSet->AddRow(pFeature);
}
break ;
//读取Measure形线目标结束
}
}
return pShpDataSet ;
}

另：我手头上有一个读取shp文件的一个小程序（C++），如果需要的朋友请留下邮箱。我会第一时间发送给你

转载的朋友请注明出处：半支烟阿杰 http:// href="http://blog.csnd.net/gisfarmer" target=_blank>blog.csnd.net/gisfarmer