Post processing of FLEXPART-WRF output

Research

In this post, I present some simple programs written in Python for post-processing the flexpart-wrf output.

It mainly contains several aspects, data merging, data processing and data visualization. I will also show some tips tp creat self-defined colormaps for nice plots.

PS: All th codes are also uploaded in my GitHub respority PyFlex

Merging datasets

In the output path, the result (e.g., flxout_d01_20160304_030000.nc ) for each run was saved as an independent file. It would be inconvenient to loop them for every post-processing function. Therefore, I zip them into a hdf5 file.

Here is the code and some instructions.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
"""pre-reading an template file for capture the geogrid information"""
test_file_path = output_path+"flxout_d01_20160304_030000.nc"
test_file      = Dataset(test_file_path)
pes = np.zeros_like(test_file['/CONC'][:,0,0,0,:,:].sum(axis = 0))
pos=Dataset(output_path+"header_d01.nc")
xx = pos['XLONG'][:]
yy = pos['XLAT'][:]

"""Creat HDF FILES"""
hdf_filename ='Chifeng_PES_72hour.hdf'
os.system("rm "+hdf_filename)
with h5py.File(hdf_filename, "w") as nf:
    dset = nf.create_dataset("test", (100,), dtype='i')
nf = h5py.File(hdf_filename, 'r+')
grp = nf.create_group("PES")
dset1 = grp.create_dataset("lon", data= xx)
dset2 = grp.create_dataset("lat", data= yy)

"""Loop the datasets for each month"""
"""Take the winter months for example"""
print "## WINTER ##"
for mo in ['01','02','12']:
	## The arrays include 6 dimensions, time, ageclass, releases (species), Z, Y, X
    ##  Since most sources are located within the footprint layer (~300 m above the ground), we only calculated the residence time of bottom vertical layer.
    pes = np.zeros_like(test_file['/CONC'][:,0,0,0,:,:].sum(axis = 0))
    t   = 0
    Dir = output_path+"winter/"+mo+"/"
    files = os.listdir(Dir)
    files = sorted(files)
    for file in files:
        filename,extname = os.path.splitext(file)
        if filename[0:4] == 'flxo':
            fc = Dataset(Dir+file)
            if fc['/CONC'].shape[0]>70.0:
                pes+=fc['/CONC'][:,0,0,0,:,:].sum(axis = 0)
                t+=1
    print mo+" "+str(t)
    pes = pes/t
    dset = grp.create_dataset(mo, data= pes)
    print 'Done' 
nf.close()

Then, we can download the “.hdf5” file for further analysis.

Colormap setting

A good colormap can really improve the representability for the figure. Here, I recommend and generate three colormaps using different approaches.

  • cmap1: generated by user-defined RGB values
  • cmap2: subsets of an existing 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
## CMAP1 from pflexible.py on https://git.nilu.no/
from matplotlib.colors import ListedColormap
color_list = [1.0000000e+00, 1.0000000e+00, 1.0000000e+00,
            9.9607843e-01, 9.1372549e-01, 1.0000000e+00,
            9.8431373e-01, 8.2352941e-01, 1.0000000e+00,
            9.6470588e-01, 7.1764706e-01, 1.0000000e+00,
            9.3333333e-01, 6.0000000e-01, 1.0000000e+00,
            8.9019608e-01, 4.4705882e-01, 1.0000000e+00,
            8.3137255e-01, 2.0000000e-01, 1.0000000e+00,
            7.5686275e-01, 0.0000000e+00, 1.0000000e+00,
            6.6274510e-01, 0.0000000e+00, 1.0000000e+00,
            5.4901961e-01, 0.0000000e+00, 1.0000000e+00,
            4.0784314e-01, 0.0000000e+00, 1.0000000e+00,
            2.4705882e-01, 0.0000000e+00, 1.0000000e+00,
            7.4509804e-02, 0.0000000e+00, 1.0000000e+00,
            0.0000000e+00, 2.8235294e-01, 1.0000000e+00,
            0.0000000e+00, 4.8627451e-01, 1.0000000e+00,
            0.0000000e+00, 6.3137255e-01, 1.0000000e+00,
            0.0000000e+00, 7.4509804e-01, 1.0000000e+00,
            0.0000000e+00, 8.4705882e-01, 1.0000000e+00,
            0.0000000e+00, 9.3725490e-01, 1.0000000e+00,
            0.0000000e+00, 1.0000000e+00, 9.7647059e-01,
            0.0000000e+00, 1.0000000e+00, 8.9411765e-01,
            0.0000000e+00, 1.0000000e+00, 8.0000000e-01,
            0.0000000e+00, 1.0000000e+00, 6.9019608e-01,
            0.0000000e+00, 1.0000000e+00, 5.6470588e-01,
            0.0000000e+00, 1.0000000e+00, 4.0000000e-01,
            0.0000000e+00, 1.0000000e+00, 0.0000000e+00,
            3.9607843e-01, 1.0000000e+00, 0.0000000e+00,
            5.6470588e-01, 1.0000000e+00, 0.0000000e+00,
            6.9019608e-01, 1.0000000e+00, 0.0000000e+00,
            7.9607843e-01, 1.0000000e+00, 0.0000000e+00,
            8.9411765e-01, 1.0000000e+00, 0.0000000e+00,
            9.7647059e-01, 1.0000000e+00, 0.0000000e+00,
            1.0000000e+00, 9.4509804e-01, 0.0000000e+00,
            1.0000000e+00, 8.7450980e-01, 0.0000000e+00,
            1.0000000e+00, 7.9215686e-01, 0.0000000e+00,
            1.0000000e+00, 7.0588235e-01, 0.0000000e+00,
            1.0000000e+00, 6.0392157e-01, 0.0000000e+00,
            1.0000000e+00, 4.8235294e-01, 0.0000000e+00,
            1.0000000e+00, 3.1372549e-01, 0.0000000e+00,
            1.0000000e+00, 0.0000000e+00, 1.4901961e-01,
            1.0000000e+00, 0.0000000e+00, 3.3333333e-01,
            1.0000000e+00, 0.0000000e+00, 4.4705882e-01,
            1.0000000e+00, 0.0000000e+00, 5.3725490e-01,
            1.0000000e+00, 0.0000000e+00, 6.1176471e-01,
            9.7647059e-01, 0.0000000e+00, 6.6666667e-01,
            8.9411765e-01, 0.0000000e+00, 6.6666667e-01,
            7.9607843e-01, 0.0000000e+00, 6.3921569e-01,
            6.9019608e-01, 0.0000000e+00, 5.9215686e-01,
            5.6470588e-01, 0.0000000e+00, 5.0980392e-01,
            3.9607843e-01, 0.0000000e+00, 3.8039216e-01]
color_list = np.reshape(color_list, (-1, 3))
name = 'flexpart_cmap'
cmap1 = ListedColormap(color_list, name)
1
2
3
4
5
6
7
8
9
10
## CMAP2 extracted from an existing colormap
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import numpy as np
def truncate_colormap(cmap, minval=0.0, maxval=1.0, n=100):
    new_cmap = colors.LinearSegmentedColormap.from_list(
        'trunc({n},{a:.2f},{b:.2f})'.format(n=cmap.name, a=minval, b=maxval),
        cmap(np.linspace(minval, maxval, n)))
    return new_cmap
cmap2 = truncate_colormap(plt.cm.gist_ncar, 0.3,0.9)

cmap3 is an interesting colormap which I clipped from an existing figure.

I cut the colorbar from the above figure, read its RGB values, and generate a new colormap

1
2
3
4
5
6
7
8
9
10
11
12
13
## CMAP3 from a current colorbar
from PIL import Image
import pandas as pd
im = Image.open('./colorbar_3.png')
rgb_im = im.convert('RGB')
r, g, b = rgb_im.getpixel((1, 1))
k = []
for i in range(18,321,1):
    k.append(rgb_im.getpixel((i,20)))
color_list = np.array(k)    
color_list = np.array([(i/255.0,j/255.0,k/255.0) for i,j,k in color_list])
name  = 'copy_cmap'
cmap3 = ListedColormap(color_list, name)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
# plotting them
from mpl_toolkits.axes_grid1 import make_axes_locatable
def sample_plot(ax,arr,cmap):
    s = ax.imshow(arr, interpolation='nearest', cmap=cmap)
    divider = make_axes_locatable(ax)
    cax = divider.new_vertical(size="5%", pad=0.3, pack_start=True)
    fig.add_axes(cax)
    fig.colorbar(s, cax=cax, orientation="horizontal")    

arr = np.linspace(0, 50, 100).reshape((10, 10))
fig, ax = plt.subplots(ncols=3)
sample_plot(ax[0],arr, cmap1)
sample_plot(ax[1],arr, cmap2)
sample_plot(ax[2],arr, cmap3)

Finally, we can read and visualize the datasets.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
# reading the hdf5 data 
PES_data = h5py.File("/Users/HYF/Downloads/Chifeng_PES_pbl_72hour.hdf")
lon   = PES_data['/PES/lon'][:]
lat   = PES_data['/PES/lat'][:]
# sample arrray of the results in January.
PES_1 = PES_data['/PES/12'][:]
PES_2 = PES_data['/PES/01'][:]
PES_3 = PES_data['/PES/02'][:]
PES_winter = (PES_1+PES_2+PES_12)/3.0     

# plot the winter retroplume with colormap1
fig = plt.figure(figsize=(4,3), frameon=True)
proj = ccrs.LambertConformal(central_latitude = 42.715, 
                             central_longitude = 118.79, 
                             standard_parallels = (30, 60))
ax   =plt.subplot(111, projection = proj)

mask_v = np.ma.masked_less_equal(PES_winter,0)
cs = ax.pcolormesh(lon,lat,mask_v,transform=ccrs.PlateCarree(),cmap =cmap1,alpha = 0.85,zorder=1,norm=matplotlib.colors.LogNorm(),
                      vmin = 1E-9) 
ax.set_extent([105,133,31,52], crs=ccrs.PlateCarree())
ax.coastlines(linewidth = 0.5,resolution='50m')
ax.add_feature(cfeature.BORDERS, linewidth=0.5)

## the xticks, yticks setting is referenced from 
## https://github.com/ARM-DOE/pyart/blob/master/pyart/graph/radarmapdisplay_cartopy.py
## I have not uploaded the original function here to make the code more tight.
fig.canvas.draw()
ax.gridlines(xlocs=xticks, ylocs=yticks, color='gray', alpha=0.5, linestyle='--',linewidth = 0.75)
ax.xaxis.set_major_formatter(LONGITUDE_FORMATTER) 
ax.yaxis.set_major_formatter(LATITUDE_FORMATTER)
lambert_xticks(ax, xticks)
lambert_yticks(ax, yticks) 


pos1 = ax.get_position()
tax = fig.add_axes([pos1.x0,pos1.y1,pos1.x1-pos1.x0,0.03])#x0,y0,long,width
tax.get_xaxis().set_visible(False)
tax.get_yaxis().set_visible(False)
tax.set_facecolor('#FFE5CE')
tax.text(0.4,0.3,"Winter",color = 'k',fontsize =9,fontweight = 'bold',transform=tax.transAxes)

#Python #Air quality #Model
Plotting topography map with hillside 气象模拟文件的后处理PartA-曲折的安装

Kommentare

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×