Skip to content

Instantly share code, notes, and snippets.

@barronh
Last active March 17, 2023 17:48
Show Gist options
  • Star 3 You must be signed in to star a gist
  • Fork 3 You must be signed in to fork a gist
  • 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
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
@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