Python绘制降雨等值线图、Python绘图

hihihan

花了3天简单研究了一下Python绘制等值线图，效果出来了图片如下：


同样数据surfer制作图如下：


感觉没那么理想，不如surfer做出来的那么好看，可能是插件算法的问题，surfer是用的克里金法（Kriging）,效果要好很多。Python我用的是 Rbf提供的几种插值，试过了，都不是很理想。

地图边界文件我用的是自己用在Surfer中的做的边界，参考文章http://bbs.06climate.com/forum.php?mod=viewthread&tid=54930
站点坐标经过百度地图转换与边界相匹配。边界文件、数据文件如下图


3天学习时间，一天多用来安装各种环境，2天时间从网上找相关代码，进行改造。主要解决了以下主要问题

1、读取数据、地图文件
2、等值线样式调试
3、地图边界加载
4、地图外白话处理
主要代码如下：

import numpy as np
import pandas as pd

from scipy.interpolate import Rbf  # 径向基函数 ： 将站点信息插到格点上 用于绘制等值线

import  matplotlib.pyplot as plt
import matplotlib.colors as  colors
import  matplotlib as mpl

from descartes import PolygonPatch
from shapely.geometry import Polygon

import cartopy.crs as ccrs
import  cartopy.io.shapereader as shpreader
from cartopy.mpl.ticker import  LongitudeFormatter, LatitudeFormatter

import boundary
import maskout


datafile ="G:/pythonProject/anshanyl.txt"
blnfile ="G:/pythonProject/anshanbln.txt"
df = pd.read_csv(datafile, encoding="utf-8")
df.columns=["lon","lat","rain","name","num","namerain"]
print(df)
lon = df["lon"]
lat= df["lat"]
rain = df["rain"]

olon=np.linspace(122,124,30)#设置网格经度
olat=np.linspace(40,41.8,30)#设置网格纬度

olon,olat=np.meshgrid(olon,olat)#网格化

func=Rbf(lon,lat,rain,function='linear')#定义径向基函数插值
rain_data_new = func(olon,olat) #插值
rain_data_new[rain_data_new <0 ] = 0
print(rain_data_new)

fig= plt.figure(figsize=(10,9),facecolor="#666666",edgecolor="Blue",frameon=False, dpi=200)
ax = fig.add_subplot(111,projection=ccrs.PlateCarree())

clevs = [0.1,10.,25.,50.,100.,250.,500]
cdict = ['#A9F090','#40B73F','#63B7FF','#0000FE','#FF00FC','#850042']

my_cmap = colors.ListedColormap(cdict)
norm = mpl.colors.BoundaryNorm(clevs,my_cmap.N)

cf = ax.contourf(olon,olat,rain_data_new,clevs,transform = ccrs.PlateCarree(),cmap=my_cmap,norm = norm)

ccf =plt.contour(olon,olat,rain_data_new, 2, colors='k',linestyles=['--'],linewidths = 0.5)#显示等值线
plt.clabel(ccf, inline=True, colors=['k', 'r', 'g', 'b'], fmt='%1.2f')#显示等值线数值
#clabel = ax.clabel(ct,fmt="%i")

pos = fig.add_axes([0.87,0.15,0.02,0.2])
plt.colorbar(cf,cax=pos)
pos.set_yticklabels((0,10,25,50,100,250,500))

ax.xaxis.set_major_formatter(LongitudeFormatter(zero_direction_label=True))
ax.yaxis.set_major_formatter(LatitudeFormatter())
ax.set_xticks(np.arange(122,124,0.5),crs = ccrs.PlateCarree())    # x轴
ax.set_yticks(np.arange(40,41.8,0.5),crs = ccrs.PlateCarree())      # y轴
#ax.gridlines()#显示网格


lines = maskout.getBJpolygon(blnfile)
ax  = boundary.addBoundary(blnfile,ax)

ext = [(122, 40), (122, 41.8), (124, 41.8), (124, 40), (122, 40)]
polygon2 = Polygon(ext,[lines])
ax.add_patch(PolygonPatch(polygon2, fc='white', ec='white', alpha=1, zorder=5 ))


plt.savefig('111.png')
plt.show()


待研究问题：1、加载各个站点信息、雨量到地图显示（散点图+pyplot的.text功能应该可以实现，待验证）2、加载地图标题等信息（很简单）

经验总结仅供参考，欢迎交流，主要文件如下：

anshanbln.txt (12.62 KB, 下载次数: 25)

2021-8-31 11:01 上传

点击文件名下载附件


boundary.py (1.65 KB, 下载次数: 23)

2021-8-31 11:02 上传

点击文件名下载附件


main.py (2.47 KB, 下载次数: 26)

2021-8-31 11:02 上传

点击文件名下载附件 maskout.py (3.01 KB, 下载次数: 28)

2021-8-31 11:02 上传

点击文件名下载附件

付亚男

这和工具没什么关系，还是代码的问题

puck66

付亚男 发表于 2021-8-31 12:13
这和工具没什么关系，还是代码的问题

你这白化方案，边缘处理的更细节啊，共享一下喽

龍_彳亍_天_下

付亚男 发表于 2021-8-31 12:13
这和工具没什么关系，还是代码的问题

这就和中央台官网的图片一模一样啊，求个脚本呀

myj0219

感谢分享！！！

feiji158

插值插的细一点呗，0.0000000001度，丝般顺滑

Rainch

收藏学习了

xiaoxiao啊

付亚男 发表于 2021-8-31 12:13
这和工具没什么关系，还是代码的问题

大佬，求一份脚本

枫枫

学习了，刚好在研究画站点降水资料，谢谢

xxx321

请问一下maskout.py是要放到哪个路径下呀，读不出来？