利用 GeoPandas 和xarray(rioxarray)库对中国PM2.5 浓度数据开展城市尺度的分区统计，得到我国所有地级市的平均PM2.5浓度水平。
相关库简介
rioxarray 是一个用于处理地理空间数据的 Python 库，它扩展了 xarray 的功能，使得 xarray 可以方便地处理带有地理参考的多维数组数据（如 NetCDF、GeoTIFF 等）。rioxarray专注于处理地理空间数据的坐标参考系统（CRS）、投影、栅格数据的裁剪、掩膜、重采样等操作，是地理数据分析和遥感数据处理的有力工具。
代码实现
导入相关库
# -*- coding: UTF-8 -*-
# 导入相关库
import pandas as pd
import geopandas as gpd
import numpy as np
import matplotlib.pyplot as plt
from shapely.geometry import Point
import xarray as xr
from shapely.geometry import box
from shapely.geometry import mapping
import rioxarray
加载数据
长时序无缝高分辨率空气污染物浓度数据集（简称LGHAP）对环境管理和地球系统科学分析具有重要意义。在当前发布的 LGHAP 数据集 （LGHAP v2） 中，提供了长达 22 年的无缝气溶胶光学深度 （AOD），近地表 PM2.5 和PM10浓度格网数据，空间分辨率为1km，覆盖中国陆地。详细数据介绍和获取方式见-zenodo。
读取中国区域的PM2.5浓度数据、裁减区域：
#读取全球PM25数据，其以nc文件存储
pm25_data = xr.open_dataset('LGHAP.Global_PM25.Y001.A2021.nc')
# 提取 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'])
PM2.5浓度分区统计
#编写自定义的分区统计函数
def stas_PM25(geometry,nc):
#创建城市几何对象的外接矩形
min_lon, min_lat, max_lon, max_lat = geometry.bounds
bbox = box(min_lon, min_lat, max_lon, max_lat)
#根据城市的几何边界对原始nc文件进行规则裁剪
clip_nc=nc.sel(lon=slice(min_lon, max_lon), lat=slice(min_lat,max_lat))
## 向NC文件中写入WGS84坐标系信息
clip_nc.rio.write_crs("epsg:4326", inplace=True)
# 指定NC文件的xy识别名称
clip_nc.rio.set_spatial_dims(x_dim="lon", y_dim="lat", inplace=True)
#把Geoseries对象转换为GeoDataframe对象
city_shp=gpd.GeoDataFrame(geometry=[geometry])
#对nc进行不规则掩膜
ds=clip_nc.rio.clip(city_shp.geometry.apply(mapping),city_shp.crs,drop=False)
#返回当前城市掩膜下的PM2.5均值，使用nanmean函数来忽略空值nan
return np.nanmean(ds['PM25'].values)
全国所有地级市PM2.5浓度分区统计
#定义用于统计的PM2.5格点数据
nc=china_pm25_data
#使用apply函数对各城市对象进行分区统计，返回对应的PM2.5浓度
china_gdf['PM25']=china_gdf.apply(lambda city:stas_PM25(city['geometry'],nc),axis=1)
结果展示
分区结果可视化
推荐阅读
欢迎关注我的公众号“AI拾贝”，原创技术文章第一时间推送。后台发送pm2.5，自动回复源码和数据。
标签：地级市,min,nc,lon,PM2.5,可视化,lat,china From： https://blog.csdn.net/m0_38065162/article/details/141226377