forked from ecjoliver/stormTracking
-
Notifications
You must be signed in to change notification settings - Fork 2
/
storm_detection.py
121 lines (102 loc) · 4.52 KB
/
storm_detection.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
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
'''
Software for the tracking of storms and high-pressure systems
'''
#
# Load required modules
#
import numpy as np
from datetime import date
from netCDF4 import Dataset
from matplotlib import pyplot as plt
import storm_functions as storm
#
# Load in slp data and lat/lon coordinates
#
dataset = 'NCEP_20CRV2C'
dataset = 'NCEP_R1'
dataset = 'NCEP_CFSR'
# Parameters
pathroot = {'NCEP_20CRV2C': '/home/oliver/data/NCEP/20CRv2c/prmsl/6hourly/', 'NCEP_R1': '/home/oliver/data/NCEP/R1/slp/', 'NCEP_CFSR': '/home/oliver/data/NCEP/CFSR/prmsl/'}
var = {'NCEP_20CRV2C': 'prmsl', 'NCEP_R1': 'slp', 'NCEP_CFSR': 'PRMSL_L101'}
# Generate date and hour vectors
yearStart = {'NCEP_20CRV2C': 1851, 'NCEP_R1': 1948, 'NCEP_CFSR': 1979}
yearEnd = {'NCEP_20CRV2C': 2014, 'NCEP_R1': 2017, 'NCEP_CFSR': 1979}
# Load lat, lon
filename = {'NCEP_20CRV2C': pathroot['NCEP_20CRV2C'] + 'prmsl.' + str(yearStart['NCEP_20CRV2C']) + '.nc',
'NCEP_R1': pathroot['NCEP_R1'] + 'slp.' + str(yearStart['NCEP_R1']) + '.nc',
'NCEP_CFSR': pathroot['NCEP_CFSR'] + 'prmsl.gdas.' + str(yearStart['NCEP_CFSR']) + '01.grb2.nc'}
fileobj = Dataset(filename[dataset], 'r')
lon = fileobj.variables['lon'][:].astype(float)
lat = fileobj.variables['lat'][:].astype(float)
fileobj.close()
# Load slp data
slp = np.zeros((0, len(lat), len(lon)))
year = np.zeros((0,))
month = np.zeros((0,))
day = np.zeros((0,))
hour = np.zeros((0,))
for yr in range(yearStart[dataset], yearEnd[dataset]+1):
if (dataset == 'NCEP_20CRV2C') + (dataset == 'NCEP_R1'):
filename = {'NCEP_20CRV2C': pathroot['NCEP_20CRV2C'] + 'prmsl.' + str(yr) + '.nc', 'NCEP_R1': pathroot['NCEP_R1'] + 'slp.' + str(yr) + '.nc'}
fileobj = Dataset(filename[dataset], 'r')
time = fileobj.variables['time'][:]
time_ordinalDays = time/24. + date(1800,1,1).toordinal()
year = np.append(year, [date.fromordinal(np.floor(time_ordinalDays[tt]).astype(int)).year for tt in range(len(time))])
month = np.append(month, [date.fromordinal(np.floor(time_ordinalDays[tt]).astype(int)).month for tt in range(len(time))])
day = np.append(day, [date.fromordinal(np.floor(time_ordinalDays[tt]).astype(int)).day for tt in range(len(time))])
hour = np.append(hour, (np.mod(time_ordinalDays, 1)*24).astype(int))
slp0 = fileobj.variables[var[dataset]][:].astype(float)
slp = np.append(slp, slp0, axis=0)
fileobj.close()
print yr, slp0.shape[0]
if dataset == 'NCEP_CFSR':
for mth in range(1, 12+1):
filename = {'NCEP_CFSR': pathroot['NCEP_CFSR'] + 'prmsl.gdas.' + str(yr) + str(mth).zfill(2) + '.grb2.nc'}
fileobj = Dataset(filename[dataset], 'r')
time = fileobj.variables['time'][:]
time_ordinalDays = time/24. + date(yr,mth,1).toordinal()
year = np.append(year, [date.fromordinal(np.floor(time_ordinalDays[tt]).astype(int)).year for tt in range(len(time))])
month = np.append(month, [date.fromordinal(np.floor(time_ordinalDays[tt]).astype(int)).month for tt in range(len(time))])
day = np.append(day, [date.fromordinal(np.floor(time_ordinalDays[tt]).astype(int)).day for tt in range(len(time))])
hour = np.append(hour, (np.mod(time_ordinalDays, 1)*24).astype(int))
slp0 = fileobj.variables[var[dataset]][:].astype(float)
slp = np.append(slp, slp0, axis=0)
fileobj.close()
print yr, mth, slp0.shape[0]
#
# Storm Detection
#
# Initialisation
lon_storms_a = []
lat_storms_a = []
amp_storms_a = []
lon_storms_c = []
lat_storms_c = []
amp_storms_c = []
# Loop over time
T = slp.shape[0]
for tt in range(T):
#
print tt, T
#
# Detect lon and lat coordinates of storms
#
lon_storms, lat_storms, amp = storm.detect_storms(slp[tt,:,:], lon, lat, res=2, Npix_min=9, cyc='anticyclonic', globe=True)
lon_storms_a.append(lon_storms)
lat_storms_a.append(lat_storms)
amp_storms_a.append(amp)
#
lon_storms, lat_storms, amp = storm.detect_storms(slp[tt,:,:], lon, lat, res=2, Npix_min=9, cyc='cyclonic', globe=True)
lon_storms_c.append(lon_storms)
lat_storms_c.append(lat_storms)
amp_storms_c.append(amp)
#
# Save as we go
#
if (np.mod(tt, 100) == 0) + (tt == T-1):
print 'Save data...'
#
# Combine storm information from all days into a list, and save
#
storms = storm.storms_list(lon_storms_a, lat_storms_a, amp_storms_a, lon_storms_c, lat_storms_c, amp_storms_c)
np.savez('storm_det_slp', storms=storms, year=year, month=month, day=day, hour=hour)