首页 > 编程语言 >【GIS开发】Esri Shapefile(.shp)矢量数据文件读取（C++、Python）

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

时间：2022-10-17 15:33:25浏览次数：74  
标签：shp plot plt GIS 文件 Shapefile double geom

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??? 如果您需要相关功能的代码定制化开发，可以留言私信作者；(✿◡‿◡) 感谢各位大佬童鞋们的支持！( ´ ▽´ )ﾉ ( ´ ▽´)っ！！！

标签：shp,plot,plt,GIS,文件,Shapefile,double,geom
From： https://blog.51cto.com/u_15800063/5760884

相关文章

  • 外向交货单屏幕增强BADI：LE_SHP_TAB_CUST_OVER
    外向交货单屏幕增强1、涉及tcode：VL01N  Vl02N VL03N2、效果图3、后台增强点4、增强点，BADI：LE_SHP_TAB_CUST_OVER5、具体方法介绍：1）ACTIVATE_TAB_PAGE2）TRANSFER_DATA_TO_......
  • 基于Logistic上市公司财务危机预警分析
    财务危机财务危机又称财务困境，是指企业由于营销、决策或不可抗拒因素的影响，使经营循环和财务循环无法正常持续或陷于停滞的状态，具体表现包括持续性亏损、无偿付能力、违约和......
  • postgis 安装错误记录
    postgres单机+postgis安装详细版https://zhuanlan.zhihu.com/p/143690286ERROR:couldnotloadlibrary"/opt/postgresql-9.3.2/lib/postgis-2.3.so":/opt/postgr......
  • 栅格转shp
    importarcpyfromarcpyimportenvenv.workspace="E:/data/temp/data"field="VALUE"inRaster="single_"outPolygons="E:/data/temp/data/shpsingle_"num=3while......
  • dedecms系统php.ini register_globals must is Off的解决方案
    访问后台，报错如下图：成因：由于register_globals设置控制PHP变量访问范围,如果开启会引起不必要的安全问题,所以这里对其进行了强制关闭,如果站长的空间不支持,可以采用以下几......
  • 外向交货单增强BADI：LE_SHP_DELIV…
    "参考销售订单创建外向交货单  CALL FUNCTION 'BAPI_OUTB_DELIVERY_CREATE_SLS'    EXPORTING    ......
  • 外向交货单增强BADI：LE_SHP_DELIVERY_PR
     "参考销售订单创建外向交货单  CALL FUNCTION 'BAPI_OUTB_DELIVERY_CREATE_SLS'    EXPORTING      SHIP_POINT  ......
  • 【GIS开发】osgEarth依赖库PROJ（Python）
    文章目录​​1、OSGeo/PROJ（C++）​​​​1.1编译sqlite3​​​​1.2编译libtiff​​​​1.3编译openssl​​​​1.4编译curl​​​​1.5编译PROJ9​​​​2、pyproj（pytho......
  • ArcGIS蜂巢六边形
    统计区域内的数量时，区域可以使用六边形。ArcMap10.4版本开始，提供了制作六边形的工具——生成细分曲面。本案例中，生成了黑龙江省范围内的细分曲面。Shape中的类型选择六边形......
  • 【GIS开发】OpenLayers入门学习（JavaScript库）
    1、简介官网地址：https://openlayers.org/源码地址：https://github.com/openlayers/openlayersOpenLayers是一个高性能、功能丰富的库，用于在Web上创建交互式地图。它......

赞助商

阅读排行

网站主页关于我们联系我们网站地图

本网站内容转载自其他媒体，侵权联系[admin##ips99.com]。

Copyright © 2020-2023 IPS99 版权所有 IPS99