-
Notifications
You must be signed in to change notification settings - Fork 2
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #12 from chen-0126/master
Add support for CMAQ and bug fixes
- Loading branch information
Showing
32 changed files
with
689 additions
and
280 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Binary file not shown.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Binary file not shown.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Binary file not shown.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,8 @@ | ||
| Read before running | ||
| ---------------------- | ||
|
|
||
| 1, Set parameters(time, directory and dimension etc.) in "for_CMAQ.py". | ||
|
|
||
|
|
||
|
|
||
|
|
Empty file.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,137 @@ | ||
| """ | ||
| # Author: Wencong Chen | ||
| # Date: 2023-06-21 | ||
| # Purpose: Convert model-ready emission file to meet CMAQ's input requirements. | ||
| (height allocation, merge files of different sectors, convert pollutants and projection) | ||
| """ | ||
|
|
||
| import time | ||
| start_time = time.clock() | ||
| #Set current working directory | ||
| from inspect import getsourcefile | ||
| dir_run = os.path.dirname(os.path.abspath(getsourcefile(lambda:0))) | ||
| if not dir_run in sys.path: | ||
| sys.path.append(dir_run) | ||
|
|
||
| import os | ||
| from emips.spatial_alloc import GridDesc | ||
| from emips import ge_data_dir | ||
| from emips.utils import SectorEnum | ||
| from emips.chem_spec import get_model_species_CMAQ | ||
| from collections import OrderedDict | ||
| import datetime | ||
| import height | ||
| import merge | ||
| import transform | ||
| import proj | ||
|
|
||
| #Choose Chemical mechanism | ||
| mechanism_name = 'cb05' | ||
| print('--------CHEMICAL MECHANISM: {}--------'.format(mechanism_name.upper())) | ||
|
|
||
| #Set year and month | ||
| year = 2019 | ||
| month = 1 | ||
|
|
||
| #Set input&output file path | ||
| dir_inter = os.path.join(r'G:\CMAQ', mechanism_name, r'merge\{0:}\{0:}{1:>02d}'.format(year, month)) | ||
|
|
||
| #Set original grid | ||
| original_proj = geolib.projinfo() | ||
| original_grid = GridDesc(original_proj, x_orig=64., x_cell=0.25, x_num=324, | ||
| y_orig=15., y_cell=0.25, y_num=180) | ||
|
|
||
| #Set target grid | ||
| ncols = 164 #the number of x dimension | ||
| nrows = 97 #the number of y dimension | ||
| nlays = 7 #the number of z dimension | ||
| xcell = 36000.0 #Grid spacing in x dimensions | ||
| ycell = 36000.0 #Grid spacing in y dimensions | ||
| proj_name = 'lcc' #Name of the projection ('lcc' is Lambert Conformal Conic projection) | ||
| xcent = 110.0 #Longitude of projection center | ||
| ycent = 34.0 #Latitude of projection center | ||
| p_alp = 25.0 #First standard parallel | ||
| p_bet = 40.0 #Second standard parallel | ||
| p_gam = 110.0 | ||
| radius = 6370000.0 #Radius of the sphere(unit: meter) | ||
|
|
||
| if ncols % 2 == 0: | ||
| xorig = (ncols-2) / 2.0 * xcell + xcell / 2.0 | ||
| else: | ||
| xorig = (ncols-1) / 2.0 * xcell | ||
| if nrows % 2 == 0: | ||
| yorig = (nrows-2) / 2.0 * ycell + ycell / 2.0 | ||
| else: | ||
| yorig = (nrows-1) / 2.0 * ycell | ||
| target_proj = geolib.projinfo(proj=proj_name, lon_0=xcent, lat_0=ycent, lat_1=p_alp, lat_2=p_bet, a=radius, b=radius) | ||
| target_grid = GridDesc(target_proj, x_orig=-xorig, x_cell=xcell, x_num=ncols, | ||
| y_orig=-yorig, y_cell=ycell, y_num=nrows) | ||
|
|
||
| #Set path of the vertical allocate file | ||
| vert_prof_file = os.path.join(ge_data_dir, 'height.txt') | ||
|
|
||
| #Set sectors that need to be processed | ||
| sectors = [SectorEnum.INDUSTRY, SectorEnum.AGRICULTURE, SectorEnum.ENERGY, | ||
| SectorEnum.RESIDENTIAL, SectorEnum.TRANSPORT] | ||
|
|
||
| #Acquisition of model species (gases and aerosol) and aerosol names based on chemical mechanisms | ||
| out_species, out_species_aer = get_model_species_CMAQ(mechanism_name) | ||
| varlist = '' | ||
| for i in out_species: | ||
| varlist = varlist + '{:<16s}'.format(i) | ||
|
|
||
| #Set global attributes | ||
| gattrs = OrderedDict() | ||
| #gattrs['Conventions'] = 'CF-1.6' | ||
| gattrs['TITLE'] = "Created using MeteoInfo, mechanism: {}".format(mechanism_name.upper()) | ||
| gattrs['IOAPI_VERSION'] = "{:<80s}".format('$Id: @(#) ioapi library version 3.1 $') | ||
| gattrs['EXEC_ID'] = "{:<80s}".format('????????????????') | ||
| gattrs['FTYPE'] = 1 | ||
| gattrs['CDATE'] = datetime.date(year=year, month=month, day=1).strftime("%Y%j") #test | ||
| gattrs['CTIME'] = 81255 #test | ||
| gattrs['WDATE'] = datetime.date(year=year, month=month, day=1).strftime("%Y%j") #test | ||
| gattrs['WTIME'] = 81255 #test | ||
| gattrs['SDATE'] = datetime.date(year=year, month=month, day=1).strftime("%Y%j") | ||
| gattrs['STIME'] = 0 | ||
| gattrs['TSTEP'] = 10000 | ||
| gattrs['NTHIK'] = 1 | ||
| gattrs['NCOLS'] = ncols | ||
| gattrs['NROWS'] = nrows | ||
| gattrs['NLAYS'] = nlays | ||
| gattrs['NVARS'] = len(out_species) | ||
| gattrs['GDTYP'] = 2 | ||
| gattrs['P_ALP'] = p_alp | ||
| gattrs['P_BET'] = p_bet | ||
| gattrs['P_GAM'] = p_gam | ||
| gattrs['XCENT'] = xcent | ||
| gattrs['YCENT'] = ycent | ||
| gattrs['XORIG'] = -xorig | ||
| gattrs['YORIG'] = -yorig | ||
| gattrs['XCELL'] = xcell | ||
| gattrs['YCELL'] = ycell | ||
| gattrs['VGTYP'] = 2 #test | ||
| gattrs['VGTOP'] = 10000.0 | ||
| gattrs['VGLVLS'] = 1.0, 0.995, 0.988, 0.98, 0.97, 0.956, 0.938, 0.893, 0.839, 0.777, 0.702, 0.582, 0.4, 0.2, 0.0 | ||
| gattrs['GDNAM'] = "{:<16s}".format('M_36_08CHINA') | ||
| gattrs['UPNAM'] = "" | ||
| gattrs['VAR-LIST'] = varlist | ||
| gattrs['FILEDESC'] = "" | ||
| gattrs['HISTORY'] = "" | ||
|
|
||
| #Run all scripts | ||
| print('Allocate according to height...') | ||
| height.run(year, month, dir_inter, original_grid, sectors, nlays, vert_prof_file) | ||
| print('Merge data from different sectors...') | ||
| merge.run(year, month, dir_inter, original_grid, sectors, nlays) | ||
| print('Distribution of particulate matter and change unit...') | ||
| transform.run(year, month, dir_inter, original_grid, out_species, out_species_aer, nlays) | ||
| print('Conversion projection...') | ||
| proj.run(year, month, dir_inter, original_grid, target_grid, out_species, out_species_aer, gattrs, mechanism_name, nlays) | ||
| print('-----------------------------') | ||
| print('--------All finished!--------') | ||
| print('-----------------------------') | ||
|
|
||
| #Calculate running time | ||
| end_time = time.clock() | ||
| time = (end_time - start_time) / 60.0 | ||
| print('Time: {:.2f}min'.format(time)) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,93 @@ | ||
| from emips.spatial_alloc import GridDesc | ||
| from collections import OrderedDict | ||
| from mipylib import dataset | ||
| import mipylib.numeric as np | ||
| import os | ||
| from emips.utils import emis_util | ||
| import emips | ||
|
|
||
| out_species_unit = ['PEC', 'POA', 'PMFINE', 'PNO3', 'PSO4', 'PMC'] | ||
|
|
||
|
|
||
| def run(year, month, dir_inter, model_grid, sectors, nlays, z_file): | ||
| """ | ||
| Allocate data to different heights. | ||
| :param year: (*int*) Year. | ||
| :param month: (*int*) Month. | ||
| :param dir_inter: (*string*) The directory where data is stored. | ||
| :param model_grid: (*GridDesc*) Model data grid describe. | ||
| :param sectors: (*GridDesc*) The sectors need to be processed. | ||
| :param nlays: (*int*) The z-dimension of the output data. | ||
| :param z_file: (*string*) The path of the vertical allocate file. | ||
| """ | ||
| print('Define dimension and global attributes...') | ||
| tdim = np.dimension(np.arange(24), 'HOUR') | ||
| ydim = np.dimension(model_grid.y_coord, 'ROW', 'Y') | ||
| xdim = np.dimension(model_grid.x_coord, 'COL', 'X') | ||
| zdim = np.dimension(np.arange(nlays), 'LAY') | ||
| dims = [tdim, zdim, ydim, xdim] | ||
|
|
||
| gattrs = OrderedDict() | ||
| gattrs['Conventions'] = 'CF-1.6' | ||
| gattrs['Tools'] = 'Created using MeteoInfo' | ||
|
|
||
| for sector in sectors: | ||
| fn = dir_inter + '\emis_{}_{}_{}_hour.nc'.format(sector.name, year, month) | ||
| print('File input: {}'.format(fn)) | ||
| dimvars = [] | ||
| if os.path.exists(fn): | ||
| f = dataset.addfile(fn) | ||
| for var in f.varnames: | ||
| if var == 'lat' or var == 'lon': | ||
| continue | ||
| else: | ||
| dimvar = dataset.DimVariable() | ||
| dimvar.name = var | ||
| dimvar.dtype = np.dtype.float | ||
| dimvar.dims = dims | ||
| dimvar.addattr('description', "EMISSION_{}".format(var)) | ||
| if var in out_species_unit: | ||
| dimvar.addattr('units', 'g/m2/s') | ||
| else: | ||
| dimvar.addattr('units', 'mole/m2/s') | ||
| dimvars.append(dimvar) | ||
|
|
||
| out_fn = dir_inter + '\emis_{}_{}_{}_hour_height.nc'.format(sector.name, year, month) | ||
| print('Create output file:{}'.format(out_fn)) | ||
| ncfile = dataset.addfile(out_fn, 'c', largefile=True) | ||
| ncfile.nc_define(dims, gattrs, dimvars) | ||
|
|
||
| data = np.zeros((tdim.length, zdim.length, ydim.length, xdim.length)) | ||
| dd = np.zeros((tdim.length, zdim.length, ydim.length, xdim.length)) | ||
| # get vertical profiles | ||
| scc = emis_util.get_scc(sector) | ||
| vertical_pro = emips.vertical_alloc.read_file(z_file, scc) | ||
| if nlays != len(vertical_pro.get_ratios()): | ||
| print('----------Warning-----------\nThe z-dimensions in the vertical profile file and the output file do not match.\n----------Warning-----------') | ||
| # read, merge and output | ||
| if round(vertical_pro.get_ratios()[0], 2) != 1.0: | ||
| print('Allocating: {}'.format(sector.name)) | ||
| else: | ||
| print('Do not need to be allocated: {}'.format(sector.name)) | ||
| print('Write data to file...') | ||
| for var in f.varnames: | ||
| if var == 'lat' or var == 'lon': | ||
| continue | ||
| else: | ||
| print(var) | ||
| dd[:, 0, :, :] = f[var][:] | ||
| if round(vertical_pro.get_ratios()[0], 2) == 1.0: | ||
| data[:, 0, :, :] = dd[:, 0, :, :] | ||
| else: | ||
| for lay in np.arange(len(vertical_pro.get_ratios())): | ||
| data[:, lay, :, :] = dd[:, 0, :, :] * vertical_pro.get_ratios()[lay] | ||
| # Turn nan to zero | ||
| data[data == np.nan] = 0 | ||
| ncfile.write(var, data) | ||
| ncfile.close() | ||
| f.close() | ||
| else: | ||
| print('File not exist: {}'.format(fn)) | ||
| continue | ||
| print('Allocate of height finished!') |
Oops, something went wrong.