【附代码+详细注释】用IDW方法插值站点数据并绘制组图
用途&效果
该程序对珠海市的臭氧多年季节平均值进行计算绘图，实现了指北针添加，shpfile添加+mask白化，subplots多个子图制作排版，IDW插值算法等功能，效果如下，代码附后。图片
实现代码
"""
@author: licaoming
@contact: 1639603805@qq.com
@software: jupyter notebook
@file: draw_zhuhai_ato_station_O3_season
@time: 2024/1/6
@function:提升20231226版本画图功能，画季节组图
@purpose:画站点数据图，插值方法IDW
"""
import pandas as pd
import numpy as np
import matplotlib.patches as mpatches
from math import radians, cos, sin, asin, sqrt
import geopandas as gpd
from shapely.geometry import Point, Polygon
import matplotlib.pyplot as plt
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import cartopy.feature as cfeat
import shapefile
from matplotlib.path import Path
#from cartopy.mpl.clip_path
from matplotlib.patches import PathPatch
import cartopy.crs as ccrs
import gma
from matplotlib.pyplot import savefig
#import cartopy as
plt.rcParams['font.sans-serif']=['SimHei']
#读取站点数据
station_data=pd.read_excel("F:/课题/珠海/20231223/station_data_2.xlsx",sheet_name="臭氧季")#季节平均
#station_data=station_data.loc=="臭氧"]#二氧化氮，臭氧
station_data
def haversine(lon1, lat1, lon2, lat2):
# IDW插值方法
R =6372.8
dLon = radians(lon2 - lon1)
dLat = radians(lat2 - lat1)
lat1 = radians(lat1)
lat2 = radians(lat2)
a = sin(dLat/2)**2 + cos(lat1)*cos(lat2)*sin(dLon/2)**2
c = 2*asin(sqrt(a))
d = R * c
return d
def IDW(x, y, z, xi, yi):
lstxyzi = []
for p in range(len(xi)):
lstdist = []
for s in range(len(x)):
d = (haversine(x, y, xi, yi))
lstdist.append(d)
sumsup = list((1 / np.power(lstdist, 2)))
suminf = np.sum(sumsup)
sumsup = np.sum(np.array(sumsup) * np.array(z))
u = sumsup / suminf
xyzi = , yi, u]
lstxyzi.append(xyzi)
return(lstxyzi)
def add_north(ax, labelsize=18, loc_x=0.1, loc_y=0.93, width=0.05, height=0.08, pad=0.1):
"""
画一个比例尺带'N'文字注释
主要参数如下
:param ax: 要画的坐标区域 Axes实例 plt.gca()获取即可
:param labelsize: 显示'N'文字的大小
:param loc_x: 以文字下部为中心的占整个ax横向比例
:param loc_y: 以文字下部为中心的占整个ax纵向比例
:param width: 指南针占ax比例宽度
:param height: 指南针占ax比例高度
:param pad: 文字符号占ax比例间隙
:return: None
"""
minx, maxx = ax.get_xlim()
miny, maxy = ax.get_ylim()
ylen = maxy - miny
xlen = maxx - minx
left =
right =
top =
center = + (top - left)*.4]
triangle = mpatches.Polygon(, color='k')
ax.text(s='N',
x=minx + xlen*loc_x,
y=miny + ylen*(loc_y - pad + height),
fontsize=labelsize,
horizontalalignment='center',
verticalalignment='bottom')
ax.add_patch(triangle)
def maskout_areas(originfig, ax, shp_path, proj=None):
# mask白化
sf = shapefile.Reader(shp_path,encoding='gbk')
vertices = []
codes = []
for shape_rec in sf.shapeRecords():
if shape_rec.record == '斗门区' or shape_rec.record == '金湾区' or shape_rec.record == '香洲区':
pts = shape_rec.shape.points
prt = list(shape_rec.shape.parts) +
for i in range(len(prt) - 1):
for j in range(prt, prt):
if proj:
vertices.append(proj.transform_point(pts, pts, ccrs.Geodetic()))
else:
vertices.append((pts, pts))
codes +=
codes += * (prt - prt - 2)
codes +=
clip = Path(vertices, codes)
clip = PathPatch(clip, transform=ax.transData)
for contour in originfig.collections:
contour.set_clip_path(clip)
#创建经纬格点
entent =
grid = 20
Lon1 = entent
Lon2 = entent
Lat1 = entent
Lat2 = entent
grid_lon_list=np.linspace(Lon1,Lon2,grid)
grid_lat_list=np.linspace(Lat1,Lat2,grid)
# 转换成网格
newlon, newlat = np.meshgrid(grid_lon_list, grid_lat_list)
#数据扁平化
grid_lon_list = newlon.flatten().tolist()
grid_lat_list = newlat.flatten().tolist()
#组图
#画季节平均值图
datayear =["春季","夏季","秋季","冬季"]
latitude =[]
longitude=[]
idw_data=[]
for dy in datayear:
#对站点数据进行IDW插值
lon=station_data["经度"].tolist()
lat=station_data["纬度"].tolist()
data = station_data
sta_data=np.array(IDW(lon,lat,data,grid_lon_list,grid_lat_list))
sta_lon=sta_data[:,0].reshape(20,20)
sta_lat=sta_data[:,1].reshape(20,20)
sta_idw=sta_data[:,2].reshape(20,20)
latitude.append(sta_lat)
longitude.append(sta_lon)
idw_data.append({dy:sta_idw})
season_data=.values()]
#画图
num=0
entent =
shp_path = r"F:/shp/zhuhaimap/zhuhaimap/zhuhai_new.shp"
#建立画布
crs=ccrs.PlateCarree()
fig,axes =plt.subplots(nrows=2,ncols=2,figsize=(30,30),sharex=True,sharey=True,subplot_kw={'projection': crs},constrained_layout=True,dpi=100)
for x, row in enumerate(axes):
for y,col in enumerate(row):
#数据
season_lon=longitude
season_lat=latitude
season_data=.values()]
season=.keys()]
#设置图片范围
axes.set_extent(entent, crs=crs)
#绘制数据图
clevs = np.arange(66,128,2)
cf = axes.contourf(season_lon,season_lat, season_data,clevs,cmap="turbo",tranform=crs)#"nipy_spectral"
#裁剪
maskout_areas(originfig=cf, ax=axes, shp_path=shp_path,proj=None)
#添加矢量地图
shp = gpd.read_file(shp_path)
zh = cfeat.ShapelyFeature(shp.geometry,ccrs.PlateCarree(), edgecolor='k',facecolor='none',)
axes.add_feature(zh, linewidth=0.6, zorder=3)
#画站点
for sta in range(len(station_data["站点"])-1):
axes.scatter(x=station_data["经度"],y=station_data["纬度"], c='k', marker='x')
axes.text(x=station_data["经度"],y=station_data["纬度"],s=station_data["站点"],fontsize=18,ha="right",va="bottom")
#画指北针
add_north(ax=axes)
#设置网格线图例和标题
axes[-1,y].set_xticks(np.arange(entent,113.7, 0.2),crs=ccrs.PlateCarree()) # 这两行画经度，范围为[-180,180]间隔为10
axes.set_yticks(np.arange(21.8,entent,0.2),crs=ccrs.PlateCarree())# 这两行画纬度，范围为[-90,90]间隔为10
axes[-1,y].xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=False))
axes.yaxis.set_major_formatter(LatitudeFormatter())
axes.tick_params(axis="both",labelsize=30)#设置x，y坐标经纬度字体大小
#设置刻度线向外
plt.rcParams['xtick.direction']='out'
plt.rcParams['ytick.direction']='out'
#设置legend图例
cax = fig.add_axes()
cbar = fig.colorbar(cf,ax=axes,orientation='horizontal',cax=cax)
cbar.ax.tick_params(which="major",direction="in",length=2,labelsize=30)#设置图例格式
cbar.set_label("单位：ug/m3",loc="right",fontdict = {'size':30})
#设置teset
season = season
axes.text(x=113.60,y=22.46,s=season,fontsize=30)
#axes.text(x=113.52,y=21.87,s="单位：ug/m3",fontsize='x-large')
num+=1
#调整间距
fig.subplots_adjust(hspace=0.1,wspace=0.1)
plt.show()
文章来源于微信公众号：自学气象人
能分享一下这个项目的更多细节吗？
页:
[1]