Skip to content

Commit b78811a

Browse files
evanmasonevanmason
committed
Fixed ZERO_CROSSING capability for RossbyWaveSpeed object
1 parent 14c5645 commit b78811a

File tree

1 file changed

+94
-75
lines changed

1 file changed

+94
-75
lines changed

‎make_eddy_tracker_list_obj.py‎

Lines changed: 94 additions & 75 deletions
Original file line numberDiff line numberDiff line change
@@ -487,6 +487,11 @@ def create_netcdf(self, directory, savedir, title,
487487
nc.j0 = np.int32(self.j0)
488488
nc.j1 = np.int32(self.j1)
489489

490+
#nc.LONMIN = grd.LONMIN
491+
#nc.LONMAX = grd.LONMAX
492+
#nc.LATMIN = grd.LATMIN
493+
#nc.LATMAX = grd.LATMAX
494+
490495
if 'ROMS' in self.DATATYPE:
491496
nc.ROMS_grid = grd.grdfile
492497
nc.model = model
@@ -512,55 +517,56 @@ def create_netcdf(self, directory, savedir, title,
512517
# Create variables
513518
nc.createVariable('track', np.int32, ('Nobs'), fill_value=self.fillval)
514519
nc.createVariable('n', np.int32, ('Nobs'), fill_value=self.fillval)
520+
515521
# Use of jday should depend on clim vs interann
516522
if self.INTERANNUAL: # AVISO or INTERANNUAL ROMS solution
517523
nc.createVariable('j1', np.int32, ('Nobs'), fill_value=self.fillval)
524+
518525
else: # climatological ROMS solution
519526
nc.createVariable('ocean_time', 'f8', ('Nobs'), fill_value=self.fillval)
527+
520528
nc.createVariable('lon', 'f4', ('Nobs'), fill_value=self.fillval)
521529
nc.createVariable('lat', 'f4', ('Nobs'), fill_value=self.fillval)
522530
nc.createVariable('A', 'f4', ('Nobs'), fill_value=self.fillval)
523531
nc.createVariable('L', 'f4', ('Nobs'), fill_value=self.fillval)
524532
nc.createVariable('U', 'f4', ('Nobs'), fill_value=self.fillval)
525533
nc.createVariable('Teke', 'f4', ('Nobs'), fill_value=self.fillval)
526534
nc.createVariable('radius_e', 'f4', ('Nobs'), fill_value=self.fillval)
535+
527536
if 'Q' in self.DIAGNOSTIC_TYPE:
528537
nc.createVariable('qparameter', 'f4', ('Nobs'), fill_value=self.fillval)
538+
529539
if 'ROMS' in self.DATATYPE:
530540
nc.createVariable('temp', 'f4', ('Nobs'), fill_value=self.fillval)
531541
nc.createVariable('salt', 'f4', ('Nobs'), fill_value=self.fillval)
532542

533-
#nc.createVariable('NLP', 'f8', ('Nobs'), fill_value=self.fillval)
534-
#nc.createVariable('bounds', np.int16, ('Nobs','four'), fill_value=self.fillval)
535-
nc.createVariable('eddy_duration', np.int16, ('Nobs'), fill_value=self.fillval)
543+
#nc.createVariable('eddy_duration', np.int16, ('Nobs'), fill_value=self.fillval)
536544

537545
# Meta data for variables
538546
nc.variables['track'].units = 'ordinal'
539547
nc.variables['track'].min_val = np.int32(0)
548+
nc.variables['track'].max_val = np.int32(0)
540549
nc.variables['track'].long_name = 'track number'
541550
nc.variables['track'].description = 'eddy identification number'
542551

543552
nc.variables['n'].units = 'ordinal'
544-
#nc.variables['n'].min_val = 4
545-
#nc.variables['n'].max_val = 293
546553
nc.variables['n'].long_name = 'observation number'
547554
# Add option here to set length of intervals.....
548555
nc.variables['n'].description = 'observation sequence number (XX day intervals)'
549556

550557
## Use of jday should depend on clim vs interann
551558
if self.INTERANNUAL: # AVISO or INTERANNUAL ROMS solution
552559
nc.variables['j1'].units = 'days'
553-
#nc.variables['j1'].min_val = 2448910
554-
#nc.variables['j1'].max_val = 2455560
555560
nc.variables['j1'].long_name = 'Julian date'
556561
nc.variables['j1'].description = 'date of this observation'
557562
nc.variables['j1'].reference = self.JDAY_REFERENCE
558563
nc.variables['j1'].reference_description = "".join(('Julian '
559564
'date on Jan 1, 1992'))
565+
560566
else: # climatological ROMS solution
561567
nc.variables['ocean_time'].units = 'ROMS ocean_time (seconds)'
562568

563-
nc.variables['eddy_duration'].units = 'days'
569+
#nc.variables['eddy_duration'].units = 'days'
564570
nc.variables['lon'].units = 'deg. longitude'
565571
nc.variables['lon'].min_val = self.LONMIN
566572
nc.variables['lon'].max_val = self.LONMAX
@@ -570,8 +576,10 @@ def create_netcdf(self, directory, savedir, title,
570576

571577
if 'Q' in self.DIAGNOSTIC_TYPE:
572578
nc.variables['A'].units = 'None, normalised vorticity (abs(xi)/f)'
579+
573580
elif 'SLA' in self.DIAGNOSTIC_TYPE:
574581
nc.variables['A'].units = 'cm'
582+
575583
nc.variables['A'].min_val = self.AMPMIN
576584
nc.variables['A'].max_val = self.AMPMAX
577585
nc.variables['A'].long_name = 'amplitude'
@@ -608,8 +616,6 @@ def create_netcdf(self, directory, savedir, title,
608616

609617
if 'Q' in self.DIAGNOSTIC_TYPE:
610618
nc.variables['qparameter'].units = 's^{-2}'
611-
#nc.variables['NLP'].units = 'None, swirl vel. / propagation vel.'
612-
#nc.variables['NLP'].long_name = 'Non-linear parameter'
613619

614620
if 'ROMS' in self.DATATYPE:
615621
nc.variables['temp'].units = 'deg. C'
@@ -681,7 +687,7 @@ def write2netcdf(self, rtime):
681687
already written tracks.
682688
Each inactive track is 'emptied' after saving
683689
"""
684-
tracks2save = np.asarray(self.get_inactive_tracks(rtime))
690+
tracks2save = np.array([self.get_inactive_tracks(rtime)])
685691
DBR = self.DAYS_BTWN_RECORDS
686692

687693
if np.any(tracks2save): # Note, this could break if all eddies become inactive at same time
@@ -695,44 +701,65 @@ def write2netcdf(self, rtime):
695701
(np.all(self.tracklist[i].ocean_time)):
696702

697703
tsize = len(self.tracklist[i].lon)
704+
698705
if (tsize >= self.TRACK_DURATION_MIN / DBR) and tsize > 1.:
699-
706+
lon = np.array([self.tracklist[i].lon])
707+
lat = np.array([self.tracklist[i].lat])
708+
amp = np.array([self.tracklist[i].amplitude])
709+
uavg = np.array([self.tracklist[i].uavg]) * 100. # to cm/s
710+
teke = np.array([self.tracklist[i].teke])
711+
radius_s = np.array([self.tracklist[i].radius_s]) * 1e-3 # to km
712+
radius_e = np.array([self.tracklist[i].radius_e]) * 1e-3 # to km
713+
n = np.arange(tsize, dtype=np.int32)
714+
track = np.full(tsize, self.ch_index)
715+
track_max_val = np.array([nc.variables['track'].max_val,
716+
np.int32(self.ch_index)]).max()
717+
#print self.tracklist[i].ocean_time
718+
#exit()
719+
#eddy_duration = np.array([self.tracklist[i].ocean_time]).ptp()
720+
700721
tend = self.ncind + tsize
701-
nc.variables['lon'][self.ncind:tend] = np.asarray(self.tracklist[i].lon)
702-
nc.variables['lat'][self.ncind:tend] = np.asarray([self.tracklist[i].lat])
703-
nc.variables['A'][self.ncind:tend] = np.array([self.tracklist[i].amplitude])
704-
self.tracklist[i].uavg *= np.array(100.) # to cm/s
705-
nc.variables['U'][self.ncind:tend] = self.tracklist[i].uavg
706-
nc.variables['Teke'][self.ncind:tend] = np.array([self.tracklist[i].teke])
707-
self.tracklist[i].radius_s *= np.array(1e-3) # to km
708-
nc.variables['L'][self.ncind:tend] = self.tracklist[i].radius_s
709-
self.tracklist[i].radius_e *= np.array(1e-3) # to km
710-
nc.variables['radius_e'][self.ncind:tend] = self.tracklist[i].radius_e
722+
nc.variables['lon'][self.ncind:tend] = lon
723+
nc.variables['lat'][self.ncind:tend] = lat
724+
nc.variables['A'][self.ncind:tend] = amp
725+
nc.variables['U'][self.ncind:tend] = uavg
726+
nc.variables['Teke'][self.ncind:tend] = teke
727+
nc.variables['L'][self.ncind:tend] = radius_s
728+
nc.variables['radius_e'][self.ncind:tend] = radius_e
729+
nc.variables['n'][self.ncind:tend] = n
730+
nc.variables['track'][self.ncind:tend] = track
731+
nc.variables['track'].max_val = track_max_val
732+
#nc.variables['eddy_duration'][self.ncind:tend] = eddy_duration
733+
711734
if 'ROMS' in self.DATATYPE:
712-
nc.variables['temp'][self.ncind:tend] = np.array([self.tracklist[i].temp])
713-
nc.variables['salt'][self.ncind:tend] = np.array([self.tracklist[i].salt])
714-
#nc.variables['bounds'][self.ncind:tend] = np.array([self.tracklist[i].bounds])
735+
temp = np.array([self.tracklist[i].temp])
736+
salt = np.array([self.tracklist[i].salt])
737+
738+
nc.variables['temp'][self.ncind:tend] = temp
739+
nc.variables['salt'][self.ncind:tend] = salt
740+
715741
if self.INTERANNUAL:
716742
# We add 1 because 'j1' is an integer in ncsavefile; julian day midnight has .5
717743
# i.e., dt.julian2num(2448909.5) -> 727485.0
718-
nc.variables['j1'][self.ncind:tend] = dt.num2julian(np.array([self.tracklist[i].ocean_time])) + 1
744+
j1 = dt.num2julian(np.array([self.tracklist[i].ocean_time])) + 1
745+
nc.variables['j1'][self.ncind:tend] = j1
746+
719747
else:
720-
nc.variables['ocean_time'][self.ncind:tend] = np.array([self.tracklist[i].ocean_time])
721-
nc.variables['n'][self.ncind:tend] = np.arange(tsize, dtype=np.int32)
722-
nc.variables['track'][self.ncind:tend] = np.full(tsize, self.ch_index)
723-
nc.variables['track'].max_val = np.int32(self.ch_index)
724-
nc.variables['eddy_duration'][self.ncind:tend] = np.array([self.tracklist[i].ocean_time]).size \
725-
* self.DAYS_BTWN_RECORDS
748+
ocean_time = np.array([self.tracklist[i].ocean_time])
749+
nc.variables['ocean_time'][self.ncind:tend] = ocean_time
750+
726751
if self.TRACK_EXTRA_VARIABLES:
727-
nc.variables['shape_error'][self.ncind:tend] = np.array([self.tracklist[i].shape_error])
752+
shape_error = np.array([self.tracklist[i].shape_error])
753+
nc.variables['shape_error'][self.ncind:tend] = shape_error
754+
728755
for j in np.arange(tend - self.ncind):
729756
jj = j + self.ncind
730-
contour_e_arr = np.asarray(self.tracklist[i].contour_e[j]).ravel()
757+
contour_e_arr = np.array([self.tracklist[i].contour_e[j]]).ravel()
758+
contour_s_arr = np.array([self.tracklist[i].contour_s[j]]).ravel()
759+
uavg_profile_arr = np.array([self.tracklist[i].uavg_profile[j]]).ravel()
731760
nc.variables['contour_e'][:contour_e_arr.size, jj] = contour_e_arr
732-
contour_s_arr = np.asarray(self.tracklist[i].contour_s[j]).ravel()
733761
nc.variables['contour_s'][:contour_s_arr.size, jj] = contour_s_arr
734-
uavg_profile_arr = np.asarray(self.tracklist[i].uavg_profile[j]).ravel()
735-
nc.variables['uavg_profile'][:uavg_profile_arr.size, jj] = uavg_profile_arr #np.asarray(self.tracklist[i].uavg_profile[j]).ravel()
762+
nc.variables['uavg_profile'][:uavg_profile_arr.size, jj] = uavg_profile_arr
736763

737764
# Flag indicating track[i] is now saved
738765
self.tracklist[i].saved2nc = True
@@ -743,6 +770,7 @@ def write2netcdf(self, rtime):
743770
# Get index to first currently active track
744771
try:
745772
lasti = self.get_active_tracks(rtime)[0]
773+
746774
except Exception:
747775
lasti = None
748776

@@ -921,35 +949,7 @@ def insert_at_index(self, xarr, ind, x):
921949

922950
return self
923951

924-
925-
926-
927-
#def set_bounds(self, centlon, centlat, radius, i, j, grd):
928-
#"""
929-
#Get indices to a bounding box around the eddy
930-
#"""
931-
#def get_angle(deg, ang):
932-
#return deg - np.rad2deg(ang)
933-
934-
#grdangle = grd.angle()[j,i]
935-
##print type(centlon), type(centlat)
936-
#a_lon, a_lat = newPosition(centlon, centlat, get_angle(0, grdangle), radius)
937-
#b_lon, b_lat = newPosition(centlon, centlat, get_angle(90, grdangle), radius)
938-
#c_lon, c_lat = newPosition(centlon, centlat, get_angle(180, grdangle), radius)
939-
#d_lon, d_lat = newPosition(centlon, centlat, get_angle(270, grdangle), radius)
940-
941-
## Get i,j of bounding box around eddy
942-
##print grd.lon().shape, grd.lat().shape, grd.shape
943-
#a_i, a_j = nearest(a_lon, a_lat, grd.lon(), grd.lat(), grd.shape)
944-
#b_i, b_j = nearest(b_lon, b_lat, grd.lon(), grd.lat(), grd.shape)
945-
#c_i, c_j = nearest(c_lon, c_lat, grd.lon(), grd.lat(), grd.shape)
946-
#d_i, d_j = nearest(d_lon, d_lat, grd.lon(), grd.lat(), grd.shape)
947-
948-
#self.imin = np.maximum(np.min([a_i, b_i, c_i, d_i]) - 5, 0) # Must not go
949-
#self.jmin = np.maximum(np.min([a_j, b_j, c_j, d_j]) - 5, 0) # below zero
950-
#self.imax = np.max([a_i, b_i, c_i, d_i]) + 5
951-
#self.jmax = np.max([a_j, b_j, c_j, d_j]) + 5
952-
#return self
952+
953953

954954
def set_bounds(self, contlon, contlat, grd):
955955
"""
@@ -1037,10 +1037,13 @@ def get_rwdistance(self, xpt, ypt, DAYS_BTWN_RECORDS):
10371037
self.distance *= 86400.
10381038
#if self.THE_DOMAIN in 'ROMS':
10391039
#self.distance *= 1.5
1040+
10401041
elif 'BlackSea' in self.THE_DOMAIN:
10411042
self.distance[:] = 15000. # e.g., Blokhina & Afanasyev, 2003
1043+
10421044
elif 'MedSea' in self.THE_DOMAIN:
10431045
self.distance[:] = 20000.
1046+
10441047
else:
10451048
Exception # Unknown THE_DOMAIN
10461049

@@ -1060,7 +1063,7 @@ def get_rwdistance(self, xpt, ypt, DAYS_BTWN_RECORDS):
10601063
'eddy at %s, %s in the %s domain'
10611064
% (lon, lat, self.THE_DOMAIN)))
10621065
c = np.abs(self._get_rlongwave_spd(xpt, ypt))[0]
1063-
print "".join(('--------- with extratropical long baroclinic',
1066+
print "".join(('--------- with extratropical long baroclinic ',
10641067
'Rossby wave phase speed of %s m/s' % c))
10651068
self.start = False
10661069

@@ -1073,8 +1076,14 @@ def _make_subset(self):
10731076
"""
10741077
pad = 1.5 # degrees
10751078
LONMIN, LONMAX, LATMIN, LATMAX = self.limits
1076-
lloi = np.logical_and(self._lon >= LONMIN - pad,
1077-
self._lon <= LONMAX + pad)
1079+
if self.ZERO_CROSSING:
1080+
ieast, iwest = (((self._lon + 360.) <= LONMAX + pad),
1081+
(self._lon > LONMIN + pad))
1082+
self._lon[ieast] += 360.
1083+
lloi = iwest + ieast
1084+
else:
1085+
lloi = np.logical_and(self._lon >= LONMIN - pad,
1086+
self._lon <= LONMAX + pad)
10781087
lloi *= np.logical_and(self._lat >= LATMIN - pad,
10791088
self._lat <= LATMAX + pad)
10801089
self._lon = self._lon[lloi]
@@ -1085,9 +1094,9 @@ def _make_subset(self):
10851094

10861095

10871096
def _make_kdtree(self):
1088-
#points = np.vstack([self._lon, self._lat]).T
10891097
points = np.vstack([self.x, self.y]).T
10901098
self._tree = spatial.cKDTree(points)
1099+
return self
10911100

10921101

10931102
def _get_defrad(self, xpt, ypt):
@@ -1104,7 +1113,7 @@ def _get_defrad(self, xpt, ypt):
11041113

11051114
def _get_rlongwave_spd(self, xpt, ypt):
11061115
"""
1107-
Get the longwave phase speed, see Chelton etal (2008) pg 446:
1116+
Get the longwave phase speed, see Chelton etal (1998) pg 446:
11081117
c = -beta * defrad**2 (this only for extratropical waves...)
11091118
"""
11101119
self.r_spd_long[:] = self._get_defrad(xpt, ypt)
@@ -1157,7 +1166,7 @@ def __init__(self, THE_DOMAIN, grd, DAYS_BTWN_RECORDS, RW_PATH=None):
11571166

11581167
def _set_east_ellipse(self):
11591168
"""
1160-
The east_ellipse is a full ellipse, but only its eastern
1169+
The *east_ellipse* is a full ellipse, but only its eastern
11611170
part is used to build the search ellipse.
11621171
"""
11631172
self.east_ellipse = patch.Ellipse((self.xpt, self.ypt),
@@ -1167,7 +1176,7 @@ def _set_east_ellipse(self):
11671176

11681177
def _set_west_ellipse(self):
11691178
"""
1170-
The east_ellipse is a full ellipse, but only its eastern
1179+
The *west_ellipse* is a full ellipse, but only its western
11711180
part is used to build the search ellipse.
11721181
"""
11731182
self.west_ellipse = patch.Ellipse((self.xpt, self.ypt),
@@ -1177,7 +1186,8 @@ def _set_west_ellipse(self):
11771186

11781187
def _set_global_ellipse(self):
11791188
"""
1180-
1189+
Set a Path object *ellipse_path* built from the eastern vertices of
1190+
*east_ellipse* and the western vertices of *west_ellipse*.
11811191
"""
11821192
self._set_east_ellipse()._set_west_ellipse()
11831193
e_verts = self.east_ellipse.get_verts()
@@ -1194,6 +1204,7 @@ def _set_global_ellipse(self):
11941204

11951205
def _set_black_sea_ellipse(self):
11961206
"""
1207+
Set *ellipse_path* for the *black_sea_ellipse*.
11971208
"""
11981209
self.black_sea_ellipse = patch.Ellipse((self.x, self.y),
11991210
2. * self.rw_c_mod, 2. * self.rw_c_mod)
@@ -1205,6 +1216,14 @@ def _set_black_sea_ellipse(self):
12051216

12061217
def set_search_ellipse(self, xpt, ypt):
12071218
"""
1219+
Set the search ellipse around a point.
1220+
1221+
args:
1222+
1223+
*xpt*: lon coordinate (Basemap projection)
1224+
1225+
*ypt*: lat coordinate (Basemap projection)
1226+
12081227
"""
12091228
self.xpt = xpt
12101229
self.ypt = ypt
@@ -1219,7 +1238,7 @@ def set_search_ellipse(self, xpt, ypt):
12191238
self.rw_c_mod *= 2.
12201239
self._set_global_ellipse()
12211240

1222-
elif 'BlackSea' in self.THE_DOMAIN:
1241+
elif 'BlackSea' in self.THE_DOMAIN:
12231242
self.rw_c_mod[:] = 1.75
12241243
self.rw_c[:] = self.rwv.get_rwdistance(x, y,
12251244
self.DAYS_BTWN_RECORDS)

0 commit comments

Comments
 (0)