-
-
Save barronh/97633470997ca77b4bb949df20d9a0a3 to your computer and use it in GitHub Desktop.
@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)
@Gwang-Jin - make sure you're aware of this.
@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!
I made the changes to the gist. Can you confirm that I made them where you did?
I made the changes to the gist. Can you confirm that I made them where you did?
Yes. That is correct!
@ 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
Thank you for your updates. With your permission, I would like to add your name at the top as a contributor.
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].
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!
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 ?
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.
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'
@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)
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.
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?
Great. I made the update.
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
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
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.
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.
I got it regarding NOx. Thank you.
@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!
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.
- Convert emissions from CMAQ-ready mass/rates and mole/rates files to flux files (moles/s to moles/m2/s).
- Set the flux file to the projected space using cdo's grid definition file
- 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.
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
-
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 -
make sure you also install other required packages for running this gist.
Pseudonetcdf
pyproj -
Now you don't have to use python-cdo interface, you can comment out like
#import cdo
#cdoo = cdo.Cdo(cdopath) -
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')
- 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')
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?
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.
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.
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.
Thank you for response.
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?