爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 5794|回复: 26

[混合编程] 气小Py,轻磅来袭!无偿接各类绘图、代码编写~~欢迎供稿~

[复制链接]

新浪微博达人勋

发表于 2023-5-6 05:27:31 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
本帖最后由 edwardli 于 2023-5-11 13:32 编辑

https://www.bilibili.com/video/B ... 999.section.playall
QQ:342205284
要画图?一起呗!~
不会画图?没关系 , 我们一起研究着画
图画成了 = 你的问题解决了  +  我的素材也有了,商量着一起出视频哈。
一举好几得



喜欢研究各种数据,爱好画各类图形,就是不缺钱。
哈哈哈哈哈哈哈哈哈哈哈哈哈哈哈哈哈哈哈哈哈哈。


本帖被以下淘专辑推荐:

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2023-5-6 07:41:52 | 显示全部楼层
已关注  学习
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-5-6 07:57:26 | 显示全部楼层
学习
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2023-5-6 08:13:20 | 显示全部楼层
233333333老李仗义
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-5-6 09:08:57 | 显示全部楼层
{:eb502:}{:eb502:}
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-5-6 10:09:24 | 显示全部楼层
工作中我们可能经常遇到批量绘图的情况,但是python运行时往往只是单线程运行,数据量、绘图量大的话CPU拉不满,绘图就可能比较慢。使用多进程运行就可以充分发挥CPU的多核性能了,这是我的小白代码,主要multiprocessing多线程用于绘制雷达图像,但在分配任务时,会有bug(数据个数不能被设置的进程数整除的话,余数的最后几个数据不会被处理)抛砖引玉,请大佬们多交流。
  1. import numpy as np
  2. import matplotlib
  3. matplotlib.use('agg')
  4. from matplotlib import pyplot as plt
  5. from multiprocessing import Process
  6. import os
  7. import cartopy.crs as ccrs
  8. import cartopy.io.shapereader as shapereader
  9. from pycwr.io import read_auto
  10. from pycwr.draw.RadarPlot import add_rings,GraphMap

  11. def get_filename(path,fileType):
  12.     file_name =[]
  13.     final_name = []
  14.     for files in os.listdir(path):  #root为目录路径 #dirs为路径下的子目录 #files为路径下的所有非子目录
  15.         if fileType in files:
  16.             file_name.append(files.replace(fileType,''))#生成不带‘.bin’后缀的文件名组成的列表
  17.             final_name = [path +item +fileType for item in file_name]#生成‘.bin’后缀的文件名组成的绝对路径列表
  18.             final_name.sort()
  19.     return final_name #输出列表

  20. def generate_heat_array(sa_radar_file,fig_dir,mfl,lon1d,lat1d,proj,extent,name_list):
  21.    
  22.     path1=[]
  23.     path2=[]
  24.     path3=[]
  25.     path4=[]
  26.     path5=[]
  27.     path6=[]
  28.     path7=[]
  29.     path8=[]
  30.     path9=[]
  31.     path10=[]
  32.     num = 0
  33.     while num <= (len(sa_radar_file)/10-1):
  34.         path1.append(sa_radar_file[10*num])
  35.         path2.append(sa_radar_file[10*num+1])
  36.         path3.append(sa_radar_file[10*num+2])
  37.         path4.append(sa_radar_file[10*num+3])
  38.         path5.append(sa_radar_file[10*num+4])
  39.         path6.append(sa_radar_file[10*num+5])
  40.         path7.append(sa_radar_file[10*num+6])
  41.         path8.append(sa_radar_file[10*num+7])
  42.         path9.append(sa_radar_file[10*num+8])
  43.         path10.append(sa_radar_file[10*num+9])
  44.         num = num+1
  45.    
  46.     #用multiprocessing实现多进程
  47.     process1=Process(target=show_heat_map,args=(path1,fig_dir,mfl,lon1d,lat1d,proj,extent,name_list))
  48.     process2=Process(target=show_heat_map,args=(path2,fig_dir,mfl,lon1d,lat1d,proj,extent,name_list))
  49.     process3=Process(target=show_heat_map,args=(path3,fig_dir,mfl,lon1d,lat1d,proj,extent,name_list))
  50.     process4=Process(target=show_heat_map,args=(path4,fig_dir,mfl,lon1d,lat1d,proj,extent,name_list))
  51.     process5=Process(target=show_heat_map,args=(path5,fig_dir,mfl,lon1d,lat1d,proj,extent,name_list))
  52.     process6=Process(target=show_heat_map,args=(path6,fig_dir,mfl,lon1d,lat1d,proj,extent,name_list))
  53.     process7=Process(target=show_heat_map,args=(path7,fig_dir,mfl,lon1d,lat1d,proj,extent,name_list))
  54.     process8=Process(target=show_heat_map,args=(path8,fig_dir,mfl,lon1d,lat1d,proj,extent,name_list))
  55.     process9=Process(target=show_heat_map,args=(path9,fig_dir,mfl,lon1d,lat1d,proj,extent,name_list))
  56.     process10=Process(target=show_heat_map,args=(path10,fig_dir,mfl,lon1d,lat1d,proj,extent,name_list))
  57.     process1.start()
  58.     process2.start()
  59.     process3.start()
  60.     process4.start()
  61.     process5.start()
  62.     process6.start()
  63.     process7.start()
  64.     process8.start()
  65.     process9.start()
  66.     process10.start()
  67.     process1.join()
  68.     process2.join()
  69.     process3.join()
  70.     process4.join()
  71.     process5.join()
  72.     process6.join()
  73.     process7.join()
  74.     process8.join()
  75.     process9.join()
  76.     process10.join()

  77. def show_heat_map(path,fig_dir,mfl,lon1d,lat1d,proj,extent,name_list):
  78.     # for key,color in zip(['dBZ'],['CN_ref']):#'V','CN_vel',
  79.     for key in ['dBZ','V']:
  80.         for level in [0,1,2,3,4,5]:#
  81.             print(level,'begin')
  82.             for nFiles in range(len(path)):
  83.                 # print(os.path.basename(path[nFiles]),'begin')
  84.                 PRD = read_auto(path[nFiles]) #新版本数据
  85.                 PRD.add_product_CR_lonlat(XLon=lon1d, YLat=lat1d)
  86.                 # XLon:np.ndarray, 1d, units:degrees
  87.                 # YLat:np.ndarray, 1d, units:degrees
  88.                 # level_height:常量,要计算的高度 units:meters
  89.                 fig = plt.figure(figsize=(6,6))
  90.                 ax = plt.axes(projection=proj , zorder = 5)
  91.                 # ax = fig.subplots(1, 1, subplot_kw={'projection': proj})
  92.                 ax.set_extent(extent, crs=proj)
  93.                 graph = GraphMap(PRD, proj)
  94.                 graph.plot_ppi_map(ax=ax, sweep_num=level, extend=extent,field_name=key) ## ,min_max=[-35,35],cmap = 'CN_vel' 0代表第一层, dBZ代表反射率产品,cmap
  95.                 # plot_lonlat_map(ax, grid_lon, grid_lat, PRD.product.CR_geo, transform=proj,)### bounds= np.arange(25,76,5),cmap = 'Greys'
  96.                 ### 设置经纬度间隔 修改Anaconda3\Lib\site-packages\pycwr\draw\RadarPlot.py
  97.                 ### 中def plot_lonlat_map 参数 parallels meridians
  98.                 add_rings(ax,rings=[50,100],color="red", linestyle='-', linewidth=0.6,)
  99.                 # ax.add_geometries(shapereader.Reader(mfl1).geometries(), crs=proj, facecolor='none',edgecolor='k',linewidth=1)
  100.                 for state, geos in zip(shapereader.Reader(mfl).records(), shapereader.Reader(mfl).geometries()):
  101.                     if state.attributes['NAME'] in name_list:
  102.                         ax.add_geometries([geos], crs=proj,facecolor='None', edgecolor='red',zorder = 5)  # 重点区域上色
  103.                 # start_time = str(PRD.scan_info.start_time.values + np.timedelta64(8, 'h'))[5:19]
  104.                 end_time = str(PRD.scan_info.end_time.values + np.timedelta64(8, 'h'))[5:19]
  105.                 ax.set_title(end_time+' Z9856 '+str(level)+key, fontsize=16)

  106.    
  107.                 plt.tight_layout()
  108.                 plt.savefig(fig_dir+os.path.basename(path[nFiles])+'_'+str(level)+'_'+key+'.png',dpi = 200)
  109.                 plt.clf()  
  110.                 # plt.close()
  111.                 print(level,end_time,' Z9856 done')
  112.             print(level,'done')
  113.    
  114. if __name__=='__main__':
  115.     #pdb.set_trace()
  116. #"""读取目标文件夹下的多个雷达基数据"""
  117.     Dir = "G:\\data\\9856\\" #目标文件夹
  118.     fig_dir = "D:\\9856_fig_all\\"
  119.     fileType = '.bz2' #雷达基数据后缀名
  120.     sa_radar_file = get_filename(Dir,fileType)
  121.     # print(sa_radar_file)
  122. #"""----------------------------------------------------------------------"""
  123.     mfl = 'G:/代码/地图文件/micaps_shp_utf8/County.shp'
  124.     mfl1 = 'G:/代码/地图文件/micaps_shp_utf8/City.shp'
  125.     name_list = ['碧江区','万山区','江口县','玉屏侗族自治县','石阡县','思南县',
  126.     '印江土家族苗族自治县','德江县','沿河土家族自治县','松桃苗族自治县',]
  127.     # extent = [108, 110, 27, 29]#限定绘图范围[108, 110, 27, 29]
  128.     extent = [107.5, 109.5, 27, 29.5]
  129.     proj = ccrs.PlateCarree()#创建投影,选择cartopy的platecarree投影
  130.     lon1d = np.arange(extent[0], extent[1], 0.01) ##lon方向0.01等间距,lon_min-1on_max范围
  131.     lat1d = np.arange(extent[2], extent[3], 0.01) ##lat方向0.01等间距,lat_min-lat_max度范围
  132.     # grid_lon, grid_lat = np.meshgrid(lon1d, lat1d, indexing="ij")
  133.     # for key,color in zip(['dBZ','ZDR','CC','PhiDP','KDP'],['CN_ref','pyart_RefDiff','pyart_RefDiff','pyart_RefDiff','pyart_RefDiff']):
  134.     # return mfl,mfl1,extent,proj,lon1d,lat1d
  135.     generate_heat_array(sa_radar_file,fig_dir,mfl,lon1d,lat1d,proj,extent,name_list)
  136.     print('finish')
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-5-6 10:16:17 | 显示全部楼层
格局很大!!!{:eb502:}{:eb502:}
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-5-6 10:45:00 | 显示全部楼层
迅速关注一波
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-5-6 11:11:07 | 显示全部楼层
李sir 这格局牛逼了。
那我等必须要不耻下问了。


密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2023-5-7 11:49:31 | 显示全部楼层
1163087557 发表于 2023-5-6 10:09
工作中我们可能经常遇到批量绘图的情况,但是python运行时往往只是单线程运行,数据量、绘图量大 ...
  1. import numpy as np
  2. import matplotlib

  3. matplotlib.use('agg')
  4. from matplotlib import pyplot as plt
  5. import multiprocessing
  6. import os
  7. import cartopy.crs as ccrs
  8. import cartopy.io.shapereader as shapereader
  9. from pycwr.io import read_auto
  10. from pycwr.draw.RadarPlot import add_rings, GraphMap


  11. def get_filename(path, fileType):
  12.     file_name = []
  13.     final_name = []
  14.     for files in os.listdir(path):  # root为目录路径 #dirs为路径下的子目录 #files为路径下的所有非子目录
  15.         if fileType in files:
  16.             file_name.append(files.replace(fileType, ''))  # 生成不带‘.bin’后缀的文件名组成的列表
  17.             final_name = [path + item + fileType for item in file_name]  # 生成‘.bin’后缀的文件名组成的绝对路径列表
  18.             final_name.sort()
  19.     return final_name  # 输出列表


  20. def show_heat_map(file, fig_dir, mfl, lon1d, lat1d, proj, extent, name_list):
  21.     # for key,color in zip(['dBZ'],['CN_ref']):#'V','CN_vel',
  22.     for key in ['dBZ', 'V']:
  23.         for level in [0, 1, 2, 3, 4, 5]:  #
  24.             print(level, 'begin')
  25.             # print(os.path.basename(path[nFiles]),'begin')
  26.             PRD = read_auto(file)  # 新版本数据
  27.             PRD.add_product_CR_lonlat(XLon=lon1d, YLat=lat1d)
  28.             # XLon:np.ndarray, 1d, units:degrees
  29.             # YLat:np.ndarray, 1d, units:degrees
  30.             # level_height:常量,要计算的高度 units:meters
  31.             fig = plt.figure(figsize=(6, 6))
  32.             ax = plt.axes(projection=proj, zorder=5)
  33.             # ax = fig.subplots(1, 1, subplot_kw={'projection': proj})
  34.             ax.set_extent(extent, crs=proj)
  35.             graph = GraphMap(PRD, proj)
  36.             graph.plot_ppi_map(ax=ax, sweep_num=level, extend=extent, field_name=key)  ## ,min_max=[-35,35],cmap = 'CN_vel' 0代表第一层, dBZ代表反射率产品,cmap
  37.             # plot_lonlat_map(ax, grid_lon, grid_lat, PRD.product.CR_geo, transform=proj,)### bounds= np.arange(25,76,5),cmap = 'Greys'
  38.             ### 设置经纬度间隔 修改Anaconda3\Lib\site-packages\pycwr\draw\RadarPlot.py
  39.             ### 中def plot_lonlat_map 参数 parallels meridians
  40.             add_rings(ax, rings=[50, 100], color="red", linestyle='-', linewidth=0.6, )
  41.             # ax.add_geometries(shapereader.Reader(mfl1).geometries(), crs=proj, facecolor='none',edgecolor='k',linewidth=1)
  42.             for state, geos in zip(shapereader.Reader(mfl).records(), shapereader.Reader(mfl).geometries()):
  43.                 if state.attributes['NAME'] in name_list:
  44.                     ax.add_geometries([geos], crs=proj, facecolor='None', edgecolor='red', zorder=5)  # 重点区域上色
  45.             # start_time = str(PRD.scan_info.start_time.values + np.timedelta64(8, 'h'))[5:19]
  46.             end_time = str(PRD.scan_info.end_time.values + np.timedelta64(8, 'h'))[5:19]
  47.             ax.set_title(end_time + ' Z9856 ' + str(level) + key, fontsize=16)

  48.             plt.tight_layout()
  49.             plt.savefig(fig_dir + os.path.basename(file) + '_' + str(level) + '_' + key + '.png', dpi=200)
  50.             plt.clf()
  51.             # plt.close()
  52.             print(level, end_time, ' Z9856 done')
  53.             print(level, 'done')


  54. def generate_heat_array(sa_radar_file, fig_dir, mfl, lon1d, lat1d, proj, extent, name_list):
  55.     pools = multiprocessing.Pool(processes=10)
  56.     for file in sa_radar_file:
  57.         pools.apply_async(show_heat_map, (file, fig_dir, mfl, lon1d, lat1d, proj, extent, name_list))
  58.     pools.close()
  59.     pools.join()


  60. if __name__ == '__main__':
  61.     # pdb.set_trace()
  62.     # """读取目标文件夹下的多个雷达基数据"""
  63.     Dir = "G:\\data\\9856\"  # 目标文件夹
  64.     fig_dir = "D:\\9856_fig_all\"
  65.     fileType = '.bz2'  # 雷达基数据后缀名
  66.     sa_radar_file = get_filename(Dir, fileType)
  67.     # print(sa_radar_file)
  68.     # """----------------------------------------------------------------------"""
  69.     mfl = 'G:/代码/地图文件/micaps_shp_utf8/County.shp'
  70.     mfl1 = 'G:/代码/地图文件/micaps_shp_utf8/City.shp'
  71.     name_list = ['碧江区', '万山区', '江口县', '玉屏侗族自治县', '石阡县', '思南县',
  72.                  '印江土家族苗族自治县', '德江县', '沿河土家族自治县', '松桃苗族自治县', ]
  73.     # extent = [108, 110, 27, 29]#限定绘图范围[108, 110, 27, 29]
  74.     extent = [107.5, 109.5, 27, 29.5]
  75.     proj = ccrs.PlateCarree()  # 创建投影,选择cartopy的platecarree投影
  76.     lon1d = np.arange(extent[0], extent[1], 0.01)  ##lon方向0.01等间距,lon_min-1on_max范围
  77.     lat1d = np.arange(extent[2], extent[3], 0.01)  ##lat方向0.01等间距,lat_min-lat_max度范围
  78.     # grid_lon, grid_lat = np.meshgrid(lon1d, lat1d, indexing="ij")
  79.     # for key,color in zip(['dBZ','ZDR','CC','PhiDP','KDP'],['CN_ref','pyart_RefDiff','pyart_RefDiff','pyart_RefDiff','pyart_RefDiff']):
  80.     # return mfl,mfl1,extent,proj,lon1d,lat1d
  81.     generate_heat_array(sa_radar_file, fig_dir, mfl, lon1d, lat1d, proj, extent, name_list)
  82.     print('finish')
复制代码

密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表