PYTHON绘制GIS叠加栅格数据示意图

魔法披风.

已于 2024-06-11 11:55:33 修改

阅读量760

点赞数 16

分类专栏： python 文章标签： python 数据可视化

于 2024-06-06 18:41:17 首次发布

本文链接：https://blog.csdn.net/Charls_Wu/article/details/139503281

版权

python 专栏收录该内容

1 篇文章 0 订阅

订阅专栏

1.需求描述

绘制指定线路所在经纬度范围的当前时刻雷达回波覆盖情况示意图，图中同时叠加或绘制以下信息：

（1）地名

（2）行政区划

（3）线路及途径点

（4）雷达回波

并且在生成的图片中附加以下信息：

（1）标题（描述线路相关信息，时间等）

（2）线路及雷达渐变填色图图例

2.PYTHON依赖

使用PYTHON实现上述需求，使用相关依赖如下：

（1）matplotlib： Python 绘图库，提供了丰富的数据可视化功能

（2）rasterio：各类标准格式的栅格数据读取及处理

（3）geopandas：地理空间数据处理

（4）cartopy：地理空间数据处理，这里主要使用其提供的各类投影定义

（5）contextily：在栅格及矢量数据可视化的基础上叠加底图信息

（6）numpy

3.实现步骤

3.1所需数据：

（1）需要可视化的栅格数据，支持PNG、TIFF等常见格式

（2）省界矢量数据，支持geojson或shp等常见格式

（3）背景地图，可以使用天地图、谷歌等支持xyz方式访问的地图服务，contextily.provider也提供众多可选地图；注意若选择地图网络响应速度慢，则会极大响应绘图时长。

（4）线路矢量数据，geojson或shp等常见格式

3.2绘制步骤

（1）创建图片绘制容器（fig）及图片绘制对象（ax）

# 参考坐标系
projection = cartopy.crs.epsg(3857)
# 创建容器及绘制对象
fig, ax = plt.subplots(figsize=(10, 10), subplot_kw={'projection': projection})

（2）设置地图绘制范围

# 读取航路矢量文件
cfp_route_gpd = geopandas.read_file(cfp_route_geojson, driver='GeoJSON')

#投影转换
cfp_route_gpd = cfp_route_gpd.to_crs(projection)

# 获取地图范围
cfp_route_bbox = cfp_route_gpd.total_bounds
cfp_route_bbox = reshape_map_extent(cfp_route_bbox)

# 设置绘制范围
ax.set_extent([cfp_route_bbox[0], cfp_route_bbox[2], cfp_route_bbox[1], cfp_route_bbox[3]], crs=projection)


def reshape_map_extent(cfp_route_bbox):
    """
    地图范围规范为正方形,并各方向扩展350KM
    """
    expand_distance = 350000
    # 使地图裁剪范围呈正方形
    x_span = cfp_route_bbox[2] - cfp_route_bbox[0]
    y_span = cfp_route_bbox[3] - cfp_route_bbox[1]
    if x_span > y_span:
        expand = (x_span - y_span) / 2
        cfp_route_bbox[3] += expand
        cfp_route_bbox[1] -= expand
    else:
        expand = (y_span - x_span) / 2
        cfp_route_bbox[2] += expand
        cfp_route_bbox[0] -= expand
    cfp_route_bbox[3] += expand_distance
    cfp_route_bbox[1] -= expand_distance
    cfp_route_bbox[2] += expand_distance
    cfp_route_bbox[0] -= expand_distance
    return cfp_route_bbox

本需求中绘制范围为指定线路bbox按照长宽大小关系处理为正方形，再在各方向延申350KM。

（3）绘制Title等附加信息

def plot_information(cfp_route_data_source_param, flight_info_param):
    # 获取当前格林威治时间
    gmt_now = datetime.datetime.utcnow()
    # 格式化时间
    formatted_time = gmt_now.strftime("%Y-%m-%d %H:%M:%S")
    topleft_title_text = f'Time: {formatted_time} ( UTC )'
    plt.title(topleft_title_text, y=1, loc='left', fontdict={
        'fontsize': 10
    })

（4）绘制栅格渐变填色图及其图例

def grid_data_clip(dataset, gdf_geometry_4326):
    """
    根据实际需要展示的地图范围，对数据进行裁剪，以提升数据渲染效率
    本例实测裁剪后绘制时间减少4秒左右
    由于本例中使用PNG作为数据源，不带地理信息数据，则先通过地图展示范围经纬度坐标计算像素坐标，通过像素坐标范围进行数据裁剪
    """
    origin_coverage_extent = [73, 134.99, 12.21, 54.2]
    grid_resolution = 0.01
    y_index_range = (origin_coverage_extent[3] - origin_coverage_extent[2]) / 0.01

    x_min_index = round((gdf_geometry_4326[0].x - origin_coverage_extent[0]) / grid_resolution)
    x_max_index = round((gdf_geometry_4326[1].x - origin_coverage_extent[0]) / grid_resolution)
    y_min_index = round((gdf_geometry_4326[0].y - origin_coverage_extent[2]) / grid_resolution)
    y_max_index = round((gdf_geometry_4326[1].y - origin_coverage_extent[2]) / grid_resolution)

    clip_window = ((y_index_range - y_max_index, y_index_range - y_min_index), (x_min_index, x_max_index))

    data = dataset.read(1, window=clip_window)  # 读取第一个波段的数据
    return data


def plot_radar(ax, map_extend):
    """
    绘制雷达
    """
    # 创建自定义颜色映射
    # 假设你有特定的边界值和对应的颜色
    boundaries = [10, 14, 20, 26, 32, 36, 40, 44, 48, 52, 56, 60, 64, 68]
    colors = [
        (0 / 255, 101 / 255, 154 / 255, 0 / 255),
        (0 / 255, 101 / 255, 154 / 255, 180 / 255),
        (0 / 255, 144 / 255, 147 / 255, 220 / 255),
        (0 / 255, 179 / 255, 125 / 255, 240 / 255),
        (117 / 255, 208 / 255, 89 / 255, 255 / 255),
        (220 / 255, 220 / 255, 30 / 255, 255 / 255),
        (244 / 255, 202 / 255, 8 / 255, 255 / 255),
        (245 / 255, 168 / 255, 24 / 255, 255 / 255),
        (236 / 255, 130 / 255, 63 / 255, 255 / 255),
        (205 / 255, 75 / 255, 75 / 255, 255 / 255),
        (182 / 255, 45 / 255, 100 / 255, 255 / 255),
        (156 / 255, 16 / 255, 109 / 255, 255 / 255),
        (125 / 255, 0 / 255, 108 / 255, 255 / 255),
        (92 / 255, 0 / 255, 100 / 255, 255 / 255)
    ]
    cm = ListedColormap(colors)
    norm = BoundaryNorm(boundaries, len(colors) - 1)

    # 读取地形栅格数据
    # todo 数据来源问题待处理
    dataset = rasterio.open("dbz.png", mode='r')

    web_mercator_extend__points = [
        (map_extend[0], map_extend[1]),
        (map_extend[2], map_extend[3]),
    ]
    gdf_web_mercator = gpd.GeoDataFrame(
        geometry=gpd.points_from_xy([x for x, y in web_mercator_extend__points], [y for x, y in web_mercator_extend__points]),
        crs="EPSG:3857")

    # 将坐标参考系统转换为WGS84 (EPSG:4326)
    gdf_geometry_4326 = gdf_web_mercator.to_crs(epsg=4326).geometry

    data = grid_data_clip(dataset, gdf_geometry_4326)
    # 创建掩膜，去除NoData区域
    data = np.ma.masked_where(data == 254, data)

    # 将栅格添加到地图中
    plot_extent = [gdf_geometry_4326[0].x, gdf_geometry_4326[1].x, gdf_geometry_4326[0].y, gdf_geometry_4326[1].y]
    im = ax.imshow(data, origin="upper", extent=plot_extent
                   , transform=ccrs.PlateCarree(), cmap=cm, norm=norm
                   , resample=True, interpolation='bicubic')

    # 添加图例
    cbar_ax = plt.axes([0.6, 0.08, 0.3, 0.015])  # [right, bottom, width, height]
    plt.colorbar(im, ax=ax, orientation="horizontal", pad=1, aspect=500, shrink=0.3, cax=cbar_ax)
    cbar_ax.text(10.2, 0.2, 'DBZ', fontsize=7)
    dataset.close()

（5）绘制地图上展示的线路图例

def plot_legend():
    route_legend_ax = plt.axes([0.16, 0.08, 0.12, 0.015])  # [right, bottom, width, height]
    y_value = 0
    x = [2.5, 3, 6, 6.5]
    y = [y_value, y_value, y_value, y_value]
    route_legend_ax.text(1, -0.18, 'Route:', fontsize=10)

    x_index = 0
    for x_value in x:
        if x_index < len(x) - 1 and x_index != 0:
            y_value = y[x_index]
            draw_concentric_circles(route_legend_ax, x_value, y_value, n_circles=2, radii=[0.13, 0.27],
                                    edge_colors=[None, (0, 0, 0, 0.8)], colors=[(0, 0, 0, 0.8), None])
        x_index += 1

    route_legend_ax.plot(x, y, '-', linewidth=1.5, color='gray')
    route_legend_ax.axis('off')


def draw_concentric_circles(ax, x, y, n_circles=4, colors=None, edge_colors=None, radii=None):
    """
    自定义绘制同心圆
    """
    if colors is None:
        colors = plt.cm.viridis(np.linspace(0, 1, n_circles))
    if radii is None:
        radii = np.linspace(5, 30, n_circles)  # 控制每个圆的半径
    for radius, color, edge_color in zip(radii, colors, edge_colors):
        circle = plt.Circle((x, y), radius, edgecolor=edge_color, linewidth=1, facecolor=(0, 0, 0, 0), color=color)
        ax.add_patch(circle)

在容器中新增一个子图并指定子图的相对位置的宽高，再将图例信息绘制到子图中，由于是线路图例，需要同时绘制线和途径点，其中途径点是绘制的点加圆环的样式，线路图例效果如下：

（6）绘制省界

# 省界矢量数据
area_gpd = geopandas.read_file("areasBorder.geojson")
area_gpd = area_gpd.to_crs(projection)

# 绘制省界
area_gpd.plot(ax=ax, transform=projection, edgecolor='black', alpha=0.5, linewidth=0.5)

此处是通过geojson文件读取数据。

（7）添加底图

# 添加注记底图
base_map_labels_service_url = 'https://ip:port/{z}/{x}/{y}.png'
contextily.add_basemap(ax
                   , source=base_map_labels_service_url)

指定地图服务并添加到图片中，这里是自定义xyz形式的底图服务地址，也可以直接通过contextily.provider获取底图服务；底图服务响应时间极大的影响图片生成时长。

（8）绘制线路及途径点

def plot_cfp_route(ax, projection, cfp_route_gpd):
    """
    绘制航线及航路点
    """
    cfp_route_gpd.plot(ax=ax, transform=projection, linewidth=1.5, edgecolor='gray')  # linestyle='--'，
    for cpt_point in cfp_route_gpd['geometry'][0].coords:
        draw_concentric_circles(ax, cpt_point[0], cpt_point[1], n_circles=2, radii=[5000, 12000],
                                edge_colors=[None, (0, 0, 0, 0.8)], colors=[(0, 0, 0, 0.8), None])

其中也用到了绘制同心圆的方法draw_concentric_circles。cfp_route_gpd是前文中通过geopandas获取的geojson文件的线路矢量数据

（9）保存图片

plt.savefig("Pic.jpg", dpi=200, bbox_inches="tight", pad_inches=0.1)

其中dpi指定图片像素大小，bbox_inches="tight",和pad_inches=0.1用于调整子图相对于绘图容器的外边距。上述数据或信息的绘制顺序影响其在生成图片中的重叠顺序！！！

4.总结

该需求中，主要注意以下几点：

（1）参考系统一，所有叠加数据坐标需要转换到统一参考系下，本例指定的是WEB墨卡托投影坐标系

（2）图像内容位置控制，图例均通过在容器中创建子图进行展示，容易控制其相对于容器的位置，包括各类文字位置等，其位置和宽高等数值需要根据实际效果进行多次调整。

（3）图像生成时间长的问题，可优化的点包括：

1）背景底图服务的响应时间，使用响应快的底图服务，如果业务需求不必须使用底图，也可以不绘制底图

2）用地名及行政区划矢量数据进行自定义绘制，而非使用地图服务，将大幅提升图片绘制效率

3）控制数据量，进行必要的数据裁剪等

最终图片生成效果如下：