第三版WRF及WRF-Chem预处理以及后处理的Python脚本.zip

import hashlib import multiprocessing import os import pickle import sys from datetime import datetime, timedelta from multiprocessing import Pool import numpy as np from loguru import logger import xarray as xr import xesmf as xe from area import area # 获取脚本所在目录和脚本名称 my_dirname, my_filename = os.path.split(os.path.abspath(sys.argv[0])) # 其他指定路径 CB05_DIR = "/home/fishercat/Build_WRF/PreChem/MEIC-2016-CB05-Section" #meic cb05文件存储路径。内部为YYYYMM的文件夹，文件夹内存放每个月都meic排放源文件 wrfinput_file = "/home/fishercat/Build_WRF/Examples/test/wrfinput_d01" #wrfinput文件位置 wrfchemi_save_dir = "/home/fishercat/Build_WRF/PreChem/MEIC_temp" #wrfchemi文件存放目录 hour_step = 12 # meic经纬度,格点的位于每个网格的中心点 lon = np.arange(70.125, 150, 0.25, dtype=np.float32) lat = np.arange(10.125, 60, 0.25, dtype=np.float32) lon, lat = np.meshgrid(lon, lat) # 排放源高度分布 emission_height_distribution = {"agriculture": [1.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000], "industry": [0.602, 0.346, 0.052, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000], "power": [0.034, 0.140, 0.349, 0.227, 0.167, 0.059, 0.024, 0.000, 0.000, 0.000, 0.000], "residential": [0.900, 0.100, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000], "transportation":[0.950, 0.050, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000]} # 排放源时间分布,BJT 0-23 emission_time_distribution = {"agriculture": [1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000], "industry": [0.504, 0.504, 0.504, 0.504, 0.504, 0.648, 0.792, 0.936, 1.080, 1.632, 1.632, 1.632, 1.632, 1.632, 1.632, 1.632, 1.632, 1.080, 0.936, 0.792, 0.648, 0.504, 0.504, 0.504], "power": [0.840, 0.708, 0.576, 0.600, 0.732, 0.804, 0.852, 0.972, 1.104, 1.200, 1.212, 1.200, 1.224, 1.236, 1.236, 1.224, 1.176, 1.104, 1.116, 1.128, 1.080, 1.008, 0.948, 0.720], "residential": [0.336, 0.312, 0.408, 0.216, 0.408, 1.080, 1.176, 1.248, 1.224, 1.584, 2.736, 2.520, 0.984, 0.696, 0.552, 0.984, 1.680, 2.328, 1.488, 0.840, 0.264, 0.360, 0.336, 0.240], "transportation":[0.480, 0.348, 0.252, 0.216, 0.192, 0.312, 0.720, 1.416, 1.536, 1.440, 1.392, 1.272, 1.212, 1.368, 1.380, 1.416, 1.428, 1.524, 1.380, 1.128, 1.068, 1.008, 0.864, 0.648]} # # cbmz_mosaic排放变量 # cbmz_mosaic_var = ["E_SO2","E_NO","E_ALD","E_HCHO","E_ORA2","E_NH3","E_HC3","E_HC5","E_HC8","E_ETH","E_CO", # "E_OL2","E_OLT","E_OLI","E_TOL","E_XYL","E_KET","E_CSL","E_ISO","E_PM25I","E_PM25J","E_SO4I", # "E_SO4J","E_NO3I","E_NO3J","E_ORGI","E_ORGJ","E_ECI","E_ECJ","E_PM_10","E_NO"] # 无机气体的摩尔质量 inorganic_gas_mole_weight = {'CO':28, 'CO2':44, 'NH3':17, 'NOx':31.6, 'SO2':64} #NOx中由10%是NO2,90%是NO,所以平均摩尔质量是46*0.1+30*0.9=31.6 # 排放的周内变化。周日为0. week_emiss_factor = { "agriculture": [1.00,1.00,1.00,1.00,1.00,1.00,1.00], "industry": [0.80,1.08,1.08,1.08,1.08,1.08,0.80], "power": [0.85,1.06,1.06,1.06,1.06,1.06,0.85], "residential": [0.80,1.08,1.08,1.08,1.08,1.08,0.80], "transportation": [0.79,1.02,1.06,1.08,1.10,1.14,0.81] } # 计算指定中心纬度的经纬度网格的面积。 def ll_area(lat,res=0.25): ''' input: lat: 经纬度网格中心点纬度组成的数组,np.2darray res: 经纬度网格分辨率/边长,单位:度 return: 面积组成的数组,单位km2. TODO: 沿海地区的经纬度网格内陆地面积应该比网格面积小，实际上的平均排放速率应该通过陆地面积计算。 以后这里应该直接返回一个numpy数组，因为meic网格和海陆分布是固定的，不需要每次临时计算 ''' startlon=0 #任意起始经度,随便是几都可以。这里选0度。 return_area = np.zeros_like(lat) isize,jsize = return_area.shape for i in range(isize): for j in range(jsize): obj = {'type':'Polygon', 'coordinates':[ #网格的四个顶点经纬度 [[startlon,lat[i,j]-res/2.0],[startlon,lat[i,j]+res/2.0],[startlon+res,lat[i,j]+res/2.0],[startlon+res,lat[i,j]-res/2.0]] ] } return_area[i,j] = area(obj) return return_area # 插值程序,从meic网格插值到wrf网格 def meic2wrf(lon_inp,lat_inp,lon,lat,emis,interp_method = 'bilinear'): ''' input: lon_inp: wrf的经度网格(wrfinput["XLONG"]) lat_inp: wrf的纬度网格(wrfinput["XLAT"]) lon: meic的经度网格,np.2darray lat: meic的纬度网格,np.2darray emis: 转化了单位后的meic排放,np.2darray interp_method: xesmf的插值方法 return: wrf投影的排放源数据 ''' grid_out = {'lon': lon_inp,'lat': lat_inp} grid_in = {'lon': lon,'lat': lat} regridder = xe.Regridder(grid_in, grid_out, interp_method,reuse_weights=False) emis_inp = regridder(emis) return emis_inp def avg_hour(iemiss,emiss_year,emiss_month): #计算本月有多少个等效小时 #等效小时： 假设每天都排放系数都相同（都是1），日内每小时的排放都相同，那么等效小时内的排放应该是这个月的总排放除以总小时数。 #这个值直接乘以各种系数就可以作为排放 ''' iemiss: 排放源种类(5种） emiss_year,emiss_month: meic排放源头的时间 ''' avg_hour_count = 0 start_time = datetime(int(emiss_year),int(emiss_month),1,0) #本月的开始时间,bjt while start_time.strftime("%m") == emiss_month: avg_hour_count += week_emiss_factor[iemiss][int(start_time.strftime("%w"))] * emission_time_distribution[iemiss][int(start_time.strftime("%H"))] start_time += timedelta(hours=1) return avg_hour_count #结束 def convert_unit(var,value,iemiss,emiss_year="2016",emiss_month="01"): ''' var: 要变化的变量名 value: 值，二维数组 iemiss: 排放源种类(5种） emiss_year,emiss_month: meic排放源头的时间 ''' #计算本月有多少个等效小时 #等效小时: 星期变化和日变化系数都是1的小时 avg_hour_count = avg_hour(iemiss,emiss_year,emiss_month) #结束 # _,len_month = calendar.monthrange(emiss_year,emiss_month) #获取排放源对应月份的天数 if var in ['CO', 'CO2', 'NH3', 'NOx', 'SO2', ]: #inorganic gas: ton/(grid.month) to mole/(km2.h) emiss = value*1e12/(ll_area(lat, 0.25)*avg_hour_count *inorganic_gas_mole_weight[var]) elif var in ['BC', 'OC', 'PM25', 'PM10', ]: # aerosol: ton/(grid.month) to ug/(m2.s) emiss = value*1e12/(ll_area(lat, 0.25)*avg_hour_count*3600) else: # organic gas: million_mole/(grid.month) to mole/(km2.h) emiss = value*1e12/(ll_area(lat, 0.25)*avg_hour_count) return emiss #返回实际上是当月等效小时数 def pickle_read(pickle_file):# 是否存在pickle文件，如果存在则读取 flag = False try: with open(pickle_file,"rb") as f: return_dict = pickle.load(f) flag = True logger.success("success load "+pickle_file) except: return_dict = {} return flag,return_dict def md5_value(file_name): #计算wrfinput文件的md5值 ''' 每一个wrfinput文件都拥有唯一的md5值 通过md5区分不同wrfinput对应的interp_meic_emission