fix WSR98D BUG

pycwr/core/NRadar.py

@@ -6,7 +6,7 @@
 import numpy as np
 import xarray as xr
 from ..configure.default_config import DEFAULT_METADATA, CINRAD_field_mapping
-from ..core.transforms import antenna_vectors_to_cartesian, cartesian_to_geographic_aeqd
+from ..core.transforms import antenna_vectors_to_cartesian, cartesian_to_geographic_aeqd, antenna_vectors_to_cartesian_cwr
 from scipy import spatial
 from ..interp.RadarInterp import radar_interp2d_var
 
@@ -92,8 +92,8 @@ def __init__(self, fields, scan_type, time, range, azimuth, elevation,latitude,
  self.fields.append(isweep_data)
  else:
  for idx, (istart, iend) in enumerate(zip(sweep_start_ray_index, sweep_end_ray_index)):
- x, y, z = antenna_vectors_to_cartesian(range[:bins_per_sweep[idx]], azimuth[istart:iend+1],\
- elevation[istart:iend+1])
+ x, y, z = antenna_vectors_to_cartesian_cwr(range[:bins_per_sweep[idx]], azimuth[istart:iend+1],\
+ elevation[istart:iend+1], altitude)
  lon, lat = cartesian_to_geographic_aeqd(x, y, longitude, latitude)
  isweep_data = xr.Dataset(coords={'azimuth': (['time', ], azimuth[istart:iend+1]),
  'elevation': (['time',], elevation[istart:iend+1]),

pycwr/core/transforms.py

@@ -103,6 +103,123 @@ def antenna_to_cartesian(ranges, azimuths, elevations):
  y = s * np.cos(theta_a)
  return x, y, z
 
+def antenna_to_cartesian_cwr(ranges, azimuths, elevations, h):
+ """
+ Return Cartesian coordinates from antenna coordinates.
+ Parameters
+ ----------
+ ranges : array
+ Distances to the center of the radar gates (bins) in meters.
+ azimuths : array
+ Azimuth angle of the radar in degrees.
+ elevations : array
+ Elevation angle of the radar in degrees.
+ h : constant
+ Altitude of the instrument, above sea level, units:m.
+ Returns
+ -------
+ x, y, z : array
+ Cartesian coordinates in meters from the radar.
+ Notes
+ -----
+ The calculation for Cartesian coordinate is adapted from equations
+ 2.28(b) and 2.28(c) of Doviak and Zrnic [1]_ assuming a
+ standard atmosphere (4/3 Earth's radius model).
+ .. math::
+ z = \\sqrt{r^2+R^2+2*r*R*sin(\\theta_e)} - R
+ s = R * arcsin(\\frac{r*cos(\\theta_e)}{R+z})
+ x = s * sin(\\theta_a)
+ y = s * cos(\\theta_a)
+ Where r is the distance from the radar to the center of the gate,
+ :math:`\\theta_a` is the azimuth angle, :math:`\\theta_e` is the
+ elevation angle, s is the arc length, and R is the effective radius
+ of the earth, taken to be 4/3 the mean radius of earth (6371 km).
+ References
+ ----------
+ .. [1] Doviak and Zrnic, Doppler Radar and Weather Observations, Second
+ Edition, 1993, p. 21.
+ """
+ theta_e = np.deg2rad(elevations) # elevation angle in radians.
+ theta_a = np.deg2rad(azimuths) # azimuth angle in radians.
+ R = 6371.0 * 1000.0 * 4.0 / 3.0 # effective radius of earth in meters.
+ r = ranges * 1.0 # distances to gates in meters.
+
+ # z = (r ** 2 + R ** 2 + 2.0 * r * R * np.sin(theta_e)) ** 0.5 - R
+ z = ((r * np.cos(theta_e)) ** 2 + \
+ (R + h + r * np.sin(theta_e)) ** 2) ** 0.5 - R
+ s = R * np.arcsin(r * np.cos(theta_e) / (R + z)) # arc length in m.
+ x = s * np.sin(theta_a)
+ y = s * np.cos(theta_a)
+
+ return x, y, z
+
+def test_xyz(ranges, azimuths, elevations, h):
+ theta_e = np.deg2rad(elevations)
+ theta_a = np.deg2rad(azimuths)
+ r = ranges * 1.0
+ R = 6371.0 * 1000.0 * 4.0 / 3.0
+ z = ((r * np.cos(theta_e)) ** 2 + \
+ (R + h + r * np.sin(theta_e)) ** 2) ** 0.5 - R
+ return ranges*np.sin(theta_a), ranges * np.cos(theta_a), z
+
+def cartesian_to_antenna_cwr(x, y, elevation , h):
+ """根据采样点距离雷达的x,y的水平距离,以及雷达仰角
+ return x, y, z
+和高度,计算该点雷达的斜距
+ ..math::
+ s = sqrt(x^2 + y^2)
+ r = sin(s/R)*(R+h)/cos(elevation)
+ R为地球半径m,h为雷达高度m,elevation为仰角degree
+ """
+ R = 6371.0 * 1000.0 * 4.0 / 3.0 # effective radius of earth in meters.
+ s = np.sqrt(x**2+y**2)
+ El = np.deg2rad(elevation)
+ ranges = np.tan(s/R)*(R+h)/np.cos(El)
+ z = (R+h)/np.cos(El + s/R)*np.cos(El) - R ##计算高度
+ az = _azimuth(x, y) ##计算方位角
+ return az, ranges, z
+
+def _azimuth(x, y):
+ '''根据某一点距离雷达x方向，y方向的距离，计算方位角，单位：弧度'''
+ az = np.pi / 2 - np.angle(x + y * 1j)
+ return np.where(az >= 0, az, 2 * np.pi + az) * 180 / np.pi
+
+def antenna_vectors_to_cartesian_cwr(ranges, azimuths, elevations, h=0, edges=False):
+ """
+ Calculate Cartesian coordinate for gates from antenna coordinate vectors.
+ Calculates the Cartesian coordinates for the gate centers or edges for
+ all gates from antenna coordinate vectors assuming a standard atmosphere
+ (4/3 Earth's radius model). See :py:func:`pyart.util.antenna_to_cartesian`
+ for details.
+ Parameters
+ ----------
+ ranges : array, 1D.
+ Distances to the center of the radar gates (bins) in meters.
+ azimuths : array, 1D.
+ Azimuth angles of the rays in degrees.
+ elevations : array, 1D.
+ Elevation angles of the rays in degrees.
+ edges : bool, optional
+ True to calculate the coordinates of the gate edges by interpolating
+ between gates and extrapolating at the boundaries. False to
+ calculate the gate centers.
+ Returns
+ -------
+ x, y, z : array, 2D
+ Cartesian coordinates in meters from the center of the radar to the
+ gate centers or edges.
+ """
+ if edges:
+ if len(ranges) != 1:
+ ranges = _interpolate_range_edges(ranges)
+ if len(elevations) != 1:
+ elevations = _interpolate_elevation_edges(elevations)
+ if len(azimuths) != 1:
+ azimuths = _interpolate_azimuth_edges(azimuths)
+ rg, azg = np.meshgrid(ranges, azimuths)
+ rg, eleg = np.meshgrid(ranges, elevations)
+ return antenna_to_cartesian_cwr(rg, azg, eleg, h)
+
 
 def antenna_vectors_to_cartesian(ranges, azimuths, elevations, edges=False):
  """

pycwr/draw/SingleRadarPlot.py

@@ -209,9 +209,9 @@ def plot_ppi(fig, ax, cx, x, y, radar_data, max_range=None, title=None, normvar=
  cmaps = plt.get_cmap(cmap)
  levels = MaxNLocator(nbins=cmap_bins).tick_values(vmin, vmax)
  norm = BoundaryNorm(levels, ncolors=cmaps.N, clip=True)
- RadarGraph._SetGrids(ax, range_cycle / 1000.)
  gci = ax.pcolormesh(x / 1000., y / 1000., radar_data, cmap=cmaps, \
- zorder=10, norm=norm)
+ zorder=0, norm=norm)
+ RadarGraph._SetGrids(ax, range_cycle / 1000.)
  ax.set_aspect("equal")
  ax.set_xlim([min_x, max_x])
  ax.set_ylim([min_y, max_y])
@@ -247,12 +247,12 @@ def _SetGrids(ax, max_range, deltaR=50):
  for i in rings:
  x0 = i * np.cos(theta)
  y0 = i * np.sin(theta)
- ax.plot(x0, y0, linestyle='-', linewidth=0.6, color='#DCDCDC')
+ ax.plot(x0, y0, linestyle='-', linewidth=0.6, color='#5B5B5B')
 
  for rad in np.arange(0, np.pi, np.pi / 6):
  ax.plot([-1 * rings[-1] * np.sin(rad), rings[-1] * np.sin(rad)], \
  [-1 * rings[-1] * np.cos(rad), rings[-1] * np.cos(rad)], \
- linestyle='-', linewidth=0.6, color='#DCDCDC')
+ linestyle='-', linewidth=0.6, color='#5B5B5B')
 
  @staticmethod
  def _FixTicks(ticks):

pycwr/io/WSR98DFile.py

@@ -83,9 +83,9 @@ def _parse_radial_single(self):
  Data_buf = self.fid.read(Momheader['Length'])
  assert (Momheader['BinLength'] == 1) | (Momheader['BinLength'] == 2), "Bin Length has problem!"
  if Momheader['BinLength'] == 1:
- dat_tmp = (np.frombuffer(Data_buf, dtype="u1", offset=dtype_98D.MomentHeaderBlockSize)).astype(int)
+ dat_tmp = (np.frombuffer(Data_buf, dtype="u1", offset=0)).astype(int)
  else:
- dat_tmp = (np.frombuffer(Data_buf, dtype="u2", offset=dtype_98D.MomentHeaderBlockSize)).astype(int)
+ dat_tmp = (np.frombuffer(Data_buf, dtype="u2", offset=0)).astype(int)
  radial_var[dtype_98D.flag2Product[Momheader['DataType']]] = np.where(dat_tmp >= 5, \
  (dat_tmp - Momheader['Offset']) /
  Momheader['Scale'],