NCUATM天氣學作業1：溫度平流、輻散、渦度、渦度平流

syn_met_01.py

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import netCDF4 as nc

# Note: Using "bin_to_nc.py" convert to NETCDF format first!!
##### C O N F I G #####
ncfile = 'output.nc'
#######################

#Constant
Re = 6378000 #m
omaga = 7.29E-5
d = 1.875 #deg

data = nc.Dataset(ncfile)

lats = data.variables['lat'][:]
lons = data.variables['lon'][:]
layers = data.variables['h'][:]
lon, lat = np.meshgrid(lons, lats)

HH = data.variables['H'][:,:,:]
UU = data.variables['U'][:,:,:]
VV = data.variables['V'][:,:,:]
TT = data.variables['T'][:,:,:]

dy = np.full(HH[0,:,:].shape, Re * d * np.pi / 180, dtype=np.float64)
dx = dy * np.cos(np.radians(lat))
f = 2 * omaga * np.sin(np.radians(lat))

def diff_x(A):
	ANS = np.zeros(A.shape)
	left = A[:-2,:]
	right = A[2:,:]
	
	#Central difference
	CEN = (right - left) / (2 * dx[1:-1,:])
	
	#Forward difference for left bondary
	LB = (A[1,:] - A[0,:]) / dx[0,:]
	
	#Backward difference for right bondary
	RB = (A[-1,:] - A[-2,:]) / dx[-1,:]
	
	ANS[0,:] = LB
	ANS[1:-1,:] = CEN
	ANS[-1,:] = RB
	
	return ANS
	
def diff_y(A):
	ANS = np.zeros(A.shape)
	up = A[:,:-2]
	dn = A[:,2:]
	
	#Central difference
	CEN = (dn - up) / (2 * dy[:,1:-1])
	
	#Forward difference for lower bondary
	LB = (A[:,1] - A[:,0]) / dy[:,0]
	
	#Backward difference for upper bondary
	UB = (A[:,-1] - A[:,-2]) / dy[:,-1]
	
	ANS[:,0] = LB
	ANS[:,1:-1] = CEN
	ANS[:,-1] = UB
	
	return ANS
	
def plot(A, title):
	map = Basemap(projection="mill", llcrnrlon=lons[0], urcrnrlon=lons[-1], llcrnrlat=lats[0], urcrnrlat=lats[-1], resolution="l")
	map.drawcoastlines()
	map.drawcountries()
	x, y = map(lon, lat)

	map.contourf(x, y, A, cmap='jet')

	cbar = plt.colorbar()
	
	plt.title(title)

	#plt.show()
	plt.savefig('out/' + title + '.png', dpi=300)
	plt.clf()
	plt.close("all")
	print(title)

for layer in range(len(layers)):
	H = HH[layer,:,:]
	U = UU[layer,:,:]
	V = VV[layer,:,:]
	T = TT[layer,:,:]

	ADV_T = (-1 * U * diff_x(T)) - (V * diff_y(T))
	DIV = diff_x(U) + diff_y(V)
	VORT = diff_x(V) - diff_y(U)
	ADV_VORTA = (-1 * U * diff_x(VORT + f)) - (V * diff_y(VORT + f))

	plot(ADV_T, f'{str(int(layers[layer]))}hPa Temperature Advection')
	plot(DIV, f'{str(int(layers[layer]))}hPa Divergence')
	plot(VORT, f'{str(int(layers[layer]))}hPa Relative Vorticity')
	plot(ADV_VORTA, f'{str(int(layers[layer]))}hPa Absolute Vorticity Advection')