Python 读取.shp文件并生成图幅编号

news2025/1/11 9:04:16

代码适用于需要处理和分析地理空间数据的场景，如城市规划、环境监测或自然资源管理，其中它可以帮助用户读取特定区域的Shapefile文件，确定其地理边界，并基于这些边界计算出按照经纬度5度间隔的图幅编号，进而用于地图制作、空间数据管理和快速数据检索。

实现思路：

代码主要用于处理地理空间数据，具体功能如下：

读取一个Shapefile（.shp）文件，获取其地理边界范围。
根据获取的边界范围，计算出按照经纬度 5° 间隔的图幅编号。
将每个国家的地图边界范围和对应的图幅编号存储起来，并打印出来。

完整代码

import csv
import geopandas as gpd


def get_map_extent(shp_file_path):
    """
    读取.shp文件并获取其边界范围。

    :param shp_file_path: .shp文件路径
    :return: 地图边界的范围 (minx, miny, maxx, maxy)
    """
    # 读取.shp文件
    gdf = gpd.read_file(shp_file_path)

    # 获取地图边界
    boundary = gdf.geometry.total_bounds

    # 返回地图范围
    return boundary


def calculate_sheet_ranges(map_extent):
    """
    根据地图边界计算图幅范围，并按照南北纬和东西经的5度分隔规则返回图幅范围的格式。

    :param map_extent: 地图边界的范围 (minx, miny, maxx, maxy)
    :return: 图幅范围列表，每个元素包含(lon_start, lat_start, lon_end, lat_end, sheet_range)
    """
    minx, miny, maxx, maxy = map_extent

    # 初始化图幅范围列表
    sheet_ranges = []

    # 计算图幅编号
    # for lon_start in range(int(minx // 5) * 5, int(maxx // 5 + 1) * 5, 5):
    #     for lat_start in range(int(miny // 5) * 5, int(maxy // 5 + 1) * 5, 5):
    for lon_start in range(int(minx // 5) * 5, int(maxx // 5 + 1) * 5, 5):
        for lat_start in range(int(miny // 5 + 1) * 5, int(maxy // 5 + 2) * 5, 5):  # GLC_FCS30D 纬度需要上调一格5度
            lon_end = lon_start + 5
            lat_end = lat_start + 5
            sheet_range = f"{'W' if lon_start < 0 else 'E'}{abs(lon_start)}" \
                          f"{'S' if lat_start < 0 else 'N'}{abs(lat_start)}"
            sheet_ranges.append(sheet_range)

    return sheet_ranges


def main():
    iso_data = []
    with open('../../ISO.csv', newline='') as csvfile:
        # 创建 CSV 读取器
        csv_reader = csv.reader(csvfile)
        data = list(csv_reader)  # 将文件内容读取到内存中

        # 初始化行数计数器
        iso_count = 0

        for country in data:
            iso_count += 1
            sids_country = country[0]
            iso_map_extent = [iso_count, sids_country]  # [1, 'BMU']

            # 示例文件路径
            shp_file_path = fr'path_to_admin_division\{sids_country}.shp'

            # 获取地图范围
            map_extent = get_map_extent(shp_file_path=shp_file_path)

            # 计算图幅范围
            sheet_ranges = calculate_sheet_ranges(map_extent=map_extent)

            # 打印结果
            print(iso_count, sids_country, map_extent)

            iso_map_extent.extend(sheet_ranges)
            iso_data.append(iso_map_extent)

            print('- Map Extent Num: ', len(sheet_ranges), iso_map_extent[2:], '\n')


if __name__ == '__main__':
    main()

实现流程：

导入库：导入csv库用于读取CSV文件，导入geopandas库用于处理地理空间数据。
定义get_map_extent函数：该函数接收一个.shp文件路径作为参数，读取该文件，并返回其地理边界范围。
定义calculate_sheet_ranges函数：该函数接收一个边界范围作为参数，计算并返回一个图幅编号列表。图幅编号的计算规则是按照经纬度5度的间隔。
定义main函数：这是程序的主要执行函数。它首先读取一个CSV文件，然后对于CSV文件中的每一个条目，执行以下步骤：
- 构造.shp文件路径。
- 调用get_map_extent函数获取地图边界范围。
- 调用calculate_sheet_ranges函数计算图幅编号。
- 打印出当前处理的国家和其地图边界范围。
- 将地图边界范围和图幅编号添加到结果列表中。
执行入口：如果该脚本作为主程序运行，将调用main函数。

代码运行示例：

预备Shapefile文件：

程序输出：

应用范围：

地理信息系统（GIS）：用于处理和分析地理空间数据。
地图制作：确定地图的边界范围和图幅编号，有助于地图的制作和布局。
空间数据分析：在进行空间数据分析时，可以利用此脚本快速获取地图边界和图幅编号，进而进行更深入的分析。
科研和教育：在地理学、城市规划、环境科学等领域的研究和教学中，用于数据处理和分析。

注意事项：

代码中的shp_file_path是一个格式化字符串，需要根据实际的文件路径进行调整。
calculate_sheet_ranges函数中的注释掉的代码块可能是用于不同的计算规则或示例。
代码中的打印语句用于调试和验证，展示了每个国家的边界范围和图幅编号。

相关函数及调用库的解释说明：

函数及调用库	描述
csv	Python标准库中的模块，用于读写CSV（逗号分隔值）文件。在此代码中，它被用来读取包含国家信息的CSV文件。
geopandas	一个开源项目，扩展了pandas库的功能，使其能够处理地理空间数据。它允许用户读取、处理和分析地理空间数据。
get_map_extent(shp_file_path)	此函数接收一个Shapefile的文件路径作为参数，使用geopandas的`read_file`方法读取Shapefile，并获取其几何边界。返回值是一个包含最小经度、最小纬度、最大经度、最大纬度的元组。
calculate_sheet_ranges(map_extent)	此函数接收一个地图边界范围作为参数，并基于这个范围计算图幅编号。按照南北纬和东西经的5度间隔规则来划分图幅，并返回一个包含这些图幅编号的列表。