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]:
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()
SNS Heatmap with values in the mesh¶
In [6]:
nx=[]
ny=[]
for y in range(129):
for x in range(159):
nx.append(x+1)
ny.append(y+1)
# make a dictionary with each column of information
d={}
d['lon'] =nx #lons[:,:].reshape(20511)
d['lat'] = ny #lats[:,:].reshape(20511)
d['conc'] = var_mean[:,:].reshape(20511)
df = pd.DataFrame(d)
#ax = sns.heatmap(toto,annot=True)
df['conc']=df.conc.astype(int)
df_nc = df.pivot("lat", "lon", "conc")
fig = plt.figure(figsize=(200,150))
mpl.rc('font', size=100)
ax = sns.heatmap(df_nc, annot=True, fmt="d")
