通过python创建矢量文件

现有四个点：(1, 1), (2, 2), (3, 3), (4, 4)

以这四个点围起来就是一个面。

如何通过python创建矢量文件。

我们以创建一个面矢量文件为例子，进行阐释。

我们可以使用geopandas、gdal库完成矢量创建。

geopandas

假设我们创建的矢量文件格式为 shp格式，代码如下：

import geopandas as gpd
from shapely.geometry import Polygon
def write_points_as_polygons_to_shp(points, output_file):
    """
    将一系列点坐标转换为以这些点为中心的小正方形多边形，并写入Shapefile文件。

    参数:
    - points: 一个包含点坐标的列表，每个点是一个元组或列表，如[(1, 1), (2, 2), ...]
    - output_file: 输出的Shapefile文件名
    返回:
    - 无，但会生成一个Shapefile文件
    """

    # 创建多边形对象
    simplify_poly_geoms = [Polygon(points)]
    gdf_poly = gpd.GeoDataFrame(geometry=simplify_poly_geoms, crs=CRS("epsg:4326"))
    # 保存为 Shapefile
    gdf_poly.to_file(output_file, driver='ESRI Shapefile')
    
if __name__ == '__main__':
    # 示例用法
    points = [(1, 1), (3, 0), (3, 3), (0, 4)]
    output_file = r'd://temp/create_shp_by_geopandas.shp'
    write_points_as_polygons_to_shp(points, output_file)

矢量在Qgis打开如下：

image-20240827095721260

gdf.to_file() 实际上是在背后调用 Fiona 来读取和写入 Shapefile 文件。

常见的矢量文件格式除了shp之外，还有geojson、GPKG。其实，只要对上述代码中的driver进行修改就可以保存为其他的shp格式。

JSON矢量创建如下：

import geopandas as gpd
from shapely.geometry import Polygon
def write_points_as_polygons_to_geojson(points, output_file):
    """
    将一系列点坐标转换为以这些点为中心的小正方形多边形，并写入Shapefile文件。

    参数:
    - points: 一个包含点坐标的列表，每个点是一个元组或列表，如[(1, 1), (2, 2), ...]
    - output_file: 输出的Shapefile文件名
    返回:
    - 无，但会生成一个Shapefile文件
    """

    # 创建多边形对象
    simplify_poly_geoms = [Polygon(points)]
    gdf_poly = gpd.GeoDataFrame(geometry=simplify_poly_geoms, crs=CRS("epsg:4326"))
    # 保存为 Shapefile
    gdf_poly.to_file(output_file, driver='GeoJSON')
    
if __name__ == '__main__':
    # 示例用法
    points =[(1, 1), (3, 0), (3, 3), (0, 4)]
    output_file = r'd://temp/create_shp_by_geopandas.geojson'
    write_points_as_polygons_to_geojson(points, output_file)

矢量在Qgis打开如下：

image-20240827095645848

gpkg矢量创建如下：

import geopandas as gpd
from shapely.geometry import Polygon
def write_points_as_polygons_to_gpkg(points, output_file):
    """
    将一系列点坐标转换为以这些点为中心的小正方形多边形，并写入Shapefile文件。

    参数:
    - points: 一个包含点坐标的列表，每个点是一个元组或列表，如[(1, 1), (2, 2), ...]
    - output_file: 输出的Shapefile文件名
    - crs_epsg: 坐标参考系统的EPSG代码，默认为4326（WGS 84）
    - square_size: 正方形的边长，默认为0.1单位

    返回:
    - 无，但会生成一个Shapefile文件
    """

    # 创建多边形对象
    simplify_poly_geoms = [Polygon(points)]
    gdf_poly = gpd.GeoDataFrame(geometry=simplify_poly_geoms, crs=CRS("epsg:4326"))
    # 保存为 Shapefile
    gdf_poly.to_file(output_file, driver='GPKG')
    
if __name__ == '__main__':    
        # 示例用法
    points = [(1, 1), (3, 0), (3, 3), (0, 4)]
    output_file = r'd://temp/create_shp_by_geopandas.gpkg'
    write_points_as_polygons_to_gpkg(points, output_file)

矢量在Qgis打开如下：

image-20240827095703990

geopandas很强大也很好用，但是它有个严重的问题，比如说是在程序打包的过程中往往出现打包错误，现不再展开，因为这个问题属于另一个故事了。

fiona

上面说到geopandas使用了fiona进行对矢量文件的读写，那么我们完全可以直接使用fiona读写矢量文件。

import fiona
from shapely.geometry import Polygon
def create_polygon_shapefile_by_fiona(coords, output_path, crs='EPSG:4326', filetype='ESRI Shapefile'):
    """
    创建一个包含单个多边形的 Shapefile 数据集。

    参数:
    - output_path: 输出 Shapefile 的文件路径。
    - coords: 一个列表，包含多边形的坐标，格式为 [(x1, y1), (x2, y2), ..., (xn, yn)]。
    - crs: 字符串，表示坐标参考系统，默认为 WGS84 (EPSG:4326)。
    """
    
    # 创建 Polygon 几何对象
    polygon = Polygon(coords)

    # 设置输出文件的模式、驱动器、模式、CRS 和模式
    schema = {
        'geometry': 'Polygon'}

    # 创建新的 Shapefile 数据集
    with fiona.open(
        output_path,
        mode='w',
        driver= filetype,
        schema=schema,
        crs=crs
    ) as dst:
        # 将 Polygon 写入 Shapefile
        try:

                dst.write({"geometry": polygon.__geo_interface__})
        except Exception as e:
            print(e)

    print(f"Polygon data has been saved to {output_path}")

if __name__ == '__main__':
    # 示例用法
    points = [(1, 1), (3, 0), (3, 3), (0, 4)]
    output_file = r'd://temp/create_shp_by_fiona.shp'
    create_polygon_shapefile_by_fiona(points, output_file)

image-20240827144527271

上面的函数中，我们把filetype设置为默认参数，这个函数还可以生成geojson、gpkg格式的矢量文件。

完整例子如下：

def create_polygon_shapefile_by_fiona(coords, output_path, crs='EPSG:4326', filetype='ESRI Shapefile'):
    """
    创建一个包含单个多边形的 Shapefile 数据集。

    参数:
    - output_path: 输出 Shapefile 的文件路径。
    - coords: 一个列表，包含多边形的坐标，格式为 [(x1, y1), (x2, y2), ..., (xn, yn)]。
    - crs: 字符串，表示坐标参考系统，默认为 WGS84 (EPSG:4326)。
    """
    from shapely.geometry import mapping
    # 创建 Polygon 几何对象
    polygon = Polygon(coords)

    # 设置输出文件的模式、驱动器、模式、CRS 和模式
    schema = {
        'geometry': 'Polygon'}

    # 创建新的 Shapefile 数据集
    with fiona.open(
        output_path,
        mode='w',
        driver= filetype,
        schema=schema,
        crs=crs
    ) as dst:
        # 将 Polygon 写入 Shapefile
        try:

                dst.write({"geometry": polygon.__geo_interface__})
        except Exception as e:
            print(e)

    print(f"Polygon data has been saved to {output_path}")


if __name__ == '__main__':
    # 示例用法
    points = [(1, 1), (3, 0), (3, 3), (0, 4)]
    output_file = r'd://temp/create_shp_by_fiona.shp'
    create_polygon_shapefile_by_fiona(points, output_file)

    output_json_file = r'd://temp/create_shp_by_fiona.geojson'
    create_polygon_shapefile_by_fiona(points, output_json_file, filetype= "GeoJSON")

    output_gpkg_file = r'd://temp/create_shp_by_fiona.gpkg'
    create_polygon_shapefile_by_fiona(points, output_gpkg_file, filetype= "GPKG")

osgeo

GDAL（Geospatial Data Abstraction Library，地理空间数据抽象库）是一个在GIS（地理信息系统）领域广泛使用的开源库，它主要用于读取和写入栅格地理空间数据格式。GDAL提供了一种统一的数据模型来处理这些格式的数据，支持多种主流的栅格数据格式，如GeoTIFF、JPEG、PNG等。

OGR（OpenGIS Simple Features Reference Implementation）是GDAL的一个子项目，专门用于处理矢量地理空间数据格式。OGR支持多种矢量数据格式，包括ESRI Shapefiles、GeoJSON、KML等，使得用户可以在不同的矢量数据格式之间进行转换和处理。

OSR模块允许用户查询、创建和修改空间参考系统，这对于确保地理空间数据在不同坐标系统之间的一致性和准确性至关重要。例如，在进行地理空间数据的投影转换或重采样时，就需要正确地设置和处理空间参考系统。

简单来说，GDAL主要针对栅格地理空间数据，OGR针对矢量地理空间数据，而OSR在GDAL/OGR中则负责空间参考系统的处理。这三者共同构成了强大的地理空间数据处理和分析工具集。

假设我们创建的矢量文件格式为 shp格式，代码如下：


from osgeo import ogr, osr


def create_polygon_shp(longitude_latitude_list, filename):
    driver = ogr.GetDriverByName("ESRI Shapefile")
    data_source = driver.CreateDataSource(filename)

    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)

    # 创建一个多边形层
    layer = data_source.CreateLayer("Polygon", srs, ogr.wkbPolygon)
    field_name = ogr.FieldDefn("data", ogr.OFTString)
    field_name.SetWidth(14)
    layer.CreateField(field_name)

    # 确保列表中的点数量至少为3，以形成一个多边形
    if len(longitude_latitude_list) < 3:
        raise ValueError("需要至少3个点来形成一个多边形")

        # 创建一个多边形几何对象
    polygon = ogr.Geometry(ogr.wkbPolygon)
    # 创建一个线性环（外环）
    ring = ogr.Geometry(ogr.wkbLinearRing)
    for point in longitude_latitude_list:
        # 向线性环中添加点
        ring.AddPoint(float(point[0]), float(point[1]))
        # 将线性环添加到多边形中
    polygon.AddGeometry(ring)

    # 创建一个要素，并设置其几何形状和属性
    feature = ogr.Feature(layer.GetLayerDefn())
    feature.SetGeometry(polygon)
    feature.SetField("data", "Example Polygon")

    # 将要素添加到层中
    layer.CreateFeature(feature)

    # 清理资源
    feature = None
    data_source = None

    print(f"Shapefile {filename} has been created with a polygon feature.")

if __name__ == '__main__':
    output_file = r'd://temp/create_shp_by_ogr.shp'
    create_point_shp(points, output_file)

矢量在Qgis打开如下：

image-20240827100252891

如果要生成其他格式的矢量，修改上述代码这一行：

driver = ogr.GetDriverByName("ESRI Shapefile")

把ESRI Shapefile 改成你想要的格式，其格式不仅限于GPKG、JSON、KML等等

ogr支持的矢量格式共有82种：

image-20240827101423689

0 ESRIC
1 FITS
2 PCIDSK
3 netCDF
4 PDS4
5 VICAR
6 JP2OpenJPEG
7 JPEG2000
8 PDF
9 MBTiles
10 BAG
11 EEDA
12 OGCAPI
13 DB2ODBC
14 ESRI Shapefile
15 MapInfo File
16 UK .NTF
17 LVBAG
18 OGR_SDTS
19 S57
20 DGN
21 OGR_VRT
22 REC
23 Memory
24 CSV
25 GML
26 GPX
27 LIBKML
28 KML
29 GeoJSON
30 GeoJSONSeq
31 ESRIJSON
32 TopoJSON
33 OGR_GMT
34 GPKG
35 SQLite
36 ODBC
37 WAsP
38 PGeo
39 MSSQLSpatial
40 PostgreSQL
41 OpenFileGDB
42 DXF
43 CAD
44 FlatGeobuf
45 Geoconcept
46 GeoRSS
47 GPSTrackMaker
48 VFK
49 PGDUMP
50 OSM
51 GPSBabel
52 OGR_PDS
53 WFS
54 OAPIF
55 Geomedia
56 EDIGEO
57 SVG
58 CouchDB
59 Cloudant
60 Idrisi
61 ARCGEN
62 XLS
63 ODS
64 XLSX
65 Elasticsearch
66 Walk
67 Carto
68 AmigoCloud
69 SXF
70 Selafin
71 JML
72 PLSCENES
73 CSW
74 VDV
75 MVT
76 NGW
77 MapML
78 TIGER
79 AVCBin
80 AVCE00
81 HTTP

小结

以上例子都是生成wgs-84坐标系的矢量，如果要生成其他坐标系的矢量需要更改EPSG的数值。

同时，以上例子都是单个面矢量的创建，至于是多个面矢量、多个孔洞面矢量该如何去创建，这里不再细说。

可以参考这篇文件。