矢量文件的读取和写入

ytkz2024-01-252024-09-09

矢量的本质是什么

矢量由点、线、面矢量组成，点、线、面是地理信息系统（GIS）中的基本元素，它们分别表示地图上的位置、路径和区域。

点（Point）：表示具有地理位置的单一位置，例如城市或地标。
线（Line）：表示两个或多个点之间的路径，例如街道或河流。
面（Polygon）：表示封闭的区域，例如湖泊或国家。

Shapefile（.shp）是一种流行的地理数据格式，用于存储地理空间信息。一个Shapefile实际上是由几个文件组成的，其中主要的有：

.shp：包含几何对象的主文件。
.shx：包含几何对象索引的索引文件。
.dbf：包含属性信息的dBase表格。

这些文件必须具有相同的前缀，并且必须存储在同一目录下。此外，还可能有其他类型的文件，例如.prj文件（包含投影信息）和.sbn/.sbx文件（包含空间索引信息）等。

怎么查看

关于是否可以使用文本编辑器（如vscode）打开和读取Shapefile，答案是不能。Shapefile文件包含的数据是以二进制形式存储的，如果直接用文本编辑器打开，你会看到一堆看似无意义的字符，而不是可读的地理数据。为了读取和操作Shapefile文件，你需要使用专门的GIS软件（如QGIS、ArcGIS）或者编程库（如Python的GDAL、geopandas库）。

写矢量文件

#!/usr/bin/env python
# -*- coding: utf-8 -*- 
from osgeo import ogr, osr

# 创建WGS84空间参考系统
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)  # WGS84的EPSG编号是4326

# 创建数据源
driver = ogr.GetDriverByName('ESRI Shapefile')
datasource = driver.CreateDataSource('output.shp')

# 创建点
layer = datasource.CreateLayer('point_layer', srs, geom_type=ogr.wkbPoint)
feature = ogr.Feature(layer.GetLayerDefn())
point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(10, 20)  # 假设这是WGS84经纬度
feature.SetGeometry(point)
layer.CreateFeature(feature)

# 创建线
layer = datasource.CreateLayer('line_layer', srs, geom_type=ogr.wkbLineString)
feature = ogr.Feature(layer.GetLayerDefn())
line = ogr.Geometry(ogr.wkbLineString)
line.AddPoint(10, 20)  # 假设这是WGS84经纬度
line.AddPoint(30, 40)  # 假设这是WGS84经纬度
feature.SetGeometry(line)
layer.CreateFeature(feature)

# 创建面
layer = datasource.CreateLayer('polygon_layer', srs, geom_type=ogr.wkbPolygon)
feature = ogr.Feature(layer.GetLayerDefn())
ring = ogr.Geometry(ogr.wkbLinearRing)
ring.AddPoint(10, 20)  # 假设这是WGS84经纬度
ring.AddPoint(10, 40)  # 假设这是WGS84经纬度
ring.AddPoint(30, 40)  # 假设这是WGS84经纬度
ring.AddPoint(30, 20)  # 假设这是WGS84经纬度
ring.AddPoint(10, 20)  # 假设这是WGS84经纬度
polygon = ogr.Geometry(ogr.wkbPolygon)
polygon.AddGeometry(ring)
feature.SetGeometry(polygon)
layer.CreateFeature(feature)

# 释放资源
feature = None
datasource = None

在实际应用中，你需要使用实际的经纬度坐标。

读取矢量

要使用GDAL库读取Shapefile文件并进行可视化，你需要使用GDAL进行读取，然后使用matplotlib进行可视化。以下是一个示例：

首先，确保已经安装了gdal和matplotlib库。

用gdal读取

#!/usr/bin/env python
# -*- coding: utf-8 -*- 
# @Time : 
# @File : read_shp.py
from osgeo import ogr
import matplotlib.pyplot as plt
# 打开Shapefile文件
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(r'D:\shape\output1.shp', 0)  # 0 means read-only. 1 means writeable.

# 检查文件是否打开成功
if dataSource is None:
    print('Could not open file')
else:
    print('Opened file successfully')

# 获取图层
layer = dataSource.GetLayer()

# 遍历图层中的所有要素
for feature in layer:
    geom = feature.GetGeometryRef()

    # 创建x和y的列表，用于存储坐标
    x = []
    y = []

    # 如果要素是一个点
    if geom.GetGeometryName() == 'POINT':
        x.append(geom.GetX())
        y.append(geom.GetY())
    # 如果要素是一个线或多边形
    elif geom.GetGeometryName() in ['LINESTRING', 'POLYGON']:
        for i in range(0, geom.GetGeometryRef(0).GetPointCount()):
            # 获取点的x和y坐标
            x.append(geom.GetGeometryRef(0).GetX(i))
            y.append(geom.GetGeometryRef(0).GetY(i))

    # 使用matplotlib创建线图
    plt.plot(x, y)

# 显示图形
plt.show()

上面的代码把shp可视化结果如下：

在arcgis可视化结果如下：

其他方式

Python 中常用的库来读取 Shapefile (.shp) 文件包括 geopandas 和 fiona。下面是两个例子：

使用 geopandas:

import geopandas as gpd

# 读取 Shapefile
gdf = gpd.read_file('path_to_your_file.shp')

# 打印 GeoDataFrame
print(gdf)

使用 fiona:

import fiona

# 打开 Shapefile
with fiona.open('path_to_your_file.shp', 'r') as shp:
    # 打印每个特征
    for feature in shp:
        print(feature)

在这两个例子中，你需要将 'path_to_your_file.shp' 替换为你的 Shapefile 文件的实际路径。

geopandas 会将 Shapefile 读取为 GeoDataFrame，这是一个可以方便地进行地理数据操作的数据结构。

fiona 会将 Shapefile 读取为一个特征的列表，每个特征都是一个字典，包含了几何形状和属性等信息。

需要注意的是，这两个库都需要你的 Python 环境中安装了相关的依赖，如 GDAL。