爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
楼主: lwqqqqq

[求助] Python处理云量信息

[复制链接]

新浪微博达人勋

发表于 2018-12-9 20:47:21 | 显示全部楼层
期待作者的经验分享
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2018-12-10 12:39:01 | 显示全部楼层
可以把代码贴上来看看,让大家给分析分析~~
另外,请问楼主的ECMWF数据是在哪里下载的?免费的吗?能否提供一下网址?
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2018-12-10 16:28:33 | 显示全部楼层
chongzika 发表于 2018-12-9 10:38
回头楼主可以分享下经验

尴尬...就用的for循环,处理一个月的数据大概需要20分钟 很慢
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2018-12-10 16:31:15 | 显示全部楼层
香草拿铁 发表于 2018-12-9 20:47
期待作者的经验分享

没有经验分享 用的for循环,尝试了一次多线程,还不如单线程处理呢
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2018-12-10 16:34:47 | 显示全部楼层
junyang517 发表于 2018-12-10 12:39
可以把代码贴上来看看,让大家给分析分析~~
另外,请问楼主的ECMWF数据是在哪里下载的?免费的吗?能否提 ...

对的 是免费的 网址如下:https://www.ecmwf.int/forecasts/datasets/archive-datasets/browse-reanalysis-datasets
代码稍后贴上
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2018-12-10 16:59:08 | 显示全部楼层
import threading
from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import time


def data_process(month):
    nc_obj = Dataset("E:\\data\\2015_" + month + ".nc")
   
    month_dict={'Jan':1,'Feb':2,'Mar':3,'Apr':4,'May':5,'Jun':6,
                'Jul':7,'Aug':8,'Sep':9,'Oct':10,'Nov':11,'Dec':12}
   
    v = [month_dict[k] for k in month_dict.keys() if k==month][0]
    list1 = [1,3,5,7,8,10,12]
    list2 = [2]
    list3 = [4,6,9,11]
        
    lats_list = range(53)
   
    lons_list = range(81)
   
    if v in list1:
        times_list = range(8*31)
    if v in list2:
        times_list = range(8*28)
    if v in list3:
        times_list = range(8*30)
        
    arr1 = np.zeros( (53,81), dtype=np.float64 )

    for lat in lats_list:
            for lon in lons_list:
                for t in times_list:
                    hcc = nc_obj.variables["hcc"][t,lat,lon]
                    lcc = nc_obj.variables["lcc"][t,lat,lon]
                    mcc = nc_obj.variables["mcc"][t,lat,lon]
                    count = 0
                    if lcc <= 0.05 and mcc <= 0.05 and hcc >= 0.05:
                        count += 1
                arr1[lat][lon] = count/len(times_list)


    m=Basemap(projection='tmerc',lat_0 = 23, lon_0 = 118,llcrnrlat=16.5,
              llcrnrlon=108,urcrnrlat=23,urcrnrlon=118)
#加载中国地图
    m.readshapefile('E:\\China_map\\region','region.shp',color='red')
   
#生成绘图背景格点场
    ny=arr1.shape[0];nx=arr1.shape[1]
    lons,lats=m.makegrid(nx,ny)
    x,y=m(lons,lats)
    yy=y[::-1]
   
#画经纬度网格
    m.drawparallels(np.arange(23,16.5,-1.5),labels=[1,0,0,0])
    m.drawmeridians(np.arange(108,118,1.5),labels=[0,0,0,1])
   
#画填色图
    cm = plt.cm.get_cmap('Blues')    #读取一个内置的色系
    shades=m.contourf(x,yy,arr1,cmap=cm)  #将数据填色绘出
    m.colorbar(shades)                      #添加colorbar

#画等值线
#curve=m.contour(x,yy,arr1,colors='black')
#plt.clabel(curve,fmt='%1.0f',colors='black')
   
#添加标题
    plt.title('2015_%d_hcc'%v)
   
#保存
    plt.savefig('2015_hcc\\%d月.png'%v,dpi=300)
#关闭绘图功能
    plt.close()


   
#单线程处理
#time1 = time.time()
#data_process("Jun")
#time2 = time.time()
#print(time2 - time1)
#print("\a")


#多线程处理
threads = []
t1 = threading.Thread(target=data_process,args=('Sep',))
threads.append(t1)
t2 = threading.Thread(target=data_process,args=('Nov',))
threads.append(t2)

time1 = time.time()

if __name__ == '__main__':
    for t in threads:
        t.setDaemon(True)
        t.start()
    t.join()
    print("\a")

time2 = time.time()
print(time2 - time1)
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2018-12-10 17:01:43 | 显示全部楼层
lwqqqqq 发表于 2018-12-10 16:34
对的 是免费的 网址如下:https://www.ecmwf.int/forecasts/datasets/archive-datasets/browse-reanalysi ...

代码已发,麻烦给点意见,看看能不能优化下
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2018-12-12 11:28:16 | 显示全部楼层
tmp = np.zeros((len(times_list),53,81), dtype=np.float64)
hcc = nc_obj.variables["hcc"]   # 不清楚你的文件结构,这部分酌情替换为nc数据读入一个三维数组
lcc = nc_obj.variables["lcc"]
mcc = nc_obj.variables["mcc"]
tmp[(lcc <= 0.05)&(mcc <= 0.05)&(hcc >= 0.05)] = 1
arr1 = np.mean(tmp,0)

以上替换for的部分,我没有测试,仅供参考
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2018-12-12 14:52:26 | 显示全部楼层
暗藏 发表于 2018-12-12 11:28
tmp = np.zeros((len(times_list),53,81), dtype=np.float64)
hcc = nc_obj.variables["hcc"]   # 不清楚 ...

好的 我来试一下 多谢了
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2018-12-13 10:18:36 | 显示全部楼层
暗藏 发表于 2018-12-12 11:28
tmp = np.zeros((len(times_list),53,81), dtype=np.float64)
hcc = nc_obj.variables["hcc"]   # 不清楚 ...

太厉害了!!!用for循环得20分钟,改了一下只要3秒,为什么这么快呢?太厉害了
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

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

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

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