Modules to load

In [1]:
from __future__ import absolute_import
import numpy as np
import scipy as sp
from netCDF4 import Dataset
import netCDF4
import os.path
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
#from matplotlib.ticker import LinearLocator, FormatStrFormatter
import numpy as np
from mpl_toolkits.basemap import Basemap
import seaborn as sns; sns.set()
import pandas as pd
import matplotlib as mpl

%matplotlib inline

Simple map

In [3]:
#Netcdf file 
path_nc  = '../../OLYMPUS/DEV/trash/DATA/PYPOSTPROC/pm_2009_yravg.nc'
#Variable 
varname = 'O3' 
#conversion
ppb_micro = 1.88 if (varname =='NO2') else 2.00 if (varname == 'O3') else 1.25 if (varname == 'NO') else 1.145 if (varname == 'CO') else 1.00 
#Vertival level 
level = 1 

nc  = path_nc 
fh       = Dataset(nc, mode='r')
var      = fh.variables[str(varname)][:,:]
var = var * ppb_micro
lons     = fh.variables['lon'][:,:]
lats     = fh.variables['lat'][:,:]
fh.close()

#Calculate the average
var_mean    = var[:,int(level),:,:].mean(0) 
var_mean    = np.ma.masked_where(var_mean < 0, var_mean)
lon_0 = lons.mean()
lat_0 = lats.mean()
lon_min =lons.min()
lon_max =lons.max()
lat_min =lats.min()
lat_max =lats.max()
plt.figure(figsize=(24,12))
m = Basemap(resolution='l',projection='tmerc',lat_0=lat_0,lon_0=lon_0,llcrnrlon=lon_min,llcrnrlat= lat_min, urcrnrlon=lon_max, urcrnrlat=lat_max)
m.drawmapboundary(fill_color='#F2F2F2')
xi, yi = m(lons, lats)
#cs = m.pcolor(xi,yi,np.squeeze(var_mean),cmap=cm.viridis,vmin=1, vmax=var.max())
cs = m.pcolor(xi,yi,np.squeeze(var_mean),cmap=cm.gist_ncar_r,vmin=1, vmax=120)
#cs = m.pcolor(xi,yi,np.squeeze(varemis_mean),cmap=cm.winter_r)
m.drawparallels(np.arange(-80,81,0.5),labels=[1,0,0,0])
m.drawmeridians(np.arange(0,360,1),labels=[0,0,0,1])
m.drawmapboundary(fill_color='white')
m.readshapefile('../../OLYMPUS/DEV/trash/DATA/PYPOSTPROC/FRA_adm_shp/FRA_adm1','areas',drawbounds=True)
cbar = m.colorbar(cs, location='right', pad="10%")

Multiple plots

In [7]:
#define a function
def plot_nc(path, lons, lats, varname, timestep) :
    level = 1
    nc  = path 
    fh       = Dataset(nc, mode='r')
    var      = fh.variables[str(varname)][:,:]
    var = var * ppb_micro
    fh.close()
    var_m    = var[timestep,int(level),:,:] #.mean(0) 
    var_m    = np.ma.masked_where(var_m < 0, var_m)    
    m = Basemap(resolution='l',projection='tmerc',lat_0=lat_0,lon_0=lon_0,llcrnrlon=lon_min,llcrnrlat= lat_min, urcrnrlon=lon_max, urcrnrlat=lat_max)
    m.drawmapboundary(fill_color='#F2F2F2')
    xi, yi = m(lons, lats)
    cs = m.pcolor(xi,yi,np.squeeze(var_m),cmap=cm.gist_ncar_r,vmin=1, vmax=var.max())
    m.drawparallels(np.arange(-80,81,0.5),labels=[1,0,0,0])
    m.drawmeridians(np.arange(0,360,1),labels=[0,0,0,1])
    m.drawmapboundary(fill_color='white')
    m.readshapefile('../../OLYMPUS/DEV/trash/DATA/PYPOSTPROC/FRA_adm_shp/FRA_adm1','areas',drawbounds=True)
    cbar = m.colorbar(cs, location='right', pad="10%")
    
    
path = '../../OLYMPUS/DEV/trash/DATA/PYPOSTPROC/pm_2009_seasavg.nc'
varname = 'O3'

plt.figure(1,figsize=(24,12))
plt.subplot(221)
plot_nc(path, lons, lats, varname, 0)
plt.title('Winter',fontsize=16)

plt.subplot(222)
plot_nc(path, lons, lats, varname, 1)
plt.title('Spring',fontsize=16)

plt.figure(2,figsize=(24,12))
plt.subplot(221)
plot_nc(path, lons, lats, varname, 2)
plt.title('Summer',fontsize=16)

plt.subplot(222)
plot_nc(path, lons, lats, varname, 3)
plt.title('Autumn',fontsize=16)
Out[7]:
Text(0.5,1,'Autumn')

Difference Map

In [6]:
path_nc  = '../../OLYMPUS/DEV/trash/DATA/PYPOSTPROC/nox_2009_yravg.nc'
path_nc2 = '../../OLYMPUS/DEV/trash/DATA/PYPOSTPROC/oxy_ref_2009_yravg.nc'

varname = 'NO2'
ppb_micro = 1.88 if (varname =='NO2') else 2.00 if (varname == 'O3') else 1.25 if (varname == 'NO') else 1.145 if (varname == 'CO') else 1.00
level = 1

nc  = path_nc 
fh       = Dataset(nc, mode='r')
var      = fh.variables[str(varname)][:,:]
var = var * ppb_micro
lons     = fh.variables['lon'][:,:]
lats     = fh.variables['lat'][:,:]
fh.close()
nc2  = path_nc2 
fh2       = Dataset(nc2, mode='r')
var2      = fh2.variables[str(varname)][:,:]
var2 = var2 * ppb_micro
fh2.close()
var_mean    = var[:,int(level),:,:].mean(0) - var2[:,int(level),:,:].mean(0) 
var_mean    = np.ma.masked_where(var_mean == 0, var_mean)
lon_0 = lons.mean()
lat_0 = lats.mean()
lon_min =lons.min()
lon_max =lons.max()
lat_min =lats.min()
lat_max =lats.max()
plt.figure(figsize=(24,12))
m = Basemap(resolution='l',projection='tmerc',lat_0=lat_0,lon_0=lon_0,llcrnrlon=lon_min,llcrnrlat= lat_min, urcrnrlon=lon_max, urcrnrlat=lat_max)
m.drawmapboundary(fill_color='#F2F2F2')
xi, yi = m(lons, lats)
#cs = m.pcolor(xi,yi,np.squeeze(var_mean),cmap=cm.viridis,vmin=1, vmax=var.max())
cs = m.pcolor(xi,yi,np.squeeze(var_mean),cmap=cm.seismic,vmin=-10, vmax=10)
#cs = m.pcolor(xi,yi,np.squeeze(varemis_mean),cmap=cm.winter_r)
m.drawparallels(np.arange(-80,81,0.5),labels=[1,0,0,0])
m.drawmeridians(np.arange(0,360,1),labels=[0,0,0,1])
m.drawmapboundary(fill_color='white')
m.readshapefile('../../OLYMPUS/DEV/trash/DATA/PYPOSTPROC/FRA_adm_shp/FRA_adm1','areas',drawbounds=True)

cbar = m.colorbar(cs, location='right', pad="10%")

3D Plot

In [7]:
#Netcdf file 
path_nc  = '../../OLYMPUS/DEV/trash/DATA/PYPOSTPROC/pm_2009_yravg.nc'
#Variable 
varname = 'O3' 
#conversion
ppb_micro = 1.88 if (varname =='NO2') else 2.00 if (varname == 'O3') else 1.25 if (varname == 'NO') else 1.145 if (varname == 'CO') else 1.00 
#Vertival level 
level = 1 

nc  = path_nc 
fh       = Dataset(nc, mode='r')
var      = fh.variables[str(varname)][:,:]
var = var * ppb_micro
lons     = fh.variables['lon'][:,:]
lats     = fh.variables['lat'][:,:]
fh.close()

var_mean    = var[:,int(level),:,:].mean(0) 
X, Y = lons, lats
X, Y = -lats, lons
Z = var_mean #np.log(var +1) #multivariate_gaussian(pos, mu, Sigma)

# Create a surface plot and projected filled contour plot under it.
fig = plt.figure(figsize=(24,20))
ax = fig.gca(projection='3d')
ax.plot_surface(X, Y, Z, rstride=3, cstride=3, linewidth=1, antialiased=True, cmap=cm.viridis)
cset = ax.contourf(X, Y, Z, zdir='z', offset=10, cmap=cm.viridis, alpha=0.5)
# Get rid of the panes
ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))

# Get rid of the spines
ax.w_xaxis.line.set_color((1.0, 1.0, 1.0, 0.0))
ax.w_yaxis.line.set_color((1.0, 1.0, 1.0, 0.0))
ax.w_zaxis.line.set_color((1.0, 1.0, 1.0, 0.0))

# Get rid of the ticks
ax.set_xticks([]) 
ax.set_yticks([]) 
ax.set_zticks([])

# Add the labels
ax.set_xlabel('Lattitude')
ax.set_ylabel('Longitude')
#ax.set_zlabel('Concentration')

ax.set_zlim(0,70)

#Change the view
#ax.view_init(27, -21)
plt.colorbar(cset,shrink=0.5, aspect=15)

plt.show()

3D difference Plot

In [5]:
path_nc  = '../../OLYMPUS/DEV/trash/DATA/PYPOSTPROC/nox_2009_yravg.nc'
path_nc2 = '../../OLYMPUS/DEV/trash/DATA/PYPOSTPROC/oxy_ref_2009_yravg.nc'

varname = 'NO2'
ppb_micro = 1.88 if (varname =='NO2') else 2.00 if (varname == 'O3') else 1.25 if (varname == 'NO') else 1.145 if (varname == 'CO') else 1.00
level = 1

nc  = path_nc 
fh       = Dataset(nc, mode='r')
var      = fh.variables[str(varname)][:,:]
var = var * ppb_micro
lons     = fh.variables['lon'][:,:]
lats     = fh.variables['lat'][:,:]
fh.close()
nc2  = path_nc2 
fh2       = Dataset(nc2, mode='r')
var2      = fh2.variables[str(varname)][:,:]
var2 = var2 * ppb_micro
fh2.close()
var_mean    = var[:,int(level),:,:].mean(0) - var2[:,int(level),:,:].mean(0) 
#var_mean    = np.ma.masked_where(var_mean == 0, var_mean)
lon_0 = lons.mean()
lat_0 = lats.mean()
lon_min =lons.min()
lon_max =lons.max()
lat_min =lats.min()
lat_max =lats.max()
plt.figure(figsize=(24,12))
X, Y = lons, lats
X, Y = -lats, lons
Z = var_mean #np.log(var +1) #multivariate_gaussian(pos, mu, Sigma)

# Create a surface plot and projected filled contour plot under it.
fig = plt.figure(figsize=(24,20))
ax = fig.gca(projection='3d')
ax.plot_surface(X, Y, Z, rstride=3, cstride=3, linewidth=1, antialiased=True, cmap=cm.seismic)
cset = ax.contourf(X, Y, Z, zdir='z', offset=10, cmap=cm.seismic, alpha=0.5)
# Get rid of the panes
ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))
ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 0.0))

# Get rid of the spines
ax.w_xaxis.line.set_color((1.0, 1.0, 1.0, 0.0))
ax.w_yaxis.line.set_color((1.0, 1.0, 1.0, 0.0))
ax.w_zaxis.line.set_color((1.0, 1.0, 1.0, 0.0))

# Get rid of the ticks
ax.set_xticks([]) 
ax.set_yticks([]) 
ax.set_zticks([])

# Add the labels
ax.set_xlabel('Lattitude')
ax.set_ylabel('Longitude')
#ax.set_zlabel('Concentration')

#ax.set_zlim(0,70)

#Change the view
ax.view_init(15, -21)
plt.colorbar(cset,shrink=0.5, aspect=15)

plt.show()
<Figure size 1728x864 with 0 Axes>