High Frequency doppler Radars¶
See background review. Also, see presentation.
In this lecture, we will explore further how to plot vector maps in python.
[1]:
import numpy as np
import scipy.io as sio
import matplotlib.pyplot as plt
import pandas as pd
from datetime import datetime, timedelta
# add mapping routines
import cartopy.crs as ccrs
import cartopy.feature as cf
import matplotlib.ticker as mticker
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
First, we will load the matlab binary file “hfr_2011.mat” using the scipy.io routine (loaded above as sio).
[2]:
# primary data set, south Oahu
hfr = sio.loadmat('/home/ocn463-data/hfr/hfr_2011.mat')
It contains:
Name, Size, Bytes, Class
Like netcdf files, matlab files can be explored by just typing in the file name.
[3]:
hfr
[3]:
{'__header__': b'MATLAB 5.0 MAT-file, Platform: GLNXA64, Created on: Wed Apr 27 11:16:15 2022',
'__version__': '1.0',
'__globals__': [],
'coast': array([[-157.6503, 21.3079],
[-157.6493, 21.3068],
[-157.6493, 21.3054],
...,
[-157.6517, 21.3093],
[-157.6507, 21.3083],
[-157.6503, 21.3079]]),
'T': array([[734504.02083291],
[734504.06249957],
[734504.10416624],
...,
[734868.89583155],
[734868.93749821],
[734868.97916488]]),
'X': array([[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
...,
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan]]),
'Y': array([[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
...,
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan],
[nan, nan, nan, ..., nan, nan, nan]]),
'README': array(['2011 HFDR currents clean and detided, T is time in matlab format, U and V are in m/s and X and Y are lon and lat in degree'],
dtype='<U122'),
'U': array([[[ nan, nan, nan, ..., -2.28903437,
-2.81524868, -3.08926222],
[ nan, nan, nan, ..., -2.2359242 ,
-3.88160247, -3.94921084],
[ nan, nan, nan, ..., -3.1498654 ,
-4.03816528, -4.8022171 ],
...,
[ nan, nan, nan, ..., nan,
nan, -0.37626221],
[ nan, nan, nan, ..., nan,
nan, -0.97906876],
[ nan, nan, nan, ..., nan,
nan, nan]],
[[ nan, nan, nan, ..., -2.55943637,
-3.42399995, -2.74154142],
[-2.41129668, nan, nan, ..., -3.37086532,
-4.44360008, -4.87301323],
[-1.6614966 , 2.53829481, nan, ..., -3.3576725 ,
-3.48378386, -4.23251635],
...,
[ nan, nan, nan, ..., nan,
nan, -0.24761281],
[ nan, nan, nan, ..., nan,
nan, 2.83584583],
[ nan, nan, nan, ..., nan,
nan, -0.79939173]],
[[-0.34274512, -2.24851056, 1.25880269, ..., -3.28266232,
-3.86356129, -3.91432103],
[-0.30126183, -1.24112729, -0.56163281, ..., -3.72570423,
-3.91486649, -3.91774765],
[-1.03515178, -1.75623236, -0.88225632, ..., -4.32275363,
-3.95083582, -4.85862864],
...,
[ nan, nan, nan, ..., nan,
nan, -0.17602402],
[ nan, nan, nan, ..., nan,
nan, 0.31123404],
[ nan, nan, nan, ..., nan,
nan, -0.83010557]],
...,
[[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
...,
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]],
[[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
...,
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]],
[[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
...,
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]]]),
'V': array([[[ nan, nan, nan, ..., -2.66567 ,
-2.22315384, -0.99844656],
[ nan, nan, nan, ..., -2.44023533,
-1.45101395, -0.70535918],
[ nan, nan, nan, ..., -2.57074761,
-2.17096201, -1.2426105 ],
...,
[ nan, nan, nan, ..., nan,
nan, -0.01555032],
[ nan, nan, nan, ..., nan,
nan, 0.67034064],
[ nan, nan, nan, ..., nan,
nan, nan]],
[[ nan, nan, nan, ..., -2.49613084,
-2.17967113, -1.73102186],
[-0.19085314, nan, nan, ..., -2.24216667,
-1.60634045, -0.68446314],
[-0.47985682, -2.6897498 , nan, ..., -3.09012476,
-2.5473986 , -2.29142956],
...,
[ nan, nan, nan, ..., nan,
nan, -0.53484962],
[ nan, nan, nan, ..., nan,
nan, -0.24008878],
[ nan, nan, nan, ..., nan,
nan, -0.53568459]],
[[-0.93326312, 0.80451525, -2.70141882, ..., -1.50553473,
-1.36389759, -0.97748875],
[-0.92073983, 1.20430801, -0.23245133, ..., -2.48700702,
-1.71913439, -1.9057203 ],
[-0.28947553, 1.65032955, 0.21173892, ..., -2.3592551 ,
-2.42616244, -2.00186308],
...,
[ nan, nan, nan, ..., nan,
nan, -0.29636747],
[ nan, nan, nan, ..., nan,
nan, -0.21941662],
[ nan, nan, nan, ..., nan,
nan, 0.00773721]],
...,
[[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
...,
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]],
[[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
...,
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]],
[[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
...,
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]]])}
This file has a lot of information, but we are specifically interested in the location (X,Y), and velocity components (U,V). In this example we load a single time (first one).
[4]:
# for south Oahu, use this
# Note that X and Y are 2D fields for the location;
# U and V are functions of time (3D), so first we
# pick a specific time
# t = 0 for first time
# t = -1 for last time
t = 0
time = hfr['T']
X = hfr['X']
Y = hfr['Y']
uvel = hfr['U'][:,:,t]
vvel = hfr['V'][:,:,t]
And now a quick vector plot using quiver.
[5]:
fig = plt.subplots(figsize=(12,12))
plt.quiver(X,Y,uvel,vvel);
We can use the same map drawing commands as before to place the land and orient the plot.
[6]:
plt.figure(figsize=(10,10))
ax = plt.axes(projection = ccrs.PlateCarree())
#ax.add_feature(cf.COASTLINE)
#ax.add_feature(cf.LAND)
land_50m = cf.NaturalEarthFeature('physical', 'land', '50m',
edgecolor='face',
facecolor=cf.COLORS['land'])
ax.add_feature(land_50m, edgecolor='gray')
ax.set_extent([-158.25, -157.5, 20.7, 21.5], crs=ccrs.PlateCarree())
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=2, color='gray', alpha = 0.25, linestyle='--')
# labels on bottom and left axes
gl.xlabels_top = False
gl.ylabels_right = False
# define the label style
gl.xlabel_style = {'size': 10, 'color': 'black'}
gl.ylabel_style = {'size': 10, 'color': 'black'}
# now we define exactly which ones to label and spruce up the labels
gl.xlocator = mticker.FixedLocator([-158.2,-158.1,-158.0,-157.9,-157.8,-157.7,-157.6])
gl.ylocator = mticker.FixedLocator([20.8,20.9,21.0,21.1,21.2,21.3,21.4])
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
plt.quiver(X,Y,uvel,vvel,
transform=ccrs.PlateCarree());
/opt/tljh/user/lib/python3.9/site-packages/cartopy/mpl/gridliner.py:451: UserWarning: The .xlabels_top attribute is deprecated. Please use .top_labels to toggle visibility instead.
warnings.warn('The .xlabels_top attribute is deprecated. Please '
/opt/tljh/user/lib/python3.9/site-packages/cartopy/mpl/gridliner.py:487: UserWarning: The .ylabels_right attribute is deprecated. Please use .right_labels to toggle visibility instead.
warnings.warn('The .ylabels_right attribute is deprecated. Please '
[7]:
# redo the plot, but instead of a single time, do over the median values over
# all time
# NOTE: time is the third dimension, so we have to include "axis=2" to compute
# the median over time; without this, nanmedian would return a single number
# nanmedian over all x,y,y
uvel = np.nanmedian(hfr['U'],axis=2)
vvel = np.nanmedian(hfr['V'],axis=2)
/opt/tljh/user/lib/python3.9/site-packages/numpy/lib/nanfunctions.py:1096: RuntimeWarning: All-NaN slice encountered
result = np.apply_along_axis(_nanmedian1d, axis, a, overwrite_input)
[8]:
plt.figure(figsize=(10,10))
ax = plt.axes(projection = ccrs.PlateCarree())
#ax.add_feature(cf.COASTLINE)
#ax.add_feature(cf.LAND)
land_50m = cf.NaturalEarthFeature('physical', 'land', '50m',
edgecolor='face',
facecolor=cf.COLORS['land'])
ax.add_feature(land_50m, edgecolor='gray')
ax.set_extent([-158.25, -157.5, 20.7, 21.5], crs=ccrs.PlateCarree())
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=2, color='gray', alpha = 0.25, linestyle='--')
# labels on bottom and left axes
gl.xlabels_top = False
gl.ylabels_right = False
# define the label style
gl.xlabel_style = {'size': 10, 'color': 'black'}
gl.ylabel_style = {'size': 10, 'color': 'black'}
# now we define exactly which ones to label and spruce up the labels
gl.xlocator = mticker.FixedLocator([-158.2,-158.1,-158.0,-157.9,-157.8,-157.7,-157.6])
gl.ylocator = mticker.FixedLocator([20.8,20.9,21.0,21.1,21.2,21.3,21.4])
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
plt.quiver(X,Y,uvel,vvel,
transform=ccrs.PlateCarree());
A note about time¶
Typical of a matlab file, time is recorded as “days since 0-0-0” or alternately hours/minutes/seconds can be included. So, if we import the time and have a look, we see the following.
[9]:
time = hfr['T']
print("The first time is", time[0], "and the last time is", time[-1])
The first time is [734504.02083291] and the last time is [734868.97916488]
Matlab has a useful function to conver this to a calendar date, but to do this in python, we need to parse the dates out.
Note this only works on pandas data series and not regular numpy arrays. Recall you can use type(time) to see what we’re dealing with here.
[10]:
type(time)
[10]:
numpy.ndarray
[11]:
# convert to a flat pandas data series
time = pd.Series(hfr['T'].squeeze())
[12]:
# and now covert to calendar date
# You don't need worry about the details here, but hopefully the
# syntax is somewhat recognizable
calendar_time = time.apply(lambda \
Time: datetime.fromordinal(int(Time)) + timedelta(days=Time%1) \
- timedelta(days = 366 ))
[13]:
calendar_time
[13]:
0 2011-01-01 00:29:59.963069
1 2011-01-01 01:29:59.963056
2 2011-01-01 02:29:59.963043
3 2011-01-01 03:29:59.963029
4 2011-01-01 04:29:59.963016
...
8755 2011-12-31 19:29:59.845656
8756 2011-12-31 20:29:59.845642
8757 2011-12-31 21:29:59.845629
8758 2011-12-31 22:29:59.845615
8759 2011-12-31 23:29:59.845602
Length: 8760, dtype: datetime64[ns]
Second, we will load the matlab binary file “hfr_2003.mat”, again using the scipy.io routine.
[14]:
# secondary data set, west Oahu
hfr = sio.loadmat('/home/ocn463-data/hfr/hfr_2003.mat')
It contains:
Name, Size, Bytes, Class
CAUTION! configuration of arrays U,V for west Oahu are different from the south Oahu data set. Note that X and Y are 1-D fields covering both latitude and longitude. U and V are 2-D functions of time at each location.
[15]:
# first we pick a specific time
# t = 0 for first time
# t = -1 for last time
t = 1
time = hfr['T']
X = hfr['X']
Y = hfr['Y']
# for west Oahu, use this for 2D array
uvel = hfr['U'][t,:]
vvel = hfr['V'][t,:]
And again a quick vector plot using quiver:
[16]:
fig = plt.subplots(figsize=(12,12))
plt.quiver(X,Y,uvel,vvel);
We can again use the map drawing commands to place the land and orient the plot.
[17]:
plt.figure(figsize=(10,10))
ax = plt.axes(projection = ccrs.PlateCarree())
#ax.add_feature(cf.COASTLINE)
#ax.add_feature(cf.LAND)
land_50m = cf.NaturalEarthFeature('physical', 'land', '50m',
edgecolor='face',
facecolor=cf.COLORS['land'])
ax.add_feature(land_50m, edgecolor='gray')
ax.set_extent([-159.1, -158.1, 20.7, 21.7], crs=ccrs.PlateCarree())
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True,
linewidth=2, color='gray', alpha = 0.25, linestyle='--')
# labels on bottom and left axes
gl.xlabels_top = False
gl.ylabels_right = False
# define the label style
gl.xlabel_style = {'size': 10, 'color': 'black'}
gl.ylabel_style = {'size': 10, 'color': 'black'}
# now we define exactly which ones to label and spruce up the labels
gl.xlocator = mticker.FixedLocator([-159.0, -158.9, -158.8,-158.7,-158.6,-158.5,-158.4,-158.3,-158.2])
gl.ylocator = mticker.FixedLocator([20.7,20.8,20.9,21.0,21.1,21.2,21.3,21.4,21.5, 21.6])
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
plt.quiver(X,Y,uvel,vvel,
transform=ccrs.PlateCarree());
STOP HERE; OLD MATLAB CODE BELOW
[year,month,day]=datevec(T);
figure() subplot(221) fi = find(month==2); s=0.04; quiver(X,Y,nanmedian(U(:,:,fi),3)s,nanmedian(V(:,:,fi),3)s,0) axis square grid fi = find(month==5); subplot(222) quiver(X,Y,nanmedian(U(:,:,fi),3)s,nanmedian(V(:,:,fi),3)s,0) axis square grid fi = find(month==7); subplot(223) grid quiver(X,Y,nanmedian(U(:,:,fi),3)s,nanmedian(V(:,:,fi),3)s,0) axis square grid fi = find(month==11); subplot(224) quiver(X,Y,nanmedian(U(:,:,fi),3)s,nanmedian(V(:,:,fi),3)s,0) axis square grid
Let’s make a color map of the two-year mean currents
figure() contourf(X, Y, nanmedian(U,3)) axis square axis([-158.4 -157.5 20.65 21.3]) hold on title(‘Two year median of surface currents from HFR’) grid
You can also plot the data from different dimensions
figure() pcolor(squeeze(V(:, 20, :))) shading flat colorbar caxis([-6 6])
figure() feather(squeeze(U(10,20,:)), squeeze(V(10,20,:)))
figure() pcolor(T, Y(:,20), squeeze(U(:,20,:))) shading flat datetick(‘x’,’ mmyy’)
figure() pcolor(Y(:,20), T, squeeze(U(:,20,:))‘) shading flat datetick(’y’,’ mmyy’) colorbar() caxis([-6 6])
figure() pcolor(X(10,:), T, squeeze(U(10,:,:))‘) shading flat datetick(’y’,’ mmyy’) colorbar() caxis([-6 6])%Let’s make a zoom where we see more energy and find the propagation speed of the perturbation observed.
axis ([-158.17 -157.66 734721.13 734750.09])
figure() pcolor(X(10,:), T, squeeze(V(10,:,:))‘) shading flat datetick(’y’,’ mmyy’) colorbar() caxis([-6 6]) %Let’s make a zoom where we see more energy and find the propagation speed of the perturbation observed.
axis ([-158.17 -157.66 734721.13 734750.09])
Select two points in the graph to find the two latitudes and the two times
ginput(2) speed = diff(ans)
The first value will give you the degrees of longitude the perturbation moved and the second value will give you the time it took to move. You should get something about 60 km/day or 0.8 m/s.
vprog = speed(1) * cosd(21) * 10e7/90/(speed(2) *86400) % in m/s
Histograms and normal distributions
Selecting one point in the HFR spatial domain
up = squeeze(U(15, 15, :)); figure() subplot(211) hist(up, 32) % you can select the number of “little bars” to separate your distribution title(‘u’)
vp = squeeze(V(15, 15, :)); subplot(212) hist(vp, 32) % you can select the number of “little bars” to separate your distribution title(‘v’)
Selecting all the data in the HFR spatial domain
figure() u0 = U(:); hist(u0, 512)
figure() v0 = V(:); hist(v0, 512)