-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathGenerateSEAPODYMFields.py
82 lines (68 loc) · 3.2 KB
/
GenerateSEAPODYMFields.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
from parcels import *
from SEAPODYM_functions import *
from netCDF4 import Dataset
import matplotlib.pyplot as plt
if __name__ == "__main__":
forcingU = 'SEAPODYM_Forcing_Data/PHYSICAL/2003Cohort_PHYS_month1.nc'
forcingV = 'SEAPODYM_Forcing_Data/PHYSICAL/2003Cohort_PHYS_month1.nc'
forcingH = 'SEAPODYM_Forcing_Data/HABITAT/2003Cohort_HABITAT_month59.nc'
filenames = {'U': forcingU, 'V': forcingV, 'H': forcingH}
variables = {'U': 'u', 'V': 'v', 'H': 'habitat'}
dimensions = {'lon': 'lon', 'lat': 'lat', 'time': 'time'}
print("Creating Grid")
grid = Grid.from_netcdf(filenames=filenames, variables=variables, dimensions=dimensions,
vmax=200)
#print("Creating Diffusion Field")
#K = Create_SEAPODYM_Diffusion_Field(grid.H, 30*24*60*60)
#grid.add_field(K)
print("Creating Diffusion Field with SEAPODYM units")
S_K = Create_SEAPODYM_Diffusion_Field(grid.H, 30*24*60*60, units='nm2_per_mon', P=3, start_age=62)
S_K.name = 'S_K'
grid.add_field(S_K)
grid.write('SEAPODYM_Grid_Comparison')
print("Loading Diffusion Field")
bfile = Dataset('DIFFUSION_last_cohort/last_cohort_diffusion.nc', 'r')
bX = bfile.variables['longitude']
bY = bfile.variables['latitude']
bT = bfile.variables['time']
bVar = bfile.variables['skipjack_diffusion_rate']
hfile = Dataset('SEAPODYM_Forcing_Data/SEAPODYM2003_HABITAT_Prepped.nc')
hX = hfile.variables['lon']
hY = hfile.variables['lat']
hT = hfile.variables['time']
hVar = hfile.variables['habitat']
#I_K = Field.from_netcdf('I_K', {'lon':'longitude', 'lat':'latitude', 'time':'time'},
# ['DIFFUSION/INTERIM-NEMO-PISCES_skipjack_diffusion_rate_20030115.nc',
# 'DIFFUSION/INTERIM-NEMO-PISCES_skipjack_diffusion_rate_20030215.nc'])
I_K = Field('I_K', data=bVar, lon=bX, lat=bY, time=bT, transpose=True)
grid.add_field(I_K)
comparison = np.zeros(np.shape(I_K.data), dtype=np.float32)
loaded = np.array(bVar)
comparison = abs(S_K.data-loaded)
print(np.min(comparison))
print(np.max(comparison))
fig = plt.figure(1)
ax = plt.axes()
plt.contourf(bX[:], bY[:], bVar[59, :, :], vmin=0, vmax=80000)
plt.title("Inna's K (final timestep)")
plt.colorbar()
fig = plt.figure(2)
ax = plt.axes()
plt.contourf(bX[:], bY[:], S_K.data[59, :, :], vmin=0, vmax=80000)
plt.title("Calculated K (final timestep)")
plt.colorbar()
fig = plt.figure(3)
ax = plt.axes()
plt.contourf(bX[:], bY[:], comparison[59, :, :], vmax=10000)
plt.show()
print(np.shape(hVar))
print(np.shape(bVar))
#example habitat to diffusion calculations
x = [147, 110, 68, 120]
y = [85, 64, 44, 60]
for i in range(len(x)):
print("Final timestep H at %s-%s = %s" % (grid.H.lon[x[i]], grid.H.lat[y[i]], grid.H.data[59, y[i], x[i]]))
print("Final timestep K at %s-%s = %s" % (grid.S_K.lon[x[i]], grid.S_K.lat[y[i]], grid.S_K.data[59, y[i], x[i]]))
print("Final timestep H file at %s-%s = %s" % (hX[x[i]], hY[y[i]], hVar[59, :, y[i], x[i]]))
print("Final timestep Inna K at %s-%s = %s" % (bX[x[i]], bY[y[i]], bVar[59, y[i], x[i]]))
#grid.write('FieldTest')