COSMIC掩星数据（2）-for循环中判断in nc_obj.variables.keys()的优化

前言

        最近研究需要全球的水汽剖面的观测数据，读了几篇文献，发现COSMIC掩星数据最为常用，而最新的wetPf2_repro2021_2010_002.tar.gz数据中直接就有0-60km的比湿，其中0-20km垂直分辨率50m，20-60km垂直分辨率100m，基本符合我的需求。因此急需解决的问题是：如何把掩星数据转化为全球格点数据。

        我计划分两步：第一步，读取每天的2000多个掩星数据文件，提取其三维和比湿（即经度'lon',纬度'lat',距离海平面的几何高度'MSL_alt',比湿'sph'），形成每天一个DataFrame表格，将其保存为pkl文件；第二步，读取每个月的30个左右的pkl文件，根据其三维，将其放到5°×5°的网格的对应高度上，然后做平均，得到月平均的全球比湿廓线的网格数据。

        现在的问题是每天的掩星数据量很大，大约200w行，我运行生成1天的pkl文件要40min，这是无法接受的，因此必须进行优化。

一、原始程序及运行时间

import numpy as np
import pandas as pd
from netCDF4 import Dataset
from satellite import function
import time
import os
time_start = time.time()  # 记录开始时间
path = 'D:/qqfile/wetPf2_repro2021_2010_001/'
txts = function.Read_filename(path)

# 先构造一个键分别为'lon','lat','MSL_alt','sph'，值为空list的字典
dic = {}
name = ['lon', 'lat', 'MSL_alt', 'sph']
for i in name:
	dic[i] = []

for file in txts:
	nc_obj = Dataset(file)
	if 'MSL_alt' in nc_obj.variables.keys():
		# missing_value: -999.0
		# current shape = (801,)
		# MSL_alt: 0-20km(50m) 20-60km(100m)
		MSL_alt = nc_obj['MSL_alt'][:].data  # Geometric height above MSL units: km [ 0. 60.]
		lon = nc_obj['lon'][:].data  # Geodetic longitude [-180.  360.]
		lat = nc_obj['lat'][:].data  # Geodetic latitude [-90.  90.]
		sph = nc_obj['sph'][:].data  # Specific humidity units: g/kg [ 0. 80.]
	
		# 创建三维定位和比湿的字典
		for j in range(len(MSL_alt)):
			dic['lon'].append(lon[j])
			dic['lat'].append(lat[j])
			dic['MSL_alt'].append(MSL_alt[j])
			dic['sph'].append(sph[j])
data = pd.DataFrame(dic)

data.to_pickle("D:/qqfile/tzzs_data.plk")
# read_pkl = pd.read_pickle("D:/qqfile/tzzs_data.plk")

time_end = time.time()  # 记录结束时间
time_sum = time_end - time_start  # 计算的时间差为程序的执行时间，单位为秒/s
print(time_sum)

其中 if 'MSL_alt' in nc_obj.variables.keys():是因为有些掩星文件是空的，不做判断会报错。然而正是因为这个判断，他越跑越慢，跑了40min才跑完。

二、优化程序及运行时间

将判断部分独立出来，先得到非空的文件路径，再对这些文件进行读取。

import numpy as np
import pandas as pd
from netCDF4 import Dataset
from satellite import function
import time
import os
time_start = time.time()  # 记录开始时间
path = 'D:/qqfile/wetPf2_repro2021_2010_001/'
txts = function.Read_filename(path)

# 先构造一个键分别为'lon','lat','MSL_alt','sph'，值为空list的字典
dic = {}
name = ['lon', 'lat', 'MSL_alt', 'sph']
for i in name:
	dic[i] = []

txts2 = []
for file in txts:
	nc_obj = Dataset(file)
	if 'MSL_alt' in nc_obj.variables.keys():
		txts2.append(file)

for file2 in txts2:
	nc_obj = Dataset(file2)
	# missing_value: -999.0
	# current shape = (801,)
	# MSL_alt: 0-20km(50m) 20-60km(100m)
	MSL_alt = nc_obj['MSL_alt'][:].data  # Geometric height above MSL units: km [ 0. 60.]
	lon = nc_obj['lon'][:].data  # Geodetic longitude [-180.  360.]
	lat = nc_obj['lat'][:].data  # Geodetic latitude [-90.  90.]
	sph = nc_obj['sph'][:].data  # Specific humidity units: g/kg [ 0. 80.]

	# 创建三维定位和比湿的字典
	for j in range(len(MSL_alt)):
		dic['lon'].append(lon[j])
		dic['lat'].append(lat[j])
		dic['MSL_alt'].append(MSL_alt[j])
		dic['sph'].append(sph[j])
data = pd.DataFrame(dic)

data.to_pickle("D:/qqfile/tzzs_data3.plk")
read_pkl = pd.read_pickle("D:/qqfile/tzzs_data3.plk")

time_end = time.time()  # 记录结束时间
time_sum = time_end - time_start  # 计算的时间差为程序的执行时间，单位为秒/s
print(time_sum)

最后用时7s！！！ Amazing啊

最后的最后，检查了一下两个pkl文件是否一致，发现完全一致！

优化太重要了。