-
Notifications
You must be signed in to change notification settings - Fork 0
/
test5.py
82 lines (67 loc) · 2.45 KB
/
test5.py
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
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
import numpy as np
from netCDF4 import Dataset,num2date
import sys
import time
file = 'uwnd_1000_2021_26_8_00Z.nc'
nc_uwndFile = Dataset(file,'r')
s = nc_uwndFile.variables['uwnd'][:]
lats = nc_uwndFile.variables['lat'][:]
lons = nc_uwndFile.variables['lon'][:]
lat = lats[:].squeeze()
lon = lons[:].squeeze()
s = s[0,0,:,:]
lonLen = len(lon)
latLen = len(lat)
dsdy_ref = np.zeros((latLen,lonLen))
dsdy = np.zeros((latLen,lonLen))
rearth = 6371221.3
dj = abs(np.radians((lat[0]-lat[1])) * rearth)
for i in range(0,lonLen):
s0 = s[0,i]
s1 = s[1,i]
if (s0 -999.99 and s1 > -999.99):
# Make sure derivative is dq/dy = [q(north) - q(south)]/dlat
dsdy_ref[0,i] = (s[0,i] - s[1,i])/dj
else:
dsdy_ref[0,i] = -999.99
for i in range(0,lonLen):
if (s[-1,i] > -999.99 and s[-2,i] > -999.99):
# Make sure derivative is dq/dy = [q(north) - q(south)]/dlat
dsdy_ref[-1,i] = (s[-2,i] - s[-1,i])/dj
else:
dsdy_ref[-1,i] = -999.99
for j in range(1,latLen-1):
for i in range(0,lonLen):
if (s[j-1,i] > -999.99 and s[j+1,i] > -999.99):
# Make sure derivative is dq/dy = [q(north) - q(south)]/dlat
dsdy_ref[j,i] = (s[j+1,i] - s[j-1,i])/(2.0*dj)
elif (s[j-1,i] < -999.99 and s[j+1,i] > -999.99 and s[j,i] > -999.99):
dsdy_ref[j,i] = (s[j+1,i] - s[j,i])/dj
elif (s[j-1,i] > -999.99 and s[j+1,i] < -999.99 and s[j,i] > -999.99):
dsdy_ref[j,i] = (s[j,i] - s[j-1,i])/dj
else :
dsdy_ref[j,i] = -999.99
# North Pole
hasNValue = s[0,:] > -999.99
hasNLeft = s[1,:] > -999.99
dsdy[0,:] = -999.99
dsdy[0,:] = np.where(hasNValue & hasNLeft, (s[0,:]-s[1,:])/dj,dsdy[0,:])
#South Pole
hasSRValue = s[-1,:] > -999.99
hasSR2Value = s[-2,:] > -999.99
dsdy[-1,:] = -999.99
dsdy[-1,:] = np.where(hasSRValue & hasSR2Value,(s[-2,:]-s[-1,:])/dj,dsdy[-1,:])
#Regular coordinates
has_value = s[1:-1, :] > -999.99
has_right = s[2:,:] > -999.99
has_left = s[:-2,:] > -999.99
dsdy[1:-1,:] = -999.99
dsdy[1:-1,:] = np.where(has_left & has_value,(s[2,:] - s[1:-1,:]) / dj, dsdy[1:-1,:])
dsdy[1:-1,:] = np.where(has_right & has_value,(s[1:-1,:] - s[:-2,:]) / dj, dsdy[1:-1,:])
dsdy[1:-1,:] = np.where(has_left & has_right,(s[2:,:] - s[:-2,:])/(2.*dj),dsdy[1:-1,:])
print(np.allclose(dsdy_ref,dsdy,1e-14))
#for j in range(0,73):
# for i in range(0,144):
#diff = dsdy_ref[j,i] - dsdy[j,i]
#print(lat[j],dsdy_ref[j,i],dsdy[j,i])
#print(diff)