Skip to content

Instantly share code, notes, and snippets.

@barronh
Last active March 17, 2023 17:48
Show Gist options
  • Save barronh/97633470997ca77b4bb949df20d9a0a3 to your computer and use it in GitHub Desktop.
Save barronh/97633470997ca77b4bb949df20d9a0a3 to your computer and use it in GitHub Desktop.
MICS to CMAQ
Display the source blob
Display the rendered blob
Raw
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
@hlbutterfly
Copy link

Following the previous post. I dig more into this. It seems like the output files for all sectors are dropping the first variable, for some reason. For instance, 'ALD2' is not in the output 'VAR-LIST' and 'NVARS' is one less, although it is showing up as a variable. So in the case of agriculture where only NH3 is valid, 'VAR-LIST' is null and 'NVARS=0'. Why this is happening?

@barronh
Copy link
Author

barronh commented Aug 30, 2021

@hlbutterfly - Try manually deleting VAR-LIST and NVARS before calling updatemeta. Does that fix it?

delattr(outf, 'VAR-LIST')
delattr(outf, 'NVARS')
outf.updatemeta()
outf.updatetflag(overwrite=True)
outf.variables.move_to_end('TFLAG', last=False)

@barronh
Copy link
Author

barronh commented Aug 30, 2021

@Gwang-Jin - make sure you're aware of this.

@hlbutterfly
Copy link

@hlbutterfly - Try manually deleting VAR-LIST and NVARS before calling updatemeta. Does that fix it?

delattr(outf, 'VAR-LIST')
delattr(outf, 'NVARS')
outf.updatemeta()
outf.updatetflag(overwrite=True)
outf.variables.move_to_end('TFLAG', last=False)

Yes! It worked! Although the two commands should not be put right before the "updatemeta" command. In that case, I got the following error message:

if update and len(varlist) != self.NVARS:
AttributeError: 'griddesc' object has no attribute 'NVARS'

I move these two commands right after outf.HISTORY = "Created from MICS.ipynb" and this worked out fine! I am getting all the required variables in the 'VAR-LIST'.

Thank you so much!

@barronh
Copy link
Author

barronh commented Aug 30, 2021

I made the changes to the gist. Can you confirm that I made them where you did?

@hlbutterfly
Copy link

I made the changes to the gist. Can you confirm that I made them where you did?

Yes. That is correct!

@hlbutterfly
Copy link

@ barronh

Just to make you aware of this.

Another potential bug that I found is that PMC and PMFINE species are double counted. In [13] is replacing PM10 with PMC and PM2.5 with PMFINE in the filename list. However, PMFINE and PMC file has already been added to the list under In [11] and ln [12] (last command). So either remove the replacing command or remove the extend command. Otherwsie, PMC and PMFINE will be doubled counted. As least, I found this mistake with my processing.

Thanks
Ling

@barronh
Copy link
Author

barronh commented Sep 3, 2021

Thank you for your updates. With your permission, I would like to add your name at the top as a contributor.

@hlbutterfly
Copy link

Thank you for your updates. With your permission, I would like to add your name at the top as a contributor.

Thanks for asking! I am fine with this but I just feel my contribution is far from being qualified as a contributor.

Another thing that I forgot to point out earlier is the calculation of PMC should be PM10 - PM2.5, instead of PM10 - OC - BC in In [12].

@barronh
Copy link
Author

barronh commented Sep 3, 2021

Good catch. Clearly a copy and paste error. I used your github handle to recognize your contribution. I put this together to be helpful, but haven't tested it like you have. I appreciate your contribution!

@atmosxiuxiu
Copy link

Hi, Dr. Henderson,
An error occurred when I executed this code( correspond to code In [31]), that is:
"CDOException: (returncode:1) cdo gridarea: Started child process "remapnn,27CN8.grid -stdatm,0 (pipe1.1)".
cdo(2) remapnn: Started child process "stdatm,0 (pipe2.1)".
cdo(3) stdatm: Set default filetype to GRIB
cdo gridarea: Open failed on >flux_regrid/area.nc<
No such file or directory"
Do you know the reason and how to fix it ?

@barronh
Copy link
Author

barronh commented Oct 25, 2021

It looks like the flux_regrid folder is not being created until after the file is made. If you make the folder, it should work. We should move the os.makedirs calls in cell 32 to cell 31.

@atmosxiuxiu
Copy link

It looks like the flux_regrid folder is not being created until after the file is made. If you make the folder, it should work. We should move the os.makedirs calls in cell 32 to cell 31.

Thank you for replying, I made the "flux_regrid" folder before executed the "cdo.gridarea( )",but there was still mistake. I was wondering if there something wrong with my CDO? And I changed the cell 5 because the original code have error when I executed, maybe related to this?
cdoo = cdo.Cdo(cdo=cdopath)==>cdoo = cdo.Cdo(cdopath)

I used the colab and executed the MICS.ipynb, the cell 31 as below:

areapath = 'flux_regrid/area.nc'
clean(areapath, overwrite=True)
**os.makedirs('flux_regrid', exist_ok=True)** # make "flux_regrid" folder

cdoo.gridarea(f'-O -remapnn,{dom}.grid -stdatm,0', output=areapath, returnCdf=True)
areaf = pnc.pncopen(areapath, format='netcdf')
projarea = areaf.variables['cell_area'][:]

The error information as below:

Could not import data from file 'flux_regrid/area.nc'
---------------------------------------------------------------------------
OSError                                   Traceback (most recent call last)
<ipython-input-46-de24e1440792> in <module>()
      3 os.makedirs('flux_regrid', exist_ok=True)
      4 
----> 5 cdoo.gridarea(f'-O -remapnn,{dom}.grid -stdatm,0', output=areapath, returnCdf=True)
      6 areaf = pnc.pncopen(areapath, format='netcdf')
      7 projarea = areaf.variables['cell_area'][:]

1 frames
/usr/lib/python3/dist-packages/cdo.py in readCdf(self, iFile)
    480     """Return a cdf handle created by the available cdf library. python-netcdf4 and scipy suported (default:scipy)"""
    481     try:
--> 482         fileObj =  self.cdf(iFile, mode='r')
    483     except:
    484       print("Could not import data from file '%s'" % iFile)

src/netCDF4/_netCDF4.pyx in netCDF4._netCDF4.Dataset.__init__()

src/netCDF4/_netCDF4.pyx in netCDF4._netCDF4._ensure_nc_success()

OSError: [Errno -101] NetCDF: HDF error: b'flux_regrid/area.nc'

@hlbutterfly
Copy link

@atmosxiuxiu
I had encountered similar errors with cdo, which took me quite long time to figure it our. Try the following command:

cdoo.remapycon(f'{dom}.grid', input=fluxpath, output=regridpath, returnCdf=False)

@Gwang-Jin
Copy link

Gwang-Jin commented Dec 14, 2021

Oh, I reproduced the same NetCDF: HDF error: b'flux_regrid/area.nc' error in[31].

When I read flux_regrid/area.nc with vim, I found

GRIB^@Ѷ^A^@^@4?b^???^A^A^@
The file flux_regrid/area.nc was not HDF or CDF but GRIB.

@Gwang-Jin
Copy link

I've solved the problem by adding a option options = '-f nc',

cdoo.gridarea( f'-O -remapnn,{dom}.grid -stdatm,0', options = '-f nc', output=areapath, returnCdf=True )

How about?

@barronh
Copy link
Author

barronh commented Dec 14, 2021

Great. I made the update.

@Gwang-Jin
Copy link

Dear Barron,

I have obtained too big value of emissions (CO : 100 times).
But I found that MICS.ipynb deal with units very seriously.
I don't know why the emissions have too big value.

Best regards,
Gwang-Jin

@Gwang-Jin
Copy link

Gwang-Jin commented Dec 30, 2021

Dear Barron,

I found why the emissions have too big value.
When converting gram [g] to mole, should the emissions be divided by some fraction?
For example,

1 moles CO to grams = 28.0101 grams.

In case of CO, MICS has unit of mass Mg/month originally

$ ncdump -h MICS_Asia_CO_2010_0.25x0.25.nc
float CO_RESIDENTIAL(time, lat, lon) ;
CO_RESIDENTIAL:long_name = "CO_RESIDENTIAL" ;
CO_RESIDENTIAL:short_name = "MICS_Asia_CO_RES" ;
CO_RESIDENTIAL:units = "Mg/month" ;

MICS.ipynb is likely to need consideration of mole-mass conversion by each chemical species.
Fortunately, only five species (CO, CO2, NOx, SO2, and NH3) need to be converted.

CO : 1 moles = 28.0101 grams
CO2 : 1 moles = 44.0095 grams
NO : 1 moles = 30.0061 grams
NO2 : 1 moles = 46.00550 grams
SO2 : 1 moles = 64.0638 grams
NH3 : 1 moles = 17.03052 grams

Best regards,
Gwang-Jin

@Gwang-Jin
Copy link

Gwang-Jin commented Dec 31, 2021

I found that MICS.ipynb has already considered mole-mass conversion regarding NH3.
How about modifying [17] :

splitfactors = {
'NOx': [e2g('NO', 0.9/30.0061), e2g('NO2', 0.1/46.00550)],
'CO': [e2g('CO', 1./28.0101)],
'SO2': [e2g('SO2', 0.97/64.0638), e2p('PSO4', 0.03)],
'NH3': [e2g('NH3', 1./17.03052)],
}

?

I am not sure that PSO4 need no mole-mass conversion. Or

'SO2': [e2g('SO2', 0.97/64.0638), e2p('PSO4', 0.03/64.0638 * 96.)],

?

After changing, I've got a reasonable CO emission [moles/s].

Thank you.

@barronh
Copy link
Author

barronh commented Dec 31, 2021

@Gwang-Jin

You're right. My code incorrectly assumed the units for gases (CO, NOx, and SO2) were in Mmol/month like the VOCs. I agree with all of your updates except for NOx. For NOx, the file emissions units are probably in "Mg/month as NO2." This is similar to using TgN instead of TgNOx. The benefit is that the mass does not depend on what the speciation split. This is very common.

So, the conversion for both NO and NO2 would use the Mw of NO2.

You're also correct that PSO4 does not need the adjustment because CMAQ uses g/s for aerosols.

I'll update the gist.

@Gwang-Jin
Copy link

I got it regarding NOx. Thank you.

@hlbutterfly
Copy link

@barronh Hi Barron! Just wondering do you have any tools or suggestions if I need to reproject gridded emission files from one Lambert Conformal Conic (LCC) projection to another LCC project with the total emissions conserved? Say I need to reproject a 27kmx27km gridded emissions into 36kmx36km grids. Cdos' remapycon seems not working and remapbil gives me close numbers (off by 0.5%). Thanks!

@barronh
Copy link
Author

barronh commented Jan 21, 2022

The best choice is using remapycon. The results still won't be perfect, but there are three basic steps to make sure it works right.

  1. Convert emissions from CMAQ-ready mass/rates and mole/rates files to flux files (moles/s to moles/m2/s).
  2. Set the flux file to the projected space using cdo's grid definition file
  3. Run remapycon

Steps 2 and 3 can be combined.

It would be better to ask this type of question thru the CMAS Center forum since it is more general than this gist.

@shumarkq
Copy link

A quick remind
If you are using different machines beyond google colab, different installation methods or versions of CDO could make "cdoo.remapycon" or "cdoo.gridarea" not working. The CDO version on Ubuntu with Colab is 1.9.3 and python3-CDO 1.3.5 here.

If you can't install old-version CDO like 1.9.3, the recent version of CDO I have tested is CDO 2.0.3 using conda.
Several changes to match original code with CDO 2.0.3

  1. install cdo bundle:
    I suggest to create a new conda env for cdo so it won't have unexpected bugs happening for conda
    conda create -n cdo
    conda install -c conda-forge cdo python-cdo xarray

  2. make sure you also install other required packages for running this gist.
    Pseudonetcdf
    pyproj

  3. Now you don't have to use python-cdo interface, you can comment out like
    #import cdo
    #cdoo = cdo.Cdo(cdopath)

  4. Change cdoo code
    change
    cdoo.remapycon(f'{dom}.grid', input=fluxpath, output=regridpath, returnCdf=False)
    to
    os.system(f'cdo remapycon,grid.{domain} {fluxpath} {regridpath}')

change
cdoo.gridarea(f'-O -remapnn,{dom}.grid -stdatm,0', options='-f nc', output=areapath, returnCdf=True)
to
os.system(f'cdo -O -f nc -gridarea -remapnn,grid.{domain} -stdatm,0 ../outputs/{domain}/flux_regrid/area.nc')

  1. Instead of creating "Projection" variable in remapped files, newer version CDO like CDO 2.0.3 here will output "crs". You also need to change code
    if k in ('time', 'lat', 'lon', 'x', 'y', 'Projection')
    to
    if k in ('time', 'lat', 'lon', 'x', 'y', 'crs')

@Gwang-Jin
Copy link

Dear @barronh ,

I'm using eta_levels in WRF :

1.000, 0.995, 0.985, 0.979, 0.973,
0.967, 0.960, 0.954, 0.949, 0.941,
0.934, 0.925, 0.917, 0.907, 0.897,
0.887, 0.878, 0.866, 0.855, 0.844,
0.832, 0.806, 0.778, 0.764, 0.749,
0.718, 0.687, 0.654, 0.623, 0.590,
0.559, 0.526, 0.495, 0.462, 0.431,
0.398, 0.367, 0.334, 0.304, 0.272,
0.244, 0.213, 0.187, 0.155, 0.120,
0.090, 0.065, 0.040, 0.015, 0.000

And I would like to make 17 vertically allocated emissions. Is it possible to allocate emissions on different vertical grid?

@Gwang-Jin
Copy link

I know I should not simply interpolate the vertical fractions to my eta_levels. But, I don't know a mathematical keyword to obtain new fractions on the new vertical grid.

@Gwang-Jin
Copy link

Gwang-Jin commented Apr 1, 2022

I've obtained the new fractions on my eta_levels by replacing layerfractions with :

layerfractions = pd.read_csv(io.StringIO("""L,SigmaTop,PPCB,PROT
1,0.995,0.00,0.06
2,0.985,0.066,0.26
3,0.979,0.067,0.34
4,0.973,0.067,0.34
5,0.967,0.1,0.0
6,0.960,0.1,0.0
7,0.954,0.1,0.0
8,0.949,0.05,0.0
9,0.941,0.05,0.0
10,0.934,0.05,0.0
11,0.925,0.05,0.0
12,0.917,0.05,0.0
13,0.907,0.05,0.0
14,0.897,0.05,0.0
15,0.887,0.05,0.0
16,0.878,0.05,0.0
17,0.866,0.05,0.0
"""), comment='#')

I distributed original fraction to new level properly, and the sum of fractions is 1.

I have one question regarding the area source of MICS-Asia RESIDENTIAL, TRANSPORT, and AGRICULTURE.
In [23], the area sources are allocated on the second vertical grids. But, the first vertical grid has also height (35m). So should be the area sources allocated on the first vertical grids?
Thank you.

@barronh
Copy link
Author

barronh commented Apr 7, 2022

I'm glad you figured it out.

You're right about the layers. When it was first made, I think that L was the index. In this version, however, the index is a zero-based value. You should change all of those to layerfractions.loc[0, '...'] = 1. I updated the notebook.

@Gwang-Jin
Copy link

Thank you for response.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment