1、簡介
1.1 什麼是Shapefile
<font color=blue>ESRI Shapefile（shp），或簡稱shapefile，是美國環境系統研究所公司（ESRI）開發的一種空間資料開放格式。該檔案格式已經成為了地理資訊軟體界的一個開放标準，這表明ESRI公司在全球的地理資訊系統市場的重要性。
- GIS 保留的資料大緻分為栅格資料和矢量資料：
【GIS開發】Esri Shapefile(.shp)矢量資料檔案讀取（C++、Python）1、簡介2、C++3、Python結語 - 而矢量資料檔案主要有如下兩種格式：
【GIS開發】Esri Shapefile(.shp)矢量資料檔案讀取（C++、Python）1、簡介2、C++3、Python結語
<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資料檔案。
【GIS開發】Esri Shapefile(.shp)矢量資料檔案讀取（C++、Python）1、簡介2、C++3、Python結語 - （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'])
- 運作結果如下：
【GIS開發】Esri Shapefile(.shp)矢量資料檔案讀取（C++、Python）1、簡介2、C++3、Python結語
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)