需要对指定格点排放源的数值进行更改,写了ncl和python两个脚本
图为关闭湖南湖北排放源前后对比
code:
需要经纬度格点数的txt文件del_d01.txt,del_d02.txt
NCL:
row = numAsciiRow("./del_d01.txt")
col = numAsciiCol("./del_d01.txt")
del_d01 = asciiread("./del_d01.txt",(/row,col/),"float")
printVarSummary(del_d01)
year=2019
a = addfile("./nothb_"+year+"/wrfchemi_12z_d01", "w")
my_e_iso = a->E_ISO
do idel =0,row-1 ;纬度的格点数范围
ilon = toint(del_d01(idel,0)-1)
ilat = toint(del_d01(idel,1)-1)
my_e_iso(:,:,ilat,ilon)=0
end do
a->E_ISO = my_e_iso
Python:
import os
import numpy as np
import Ngl,Nio
# emission year
yyyy_emi = 2015
#-- read the ascii data
f_grid = open("./del_d01.txt",'r')
data = f_grid.readlines() #-- data: type list
nrows = len(data)
#-- assign lists to append elements
lat0 = []
lon0 = []
for i in data:
line = i.strip()
# print(line)
cols = line.split()
lon0.append(cols[0])
lat0.append(cols[1])
#-- convert string to float
lat = np.array(lat0).astype(int)
lon = np.array(lon0).astype(int)
######################
#-- nc data file name
fname = "./nothb_"+str(yyyy_emi)+"/wrfchemi_00z_d01"
#-- open file
f = Nio.open_file(fname, "w")
#-- read emission arrays
my_e_iso = f.variables["E_ISO"]
#-- convert data change thb to 0
for i in list(range(nrows)):
ilat=lat[i]
ilon=lon[i]
my_e_iso[:,:,ilat,ilon] =0
#-- assign values --> write data to file
f.variables["E_ISO"].assign_value(my_e_iso)
#-- close output stream (not necessary)
f.close()