【GIS开发】Esri Shapefile(.shp)矢量数据文件读取（C++、Python）

1、简介

1.1 什么是Shapefile

<font color=blue>ESRI Shapefile（shp），或简称shapefile，是美国环境系统研究所公司（ESRI）开发的一种空间数据开放格式。该文件格式已经成为了地理信息软件界的一个开放标准，这表明ESRI公司在全球的地理信息系统市场的重要性。

• GIS 保留的数据大致分为栅格数据和矢量数据：
• 而矢量数据文件主要有如下两种格式：

<font color=purple>Shapefile属于一种矢量图形格式，它能够保存几何图形的位置及相关属性。用于描述几何体对象：点，折线与多边形。例如，Shapefile文件可以存储井、河流、湖泊等空间对象的几何位置。除了几何位置，shp文件也可以存储这些空间对象的属性，例如一条河流的名字，一个城市的温度等等。

Shapefile 是一种用于存储地理要素的几何位置和属性信息的非拓扑简单格式。 shapefile 中的地理要素可表示为点、线或面（区域）。 包含 shapefile 的工作空间还可以包含 dBASE 表，它们用于存储可连接到 shapefile 的要素的附加属性。

<font color=green>Shapefile文件指的是一种文件存储的方法，实际上该种文件格式是由多个文件组成的。其中，要组成一个Shapefile，有三个文件是必不可少的，它们分别是".shp", ".shx"与 ".dbf"文件。表示同一数据的一组文件其文件名前缀应该相同。例如，存储一个关于湖的几何与属性数据，就必须有lake.shp，lake.shx与lake.dbf三个文件。而其中“真正”的Shapefile的后缀为shp，然而仅有这个文件数据是不完整的，必须要把其他两个附带上才能构成一组完整的地理数据。除了这三个必须的文件以外，还有八个可选的文件，使用它们可以增强空间数据的表达能力。所有的文件名都必须遵循MS DOS的8.3文件名标准（文件前缀名8个字符，后缀名3个字符，如shapefil.shp），以方便与一些老的应用程序保持兼容性，尽管现在许多新的程序都能够支持长文件名。此外，所有的文件都必须位于同一个目录之中。

1.2 格式分析

shapefile文件主要由四个文件组成：

• 必须的文件: (几何与属性是一对一关系，这种关系基于记录编号。dBASE 文件中的属性记录必须与主文件中的记录采用相同的顺序。)

.shp - 用于存储要素几何的主文件；必需文件。
.shx - 用于存储要素几何索引的索引文件；必需文件。
.dbf - 用于存储要素属性信息的 dBASE 表；必需文件。
.prj - 用于存储坐标系信息的文件；由 ArcGIS 使用。

• 其他可选的文件：

.sbn 和 .sbx - 用于存储要素空间索引的文件。
.fbn 和 .fbx - 用于存储只读 shapefile 的要素空间索引的文件。
.ain 和 .aih - 用于存储某个表中或专题属性表中活动字段属性索引的文件。
.atx - .atx 文件针对在 ArcCatalog 中创建的各个 Shapefile 或 dBASE 属性索引而创建。ArcGIS 不使用 shapefile 和 dBASE 文件的 ArcView GIS 3.x 属性索引。已为 shapefile 和 dBASE 文件开发出新的属性索引建立模型。
.ixs - 读/写 shapefile 的地理编码索引。
.mxs - 读/写 shapefile（ODB 格式）的地理编码索引。

.xml - ArcGIS 的元数据 - 用于存储 shapefile 的相关信息。
.cpg - 可选文件，指定用于标识要使用的字符集的代码页。

<font color=grey>在每个.shp,.shx与.dbf文件之中，图形在每个文件的排序是一致的。也就是说，.shp的第一条记录与.shx及.dbf之中的第一条记录相对应，如此类推。此外，在.shp与.shx之中，有许多字段的字节序是不一样的。因此用户在编写读取这些文件格式的程序时，必须十分小心地处理不同文件的不同字节序。

<font color=black>Shapefile通常以X与Y的方式来处理地理坐标，一般X对应经度，Y对应纬度，用户必须注意X，Y的顺序。

<font color=grey>组成 shapefile 的每个文件均被限制为 2 GB。 因此，.dbf 文件不能超过 2 GB，.shp 文件也不能超过 2 GB（只有这两个文件的容量会很大）。 所有组成文件的总大小可以超过 2 GB。

NULL = 0 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

1.3 资源下载

可以从OpenStreetMap下载一些地图数据。 https://www.openstreetmap.org https://download.geofabrik.de/ https://extract.bbbike.org/

2、C++

2.1 shapelib

http://shapelib.maptools.org/ https://github.com/OSGeo/shapelib The Shapefile C Library provides the ability to write simple C programs for reading, writing and updating (to a limited extent) ESRI Shapefiles, and the associated attribute file (.dbf).

If you don't know, you probably don't need this library. The Shapefile format is a working and interchange format promulagated by ESRI for simple vector data with attributes. An excellent white paper on the shapefile format is available from ESRI, but it is .pdf format, so you will need Adobe Acrobat to browse it.

The file format actually consists of three files.

XXX.shp - holds the actual vertices.
XXX.shx - hold index data pointing to the structures in the .shp file.
XXX.dbf - holds the attributes in xBase (dBase) format.

BBBike 是一款免费的OpenStreetMap数据下载服务器！这个服务器能够OpenStreetMap中提取全球200多个地区的不同格式数据。

点击下载全球200多个地区的目录，每个城市的大小为2-50MB。支持的格式有: OSM，PBF，GeoJSON，SQLite，Garmin (style OSM，Cycle，Leisure，Onroad，Openfietslite，OpenSeaMap，opentopmap，BBBike) ，Osmand，mapsforge，Navit，maps.me，SVG，和 Esri Shapefile。

• SHPObject ：shape文件中包含的要素对象

typedef struct
  {
    int		nSHPType;	//Shape Type (SHPT_* - see list above) 要素类型

    int		nShapeId; 	//Shape Number (-1 is unknown/unassigned) ID

    int		nParts;		//# of Parts (0 implies single part with no info) 要素有几个部分组成
    int		*panPartStart;  //Start Vertex of part 要素部分的起始点
    int		*panPartType;	//Part Type (SHPP_RING if not SHPT_MULTIPATCH) 各个部分的类型
    
    int		nVertices;	//Vertex list 
    double	*padfX;		
    double	*padfY;
    double	*padfZ;		//(all zero if not provided)
    double	*padfM;		//(all zero if not provided)

    double	dfXMin;		//Bounds in X, Y, Z and M dimensions
    double	dfYMin;
    double	dfZMin;
    double	dfMMin;

    double	dfXMax;
    double	dfYMax;
    double	dfZMax;
    double	dfMMax;
  } SHPObject;

• （1）在网站BBBike上框选需要的区域，下载对应的shp数据文件。
• （2）解压shp文件如下：

• （3）编写代码（C++、OpenGL、glut）绘制shp数据文件如下：

只绘制buildings.shp: 同时绘制buildings.shp、landuse.shp、railways.shp、waterways.shp等: 部分代码如下：

#include <iostream>
#include <shapefil.h>

int main()
{
	const char* filename = "C:\\Users\\tomcat\\Desktop\\planet_103.95105,1.31994_103.96025,1.32555-shp\\buildings.shp";
	SHPHandle hShp = SHPOpen(filename,"r");

	int nShapeType, nVertices;
	int nEntities = 0;
	double* minB = new double[4];
	double* maxB = new double[4];
	SHPGetInfo(hShp, &nEntities, &nShapeType, minB, maxB);
	printf("ShapeType:%d\n", nShapeType);
	printf("Entities:%d\n", nEntities);

	SHPClose(hShp);
	system("Pause");
}

SHPHandle OpenShapeFile(char* fileName, bool useBoundingBox = true)
{    
    SHPHandle hSHP=SHPOpen(fileName, "rb" );
	if (hSHP == nullptr) {
		return hSHP;
	}
	
	if (useBoundingBox) {
		m_sBoundingBox.fMaxX = hSHP->adBoundsMax[0];
		m_sBoundingBox.fMaxY = hSHP->adBoundsMax[1];
		m_sBoundingBox.fMinX = hSHP->adBoundsMin[0];
		m_sBoundingBox.fMinY = hSHP->adBoundsMin[1];
	}

	int nShapeType, nEntities;
	double adfMinBound[4], adfMaxBound[4];
	SHPGetInfo(hSHP, &nEntities, &nShapeType, adfMinBound, adfMaxBound);

}

再来一组世界地图的shp文件的绘制，如下：

2.2 gdal

ogr库是gdal的一部分，定义并实现了空间参考类OGRSpatialReference，这个库不仅能自定义地图投影参数，还实现了wkt和proj4两个字符串接口，用起来很方便。

• （1）与proj4的接口函数

OGRErr importFromProj4( const char * );//根据proj4字符串初始化
OGRErr exportToProj4( char ** ) const;//输出proj4字符串

• （2）与wkt字符串的接口函数

OGRErr importFromWkt( char ** );//根据现有wkt初始化
OGRErr exportToWkt( char ** ) const;//输出wkt字符串

• （3）常用的用户自定义地图投影函数

/** Universal Transverse Mercator */
OGRErr SetUTM( int nZone, int bNorth = TRUE );

/** Transverse Mercator */
OGRErr      SetTM( double dfCenterLat, double dfCenterLong,
				   double dfScale,
				   double dfFalseEasting, double dfFalseNorthing );

OGRErr      SetGeogCS( const char * pszGeogName,
					   const char * pszDatumName,
					   const char * pszEllipsoidName,
					   double dfSemiMajor, double dfInvFlattening,
					   const char * pszPMName = nullptr,
					   double dfPMOffset = 0.0,
					   const char * pszUnits = nullptr,
					   double dfConvertToRadians = 0.0 );

OGRErr      SetTOWGS84( double, double, double,
						double = 0.0, double = 0.0, double = 0.0,
						double = 0.0 );

/** 或者直接从epsg编码一次性得到所有参数 */
OGRErr      importFromEPSG( int );

https://cpp.hotexamples.com/examples/-/OGRSpatialReference/exportToProj4/cpp-ogrspatialreference-exporttoproj4-method-examples.html

C++读取shapefile中的.prj文件，代码如下：

#include <iostream>
#include <fstream>

#include "ogr_spatialref.h"
#pragma comment(lib, "gdal_i.lib")
   
std::string result;
std::ifstream ifs("C:\\Users\\tomcat\\Desktop\\planet_103.95105,1.31994_103.96025,1.32555-shp\\shape\\buildings.prj", std::ios::in | std::ios::binary);
if (!ifs) {
	std::cout << "error" << std::endl;
	return ;
}

std::string str((std::istreambuf_iterator<char>(ifs)),
std::istreambuf_iterator<char>());

OGRSpatialReference* theSrs = new OGRSpatialReference(str.c_str());
OGRErr err = theSrs->Validate();

char* theProj4;
if (theSrs && theSrs->Validate() == OGRERR_NONE) {
	theSrs->morphFromESRI();
	
	if (theSrs->exportToProj4(&theProj4) == OGRERR_NONE) {
		result = theProj4;
	}
}

输入（.prj文件，wkt字符串）：

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

输出（proj4字符串）：

+proj=tmerc +lat_0=0 +lon_0=114 +k=1 +x_0=38500000 +y_0=0 +ellps=GRS80 +units=m +no_defs

3、Python

https://www.lfd.uci.edu/~gohlke/pythonlibs

3.1 pyshp

https://pypi.org/project/pyshp/

The Python Shapefile Library (PyShp) reads and writes ESRI Shapefiles in pure Python.

pip install pyshp

import os
import shapefile  # 使用pyshp
import matplotlib.pyplot as plt

def plot_polygon(poly,symbol='r-',**kwargs):
    x, y = zip(*poly)
    plt.plot(x,y,symbol,**kwargs)

def plot_layer(filename,symbol,layer_index=0,**kwargs):
    sf = shapefile.Reader(filename)
    shapes = sf.shapes()
    for i in range(len(shapes)):
        shapeType = shapes[i].shapeType
        points = shapes[i].points
        plot_polygon(points,symbol,**kwargs)

os.chdir(r"C:/Users/tomcat/Desktop/planet_103.95105,1.31994_103.96025,1.32555-shp/shape")

plot_layer('buildings.shp','g-',markersize=1)
plot_layer('landuse.shp','bo',markersize=1)
plot_layer('roads.shp','b-',markersize=1)
plot_layer('waterways.shp','r-',markersize=1)

plt.axis('equal')
plt.gca().get_xaxis().set_ticks([])
plt.gca().get_yaxis().set_ticks([])
plt.title("draw by pyshp")
plt.show()

3.2 osgeo

pip install gdal

• 读取shp文件代码如下：

import os
import matplotlib.pyplot as plt
from osgeo import ogr

def plot_point(point,symbol='ko',**kwargs):
    x,y=point.GetX(),point.GetY()
    plt.plot(x,y,symbol,**kwargs)

def plot_line(line,symbol='g-',**kwargs):
    x,y=zip(*line.GetPoints())
    plt.plot(x,y,symbol,**kwargs)

def plot_polygon(poly,symbol='r-',**kwargs):
    for i in range(poly.GetGeometryCount()):
        subgeom=poly.GetGeometryRef(i)
        x,y=zip(*subgeom.GetPoints())
        plt.plot(x,y,symbol,**kwargs)

def plot_layer(filename,symbol,layer_index=0,**kwargs):
    ds=ogr.Open(filename)
    for row in ds.GetLayer(layer_index):
        geom=row.geometry()
        geom_type=geom.GetGeometryType()

        if geom_type==ogr.wkbPoint:
            plot_point(geom,symbol,**kwargs)
        elif geom_type==ogr.wkbMultiPoint:
            for i in range(geom.GetGeometryCount()):
                subgeom=geom.GetGeometryRef(i)
                plot_point(subgeom,symbol,**kwargs)

        elif geom_type==ogr.wkbLineString:
            plot_line(geom,symbol,**kwargs)
        elif geom_type==ogr.wkbMultiLineString:
            for i in range(geom.GetGeometryCount()):
                subgeom=geom.GetGeometryRef(i)
                plot_line(subgeom,symbol,**kwargs)

        elif geom_type == ogr.wkbPolygon:
            plot_polygon(geom,symbol,**kwargs)
        elif geom_type==ogr.wkbMultiPolygon:
            for i in range(geom.GetGeometryCount()):
                subgeom=geom.GetGeometryRef(i)
                plot_polygon(subgeom,symbol,**kwargs)

os.chdir(r"C:/Users/tomcat/Desktop/planet_103.95105,1.31994_103.96025,1.32555-shp/shape")

plot_layer('buildings.shp','g-',markersize=1)
plot_layer('landuse.shp','bo',markersize=1)
plot_layer('roads.shp','b-',markersize=1)
plot_layer('waterways.shp','r-',markersize=1)

plt.axis('equal')
plt.gca().get_xaxis().set_ticks([])
plt.gca().get_yaxis().set_ticks([])
plt.title("draw by gdal")
plt.show()

• 读取.prj文件，代码如下：

#!/usr/bin/env python 

from osgeo import osr
import sys, os

os.chdir(r"C:/Users/tomcat/Desktop/planet_103.95105,1.31994_103.96025,1.32555-shp/shape")

def main(prj_file): 
    prj_text = open(prj_file, 'r').read()
    srs = osr.SpatialReference()
    if srs.ImportFromWkt(prj_text):
        raise ValueError("Error importing PRJ information from: %s" % prj_file) 
    print(srs.ExportToProj4()) 
    # print(srs.ExportToWkt()) 

if __name__=="__main__": 
    # main(sys.argv[1]) 
    main("buildings.prj")

3.3 geopandas

https://geopandas.org/en/stable//

geopandas沿用了pandas的数据类型，也就是具有GeoSeries，GeoDataFrame两种主要的数据结构，继承了pandas数据结构中大部分操作方法。

import os
import geopandas
import matplotlib.pyplot as plt

os.chdir(r"C:/Users/tomcat/Desktop/planet_103.95105,1.31994_103.96025,1.32555-shp/shape")

pic = geopandas.read_file("buildings.shp")
print(pic.crs)      # 查看数据对应的投影信息
print(pic.head())   # 查看前5行数据
pic.plot()
plt.title("draw by geopandas")
plt.show()

3.4 fiona

Fiona中对数据的描述模型和GDAL中的不同： （1）GDAL中对于矢量数据采用数据源（DataSource）- 图层（Layer）- 要素（Feature）- 属性和几何体（Attributes and Geometry）。 （2）Fiona采用Python中内置的数据结构表示矢量数据，一个要素以GeoJSON表示，使用Python内置的字典（dict）结构组织；一个图层包含在一个集合中（Collection）。可以对该集合进行迭代遍历，得到其中的要素。

• 安装库如下：

pip install fiona

• 测试代码如下：

import fiona

with fiona.open(r'C:\Users\tomcat\Desktop\100000_full\100000_full.shp', encoding='utf-8') as c:

    print(f'bounds: {c.bounds}')
    print(f'crs: {c.crs}')
    print(f'driver: {c.driver}')
    print(f'encoding: {c.encoding}')

    fields = c.schema['properties']
    print('properties: ')
    for k, v in fields.items():
        print(f'{k} -> {v}')
  
    for f in c.items():
        print(f[1]['properties']['name'])

• 运行结果如下：

3.5 shapely

shapely是一个python包，用于设置平面特征的理论分析和操作（通过python的 ctypes 模块）来自著名和广泛部署的地理类库的功能。

pip install pyproj
pip install shapely

计算距离和面积，如下:

from pyproj import Geod
from shapely.geometry import Point, LineString, Polygon
 
# 计算距离
line_string = LineString([Point(123.429056, 41.833526), Point(123.435235, 41.773191)])
geod = Geod(ellps="WGS84")
length = geod.geometry_length(line_string)
print(length)

# 计算封闭区域面积、周长
geod = Geod(ellps="WGS84")
poly_area, poly_perimeter = geod.geometry_area_perimeter(Polygon([
    (123.432746, 41.792136), (123.435922, 41.784584),
    (123.450857, 41.784392), (123.447681, 41.79252),
]))
print(poly_area, poly_perimeter)

结语

如果您觉得该方法或代码有一点点用处，可以给作者点个赞，或打赏杯咖啡；╮(￣▽￣)╭ 如果您感觉方法或代码不咋地//(ㄒoㄒ)//，就在评论处留言，作者继续改进；o_O??? 如果您需要相关功能的代码定制化开发，可以留言私信作者；(✿◡‿◡) 感谢各位大佬童鞋们的支持！( ´ ▽´ )ﾉ ( ´ ▽´)っ！！！