3

我正在尝试计算相对涡量，即 dV/dX - dU/dY，我目前正在 numpy 上使用梯度函数。下面是我的代码。我想知道是否有更好的方法来做到这一点，而不是在我想做 dU/dY 时尝试重塑数组。给定两个带有数字的矩阵只说 U 和 Y 并且我想对 U wrt Y 进行微分，有没有更好的方法来进行微分。

import numpy as np
import netCDF4
import matplotlib.pyplot as plt
from numpy import *
import decimal
from netCDF4 import Dataset

ncfile= Dataset('test.nc','r')

#--------------------Reading in Variables---------------------------------#

lon      = ncfile.variables['lon'][:]
lat      = ncfile.variables['lat'][:]
UWind850 = ncfile.variables['U'][:,22,:,:] (time, level,lat,lon)
VWind850 = ncfile.variables['V'][:,22,:,:] (time, level,lat,lon)
time     = ncfile.variables['time'][:] 
MSLP     = ncfile.variables['PSL'][:]

# Variable[time,Longitude,Latitude]
#These values are equivalent to I,J and L in the netCDF file

t = 30  #time
x = 300 #longitude 
y = 240 #latitude


#-----------------------Calculating Vorticity-----------------------------#

dX = np.gradient(lon) #shape 300
dY = np.gradient(lat) #shape 240

#VWind850.shape (30,240,300)
#UWind850.shape (30,240,300)

dV = (np.gradient(VWind850))

#dV.shape(3,30,240,300) --The extra "3" dimension is caused by the gradient because 
#Its creating a Matrix for gradients by (time,latitude,longitude)

Vgradient =  dV[2]/dX



UWindTemp = np.reshape(UWind850,(30,300,240)) # I am reshaping so I can divide by dY



dU = (np.gradient(UWindTemp))
Ugradient =  dU[2]/dY
Ugradient = np.reshape(Ugradient,(30,240,300)) # Taking it back to normal


VORT= Vgradient - Ugradient 

# VORT.shape(time, latitude, longitude)
#-------------------------------------------------------------------------#
4

0 回答 0