【可视化】基于Cartopy绘制中国PM2.5浓度分布填色图

晨昏三变

已于 2024-08-11 15:44:16 修改

阅读量891

点赞数 24

分类专栏：可视化文章标签： python matplotlib 信息可视化

于 2024-08-11 15:33:24 首次发布

本文链接：https://blog.csdn.net/m0_38065162/article/details/141106354

版权

可视化专栏收录该内容

7 篇文章 0 订阅

订阅专栏

【可视化】基于Cartopy绘制中国PM2.5浓度分布填色图

结合Geopandas和Cartopy两个强大的地理数据处理库，对中国的PM2.5浓度进行可视化。通过一个简单案例，帮助大家熟悉Python绘制地图的方法。其中利用Geopandas处理和分析地理空间数据，并利用Cartopy生成平面填色地图，最终将展示中国全年平均PM2.5浓度的空间分布。

代码实现

导入相关库

# -*- coding: UTF-8 -*-
# 导入相关库
import xarray as xr
import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from shapely.geometry import box
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter

加载数据

长时序无缝高分辨率空气污染物浓度数据集（简称LGHAP）对环境管理和地球系统科学分析具有重要意义。在当前发布的 LGHAP 数据集（LGHAP v2）中，提供了长达 22 年的无缝气溶胶光学深度（AOD），近地表 PM2.5 和PM10浓度格网数据，空间分辨率为1km，覆盖中国陆地。详细数据介绍和获取方式见zenodo。

读取全球年均1 km无缝PM2.5网格数据集(LGHAP v2)

#读取全球PM25数据，其以nc文件存储
pm25_data = xr.open_dataset('LGHAP.Global_PM25.Y001.A2021.nc',engine ='netcdf4')

pm25_data.head(10)
# 提取 PM2.5 浓度值和坐标
pm25 = pm25_data['PM25']  
lat = pm25_data['lat']  
lon = pm25_data['lon']

裁剪中国区划范围

#读取中国城市区划shp
china_gdf=gpd.read_file('市.shp')
# 使用 unary_union 将中国的多个行政区划合并为一个
china_shape = china_gdf.geometry.union_all()
#根据行政区划对象创建裁剪范围
min_lon, min_lat, max_lon, max_lat = china_shape.bounds
bbox = box(min_lon, min_lat, max_lon, max_lat)
# 裁剪 PM2.5 数据
china_pm25_data = pm25_data.sel(lon=slice(min_lon, max_lon), lat=slice(min_lat,max_lat))
china_pm25_values = china_pm25_data['PM25']
lon2d, lat2d = np.meshgrid(china_pm25_data['lon'], china_pm25_data['lat'])

读取中国的行政区划和国界线shp

#读取国界线
china_boundary = gpd.read_file('国界线_WGS84.shp') 
#读取行政区划
city_boundry=gpd.read_file('市.shp')

基于Cartopy绘制中国PM2.5浓度分布填色图

#新建绘图框，clear=True在每次新建绘图时删除前置图，节省内存
fig = plt.figure(figsize=(15,15),dpi=300,clear=True)
#创建地图
ax = fig.add_subplot(111,projection=ccrs.PlateCarree())
ax.set_extent([min_lon, max_lon, min_lat, max_lat],crs=ccrs.PlateCarree())
# 添加河流，湖泊
ax.add_feature(cfeature.RIVERS.with_scale('50m'))
ax.add_feature(cfeature.LAKES.with_scale('50m'))
#绘制行政区划
ax.add_geometries(city_boundry["geometry"],crs=ccrs.PlateCarree(),
                fc="None",ec="black",linewidth=.3)
#绘制国界线
ax.add_geometries(china_boundary["geometry"],crs=ccrs.PlateCarree(),
                fc="None",ec="black",linewidth=.7)
#绘制PM2.5浓度
cf=ax.contourf(np.transpose(lon2d), np.transpose(lat2d),china_pm25_values,levels=np.arange(0, 80), cmap='Spectral_r',transform=ccrs.PlateCarree())
#添加colorbar
plt.colorbar(cf,ax=ax, extend='both',orientation='horizontal',label='PM2.5(μg/m3)',pad=0.05) 
#绘制南海子图
sub_ax = fig.add_axes([0.72, 0.3, 0.15, 0.15],
                      projection=ccrs.PlateCarree())
#设置子图范围
sub_ax.set_extent([105,125,0,25], crs=ccrs.PlateCarree())
#绘制子图
sub_ax.add_geometries(city_boundry["geometry"],crs=ccrs.PlateCarree(),
                fc="None",ec="black",linewidth=.3)
sub_ax.add_geometries(china_boundary["geometry"],crs=ccrs.PlateCarree(),
                fc="None",ec="black",linewidth=.5)
sub_ax.contourf(np.transpose(lon2d), np.transpose(lat2d),china_pm25_values,levels=np.arange(0, 80), cmap='Spectral_r',transform=ccrs.PlateCarree())
#添加经纬度刻度
x_extent = [ 80,90,100,110,120,130]
y_extent = [ 5,15,25,35,45]
ax.set_xticks(x_extent, crs=ccrs.PlateCarree())
ax.set_yticks(y_extent, crs=ccrs.PlateCarree())

# 利用Formatter格式化刻度标签
lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
#展示结果
plt.show()