Hi, everyone. I am working on a research focusing on the pollutants from fire. However, I want to edit the emission file to simulate fire in different place and compare the results. Do you know what tool I can use to modify emission file. Especially the value in specific location.
There are lots of tools and ad hoc methods – and everyone will have their own preference on how to approach this. I’m going to detail some very simplistic ad hoc methods.
The approach will differ if you intend to add emissions in a 2D or inline emissions file. Below are super simple approaches that assume you know “where” the emissions go. At the end, I added a description of how to know the “where” in grid indexes and projected X/Y.
These are just suggestions – no warranty expressed or implied. Test the results and visualize them to make sure you are getting what you intended.
2D Emissions File
Below is a super simple python script that assumes you want to add 10 moles/s of NOx at a specific location to an existing 2D file. Then, run the code below with your own modifications.
#!/usr/bin/env python
# requires you have netcdf4-python installed (e.g., pip install pyproj)
import netCDF4
epath = '/path/to/existing/gridded_2d_emissions_file.nc'
f = netCDF4.Dataset(epath, mode='r+s')
t=17 # 17Z
k=0 # LAY=1 because indices are 0-based in C-type languages
j=100 # ROW=101
i=100 # COL=101
f.variables['NO'][t, k, j, i] += 9
f.variables['NO2'][t, k, j, i] += 1
# ... other species and values...
f.close()
This 2D approach would also work for a 3D gridded file with obvious modifications. To get the existing file, you can copy one you already have or generate an empty file or using m3fake from the IOAPI package.
Inline Emissions File
Another approach would be to create an inline point source file for just one point. To do that, I’d recommend a three step process.
- Use
ncksfrom the nco package to generate Common Data Language templates from an existing set of ptfire files (stack group and inlin) - Manually edit those template files to match your intention.
- Use ncgen from the netcdf package to transform the templates into real files
Below is a bash script to generate the templates and then another to
# Requires NCO package and ncgen.
# Get the first fire location from a stack group file as a template
ncks -a -dROW,0 /path/to/existing/ptfire_stack_group_file.nc -o stack_group_one_fire.nc
# Get the first (same) fire location from an inlin emission file as a template
ncks -a -dROW,0 /path/toexisting/ptfire_inline_file.nc -o inline_one_fire.nc
ncdump stack_group_one_fire.nc > stack_group_one_fire.cdl
ncdump inline_one_fire.nc > inline_one_fire.cdl
Now manually edit the cdl files to match what you want.
- In both files:
a. edit NROWS=1 instead of the original number of columns,
b. edit SDATE/STIME and TFLAG values. - In stack group file, make sure
a. the XLOCA/YLOCA variables are in meters from origin
b. the ROW/COL values are j+1/i+1 because Fortran uses 1-based indexing
b. the stack parameters are appropriate for the fire you want - In the inline emissions file, make sure emissions you add are in units consistent with the variable
Then run this script to generate the files.
ncgen -o stack_group_one_fire.nc stack_group_one_fire.cdl
ncgen -o inline_one_fire.nc inline_one_fire.cdl
Because the ncks/ncgen approach is so manual, I likely haven’t been able to perfectly describe the manual editing process. So, it may take you a bit of iteration to get it fully working.
Coordinates
Both the examples assume that you know the j/i and or the x/y. The script below would help you know which j and i to use for a gridded file or the x/y used in the inline file. It requires values from the GRIDDESC file. I’ve included examples that work with a Lambert Conic Conformal projection.
#!/usr/bin/env python
# requires you have pyproj installed (e.g., pip install pyproj)
import pyproj
# Example GRIDDESC values for for EPA's 12US2 domain
# you'll need to edit this to match your GRIDDESC values.
GDTYP = 2
YCENT = 40.
XCENT = -97.
P_ALP = 33.
P_BET = 45.
XORIG = -2412000.
YORIG = -1620000.
XCELL = 12000.
YCELL = 12000.
# Conversion to proj4 using parameters from GRIDDESC
proj4str = f'+proj=lcc +lat_0={YCENT} +lon_0={XCENT} +lat_1={P_ALP} +lat_2={P_BET} +x_0=0 +y_0=0 +R=6370000 +units=m +no_defs'
proj = pyproj.Proj(proj4str)
lon = -98
lat = 40
x, y = proj(lon, lat)
i = int((x - XORIG) / XCELL)
j = int((y - YORIG) / YCELL)
print(f'lon={lon:.6f}; lat={lat:.6f}')
print(f'x={x:.2f}; y={y:.2f}')
print(f'i={i:d}; j={j:d}')
print(f'COL={i + 1:d}; ROW={j + 1:d}')
3 Likes
Barron laid out some excellent simple options to generate and modify emissions files.
If you are starting from an emissions modeling platform it is also an option to make inventory-level modifications, such as changes to locations or emissions values, and rerun SMOKE to generate new model-ready emissions inputs.
Complete packages of SMOKE inputs, inventories, run scripts, and executables for EQUATES are available in the SMOKE_INPUTS folder of the EQUATES Google Drive linked here:
Packages for the EMPs, including the 2020 - 2022 platforms, are available in the data files sections on the Emissions Modeling Platform website:
2 Likes
Thank you so much for your reply. The content is very helpful!
HI, Barron. I have tried the second method. When I was trying to edit the CDL files with vscode or vim. The file has some strange characters. Do you know the correct way to edit the CDL file?
Thank you, James. I have never used SMOKE before so I will try Barron’s method first. And it seems that I still need to consider learning SMOKE
Using ncks directly might have given you different results, so I edited my response above. It now uses ncks and then ncdump. The contents of CDL is common data language and it should look like:
netcdf stack_group_one_fire {
dimensions:
TSTEP = UNLIMITED ; // (1 currently)
LAY = 1 ;
ROW = 1 ;
COL = 1 ;
VAR = 17 ;
DATE-TIME = 2 ;
variables:
float ACRESBURNED(TSTEP, LAY, ROW, COL) ;
ACRESBURNED:long_name = "ACRESBURNED " ;
ACRESBURNED:units = "acres/day " ;
ACRESBURNED:var_desc = "number of acres burned for a fire in one day " ;
int COL(TSTEP, LAY, ROW, COL) ;
COL:long_name = "COL " ;
COL:units = "none " ;
COL:var_desc = "Grid column number " ;
int IFIP(TSTEP, LAY, ROW, COL) ;
IFIP:long_name = "IFIP " ;
IFIP:units = "none " ;
IFIP:var_desc = "FIPS CODE " ;
int ISTACK(TSTEP, LAY, ROW, COL) ;
ISTACK:long_name = "ISTACK " ;
ISTACK:units = "none " ;
ISTACK:var_desc = "Stack group number
...
// global attributes:
:IOAPI_VERSION = "ioapi-3.2: $Id: init3.F90 185 2020-08-28 16:49:45Z coats $ " ;
:EXEC_ID = "???????????????? " ;
...
data:
ACRESBURNED =
36.36 ;
COL =
99 ;
IFIP =
4005 ;
ISTACK =
1 ;
LATITUDE =
36.566 ;
LMAJOR =
1 ;
...
}
I have put ... in places to represent more data that I am not showing. Note that there are four sections: dimensions, variables, global properties, and data. The contents of each section will make sense after you stare for a while. You’re going to be edit the global properties (NROWS=1, SDATE, STIME) and data (TFLAG, stack properties, emission values).
I can’t tell you what the data should be… You can review existing fires in your files to get a sense for the data. I _can_show you what variables you will encounter. Below are variables and their var_desc from the stack group file. The var_desc attribute gives you a sense of the contents. You’ll also see units for each variable (not shown here).
The stack_group file should have the following variables:
TFLAG:var_desc = "Timestep-valid flags: (1) YYYYDDD or (2) HHMMSS " ;
ISTACK:var_desc = "Stack group number " ;
LATITUDE:var_desc = "Latitude " ;
LONGITUDE:var_desc = "Longitude " ;
STKDM:var_desc = "Inside stack diameter " ;
STKHT:var_desc = "Stack height above ground surface " ;
STKTK:var_desc = "Stack exit temperature " ;
STKVE:var_desc = "Stack exit velocity " ;
STKFLW:var_desc = "Stack exit flow rate " ;
STKCNT:var_desc = "Number of stacks in group " ;
ROW:var_desc = "Grid row number " ;
COL:var_desc = "Grid column number " ;
XLOCA:var_desc = "Projection x coordinate " ;
YLOCA:var_desc = "Projection y coordinate " ;
IFIP:var_desc = "FIPS CODE " ;
LMAJOR:var_desc = "1= MAJOR SOURCE in domain, 0=otherwise " ;
LPING:var_desc = "1=PING SOURCE in domain, 0=otherwise " ;
ACRESBURNED:var_desc = "number of acres burned for a fire in one day " ;
The emission variables are more straight forward. But again, if you’re making an experiment, you have decide what the right values are. The inline file should have emission variables with varying units by type and a heat-flux variable (HFLUX).
TFLAG:units = "<YYYYDDD,HHMMSS>" ;
FORM:units = "moles/s " ;
FORM_PRIMARY:units = "moles/s " ;
NMOG:units = "g/s " ;
MEOH:units = "moles/s " ;
BENZ:units = "moles/s " ;
ALD2:units = "moles/s " ;
ALD2_PRIMARY:units = "moles/s " ;
NAPH:units = "moles/s " ;
ETHYLBENZ:units = "moles/s " ;
BUTADIENE13:units = "moles/s " ;
ACROLEIN:units = "moles/s " ;
CO:units = "moles/s " ;
HFLUX:units = "BTU/hr " ;
NH3:units = "moles/s " ;
NH3_FERT:units = "moles/s " ;
AACD:units = "moles/s " ;
ACET:units = "moles/s " ;
ALDX:units = "moles/s " ;
APIN:units = "moles/s " ;
CH4:units = "moles/s " ;
ETH:units = "moles/s " ;
ETHA:units = "moles/s " ;
ETHY:units = "moles/s " ;
ETOH:units = "moles/s " ;
FACD:units = "moles/s " ;
GLY:units = "moles/s " ;
GLYD:units = "moles/s " ;
IOLE:units = "moles/s " ;
ISOP:units = "moles/s " ;
ISPD:units = "moles/s " ;
IVOC:units = "moles/s " ;
KET:units = "moles/s " ;
MGLY:units = "moles/s " ;
NVOL:units = "moles/s " ;
OLE:units = "moles/s " ;
PACD:units = "moles/s " ;
PAR:units = "moles/s " ;
PRPA:units = "moles/s " ;
SOAALK:units = "moles/s " ;
TERP:units = "moles/s " ;
TOL:units = "moles/s " ;
UNR:units = "moles/s " ;
XYLMN:units = "moles/s " ;
HONO:units = "moles/s " ;
NO:units = "moles/s " ;
NO2:units = "moles/s " ;
PAL:units = "g/s " ;
PCA:units = "g/s " ;
PCL:units = "g/s " ;
PEC:units = "g/s " ;
PFE:units = "g/s " ;
PH2O:units = "g/s " ;
PK:units = "g/s " ;
PMG:units = "g/s " ;
PMN:units = "g/s " ;
PMOTHR:units = "g/s " ;
PNA:units = "g/s " ;
PNCOM:units = "g/s " ;
PNH4:units = "g/s " ;
PNO3:units = "g/s " ;
POC:units = "g/s " ;
PSI:units = "g/s " ;
PSO4:units = "g/s " ;
PTI:units = "g/s " ;
PMC:units = "g/s " ;
SO2:units = "moles/s " ;
SULF:units = "moles/s " ;
VOC_INV:units = "g/s " ;
There will be a lot you need to work through.
Hi, Barron. I can see the value name but the value is all ^@ kind of character in my editor. That is my problem.
Yours,
Cheb
Did you recreate the files using the updated recommendation:
# Requires NCO package and ncgen.
# Get the first fire location from a stack group file as a template
ncks -a -dROW,0 /path/to/existing/ptfire_stack_group_file.nc -o stack_group_one_fire.nc
# Get the first (same) fire location from an inlin emission file as a template
ncks -a -dROW,0 /path/toexisting/ptfire_inline_file.nc -o inline_one_fire.nc
ncdump stack_group_one_fire.nc > stack_group_one_fire.cdl
ncdump inline_one_fire.nc > inline_one_fire.cdl
Oh! I am sorry that I missed one step. That is my bad. I get a txt file that is actually can be edit now.
Thank you!