Python3.CMIP6数据边下载边裁剪

墨家大宝

师姐想下载亿些CMIP6数据，全球日数据太大了硬盘不够，而且只想要中国范围的月平均数据。她推给我一篇介绍在shell脚本中用wget+cdo实现类似需求的方式，而我并不懂shell和cdo，大概看了一下后觉得用Python应该可以做得更好一些。

在我理解中，shell脚本中用wget+cdo的实现方法为：循环使用wget下载原始数据到本地硬盘上，每下载完一条后通过cdo命令将数据处理成所需要的另存到硬盘，再将原始数据删除。这样做思路清晰，但有两个地方我认为可以改进：

• 多线程下载
• 原始数据直接在内存中就处理，省去硬盘写入、读取、删除的开销，硬盘上只保存最终所需要的数据
  

思路为： urllib 多线程下载，下载数据用 io.BytesIO 保存在内存， xarray 读取下载数据并处理保存到硬盘。

1. from urllib.request import Request, urlopen
   
2. from io import BytesIO
   
3. from multiprocessing.dummy import Pool
   
4. from os.path import exists
   
5. import json
   
6. from datetime import datetime
   
7. 
   
8. 
   
9. import xarray
   
10. 
    
11. 
    
12. from get_url import nc_urls
    
13. 
    
14. 
    
15. threads = 10
    
16. lon_w, lon_e = 70, 140
    
17. lat_s, lat_n = 15, 55
    
18. save_to = r'./cmip6_china_monthmean'
    
19. 
    
20. 
    
21. base_headers = {
    
22.     'user-agent': 'Mozilla/5.0 (Windows NT 10.0; Win64; x64)  '
    
23.                   'AppleWebKit/537.36 (KHTML, like Gecko)  '
    
24.                   'Chrome/74.0.3729.131 Safari/537.36',
    
25.     'Connection': 'Keep-Alive'
    
26. }
    
27. pool = Pool(threads)
    
28. lon_slice = slice(lon_w, lon_e)
    
29. lat_slice = slice(lat_s, lat_n)
    
30. def worker(info):
    
31.     nc_url, pos_start, pos_end, n, cache = info
    
32.     headers = {**base_headers, 'range': f'bytes={pos_start}-{pos_end}'}
    
33.     while True:
    
34.         try:
    
35.             cache[n] = urlopen(Request(nc_url, headers=headers), timeout=5).read()
    
36.             return
    
37.         except Exception as e:
    
38.             continue
    
39. 
    
40. done_json = f'{save_to}/done.json'
    
41. if exists(done_json):
    
42.     with open(done_json) as f:
    
43.         done = json.load(f)
    
44. else:
    
45.     done = {}
    
46. 
    
47. for nc_name, nc_url in nc_urls.items():
    
48.     if nc_name in done:
    
49.         print(f'{nc_name}\n已处理过，跳过')
    
50.         continue
    
51.     cache = [None] * threads
    
52.     print(f'{nc_name}\n{datetime.now()}开始下载')
    
53.     while True:
    
54.         try:
    
55.             nc_bytes = int(urlopen(Request(nc_url, headers=base_headers), timeout=5).headers['content-length'])
    
56.             break
    
57.         except Exception as e:
    
58.             continue
    
59.     subsize = nc_bytes // threads
    
60.     worker_info = [[nc_url, x*subsize, (x+1)*subsize-1, x, cache] for x in range(threads)]
    
61.     worker_info[-1][2] = nc_bytes - 1
    
62.     pool.map(worker, worker_info)
    
63.     print(f'{datetime.now()}下载完成')
    
64.     f = BytesIO(b''.join(cache))
    
65.     del cache
    
66.     f.seek(0)
    
67.     ds = xarray.open_dataset(f).sel(lon=lon_slice, lat=lat_slice)
    
68.     ds = ds.groupby(ds.time.astype('datetime64[M]')).mean()
    
69.     ds.to_netcdf(f'{save_to}/{nc_name[:-3]}_china_monthmean.nc')
    
70.     print(f'{datetime.now()}处理完成')
    
71.     del ds, f
    
72.     done[nc_name] = nc_url
    
73.     with open(done_json, 'w') as f:
    
74.         json.dump(done, f, indent=4)
    
75. 
    
76. pool.close()
    
77. pool.join()

复制代码

其中 get_url.py 是用来从网站提供的sh脚本中提取下载链接的，内容如下：

1. from os.path import basename, exists
   
2. from glob import glob
   
3. import re
   
4. import json
   
5. 
   
6. wget_sh_glob = './pr_ssp245/wget*.sh'
   
7. nc_urls_json = './pr_ssp245/nc_urls.json'
   
8. 
   
9. if exists(nc_urls_json):
   
10.     with open(nc_urls_json) as f:
    
11.         nc_urls = json.load(f)
    
12. else:
    
13.     nc_urls = {}
    
14.     for sh_path in glob(wget_sh_glob):
    
15.         with open(sh_path) as f:
    
16.             txt = f.read()
    
17.         for nc_url in re.findall(r'(http://\S+?.nc)\'', txt):
    
18.             nc_name = basename(nc_url)
    
19.             if nc_name in nc_urls:
    
20.                 print(f'{nc_name}')
    
21.             else:
    
22.                 nc_urls[nc_name] = nc_url
    
23.     with open(nc_urls_json, 'w') as f:
    
24.         json.dump(nc_urls, f, indent=4)

复制代码

效果如图：

收获如下：

• 虽然不太优雅，但第一次成功尝试多线程，以后有的是进步空间
• 第一次尝试直接读取内存中的二进制nc文件（需要安装 h5netcdf 包）
  

存在的问题：

• 有些数据从如1850年开始计时，所以需要安装`cftime`包
• 下载过程中没有失败跳过判断（需求多的话我再加上吧）
• 第一次用 io.BytesIO ，尝试过多线程直接按照各自位置向 io.BytesIO 中存，但失败了，所以妥协了先用字节串列表存数据，再拼接到 io.BytesIO 中，希望擅长多线程或者相关“内存临时文件”编程的前辈指教一下
  

不想去气象局

CMIP6本身就有月资料的呀。

wjy_ecnu

CMIP6的wget脚本好像是一次超过1000个文件就只能下载前1000个 这个会不会影响

墨家大宝

wjy_ecnu 发表于 2020-8-29 19:33
CMIP6的wget脚本好像是一次超过1000个文件就只能下载前1000个 这个会不会影响

我不清楚wget脚本干了啥，要是1000个之后的下载链接在wget脚本中也有的话应该都可以用正则表达式提取出来然后下载处理

墨家大宝

不想去气象局 发表于 2020-8-29 19:27
CMIP6本身就有月资料的呀。

我其实不了解CMIP6，就是师姐有这么个需求。是不是有个降尺度啥的没有月资料？

随缘

师姐必须以身相许

吕青

随缘 发表于 2020-8-29 21:26
师姐必须以身相许

哈哈哈哈哈哈哈哈

井中月

墨家大宝 发表于 2020-8-29 19:43
我不清楚wget脚本干了啥，要是1000个之后的下载链接在wget脚本中也有的话应该都可以用正则表达式提取出来 ...

CMIP资料的那个链接离开那个wget脚本下不了，那个脚本里有类似权限认证的东西

墨家大宝

井中月 发表于 2020-8-31 09:58
CMIP资料的那个链接离开那个wget脚本下不了，那个脚本里有类似权限认证的东西

那可能数据来源的网站不一样吧。我猜你的wget脚本权限认证可能也就是在request的header甚至是在url后面加id、密码之类的参数吧

Lancelot

1. 这么大的内存怎么的硬盘空间那么小
2. CMIP6 wget的脚本经过处理可以绕过认证登录
3. 你师姐下的数据有逐月数据