fix vcs problems

meta.yaml

@@ -1,6 +1,6 @@
 package:
  name: pycwr
- version: 0.2.7
+ version: 0.2.8
 
 build:
  skip: True # [py<35]

pycwr/GraphicalInterface/RadarInterface.py

@@ -27,19 +27,19 @@
 
 
 class LineBuilder:
- def __init__(self, fig, ax, radar_data, product):
+ def __init__(self, fig, ax, radar_data, product, map_bool):
  self.ax = ax
  self.xs = []
  self.ys = []
  self.fig = fig
+ self.map = map_bool
  self.cid = self.fig.canvas.mpl_connect('button_press_event', self)
  self.cursor = self.fig.canvas.mpl_connect('motion_notify_event', self.mouse_move)
  self.radar_dat = radar_data
  self.product = product
 
  def __call__(self, event):
  if len(self.xs) < 2:
- print("click", event.xdata, event.ydata)
  self.xs.append(event.xdata)
  self.ys.append(event.ydata)
  if len(self.xs) == 1:
@@ -52,8 +52,13 @@ def __call__(self, event):
  cv = FigureCanvas(Figure(figsize=(8, 6)))
  ax = cv.figure.add_axes([0.1, 0.3, 0.8, 0.6])
  cax = cv.figure.add_axes([0.1, 0.1, 0.8, 0.06])
- VerticalSection.GUI_section(cv.figure, ax, cax, self.radar_dat, [self.xs[0]*1000, self.ys[0]*1000],\
+ if not self.map:
+ VerticalSection.GUI_section(cv.figure, ax, cax, self.radar_dat, [self.xs[0]*1000, self.ys[0]*1000],\
  [self.xs[1]*1000, self.ys[1]*1000], field_name[self.product])
+ else:
+ VerticalSection.GUI_section_map(cv.figure, ax, cax, self.radar_dat,
+ [self.xs[0], self.ys[0]], \
+ [self.xs[1], self.ys[1]], field_name[self.product])
  cv.show()
  self.fig.canvas.draw()
  else:
@@ -182,7 +187,8 @@ def on_actionvertical_changed(self):
  """垂直剖面的绘制"""
  if self.actionvertical.isChecked():
  try:
- self.linebuilder = LineBuilder(self.fig, self.ax, self.radar_dat, self.find_var_in_groupBox())
+ self.linebuilder = LineBuilder(self.fig, self.ax, self.radar_dat, self.find_var_in_groupBox(),\
+ self.actionwithmap.isChecked())
  self.clickevent = True
  except AttributeError:
  pass
@@ -199,7 +205,6 @@ def on_actionwithmap_changed(self):
  """
  Slot documentation goes here.
  """
-
  pass
 
  @pyqtSlot()
@@ -401,6 +406,16 @@ def plot_graph_PPI(self, radar, level, product, map, continuously):
  self.ax.tick_params(axis="x", which="both", direction='in')
  self.MplWidget.canvas.draw()
 
+ if self.actionvertical.isChecked(): #尝试重新绑定
+ try:
+ self.fig.canvas.mpl_disconnect(self.linebuilder.cid)
+ self.fig.canvas.mpl_disconnect(self.linebuilder.cursor)
+ self.linebuilder = LineBuilder(self.fig, self.ax, self.radar_dat, self.find_var_in_groupBox(), \
+ self.actionwithmap.isChecked())
+ self.clickevent = True
+ except AttributeError:
+ pass
+
  @pyqtSlot()
  def on_pushButton_clicked(self):
  """

pycwr/configure/default_config.py

@@ -232,6 +232,13 @@
  'valid_max': 0,
  'valid_min': 1,
  'coordinates': 'elevation azimuth range'},
+ "beam_width":{
+ 'units': 'degrees',
+ 'standard_name': 'clutter_phase_alignment',
+ 'long_name': 'Antenna beam width polarization',
+ 'valid_max': 0,
+ 'valid_min': 5,
+ 'coordinates': 'sweep'},
 }
 
 # CINRAD files

pycwr/core/NRadar.py

@@ -6,9 +6,9 @@
 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, antenna_vectors_to_cartesian_cwr
-from scipy import spatial
-from ..interp.RadarInterp import radar_interp2d_var
+from ..core.transforms import cartesian_to_geographic_aeqd,\
+ antenna_vectors_to_cartesian_cwr, antenna_vectors_to_cartesian_rhi, cartesian_to_antenna_cwr,\
+ antenna_vectors_to_cartesian_vcs
 
 class PRD(object):
  """
@@ -122,7 +122,8 @@ def __init__(self, fields, scan_type, time, range, azimuth, elevation,latitude,
  "nyquist_velocity":(['sweep',], nyquist_velocity),
  "unambiguous_range":(['sweep',], unambiguous_range),
  "rays_per_sweep": (['sweep',], sweep_end_ray_index-sweep_start_ray_index+1),
- "fixed_angle": (["sweep",], fixed_angle)},
+ "fixed_angle": (["sweep",], fixed_angle),
+ "beam_width":(["sweep",], 360./(sweep_end_ray_index-sweep_start_ray_index+1))},
  coords={"sweep": np.arange(nsweeps, dtype=int)})
  self.scan_info['latitude'].attrs = DEFAULT_METADATA['latitude']
  self.scan_info['longitude'].attrs = DEFAULT_METADATA['longitude']
@@ -135,45 +136,95 @@ def __init__(self, fields, scan_type, time, range, azimuth, elevation,latitude,
  self.scan_info['fixed_angle'].attrs = DEFAULT_METADATA['fixed_angle']
  self.scan_info['start_time'].attrs = DEFAULT_METADATA['start_time']
  self.scan_info['end_time'].attrs = DEFAULT_METADATA['end_time']
+ self.scan_info['beam_width'].attrs = DEFAULT_METADATA['beam_width']
  self.nsweeps = nsweeps
  self.nrays = nrays
  self.sitename = sitename
 
- def get_vertical_section(self, start_point, end_point, field_name):
+ def ordered_az(self, inplace=False):
  """
- :param start_point: units:m
- :param end_point: units:m
+ regrid radar object by azimuth
  :return:
  """
+ if inplace:
+ for isweep in self.scan_info.sweep.values:
+ self.fields[isweep] = self.fields[isweep].swap_dims({"time": "azimuth"}).sortby("azimuth") ##对数据重新排序
+ return None
+ else:
+ prd_dat = PRD_AZ()
+ prd_dat.scan_info = self.scan_info
+ for isweep in self.scan_info.sweep.values:
+ prd_dat.fields.append(self.fields[isweep].swap_dims({"time": "azimuth"}).sortby("azimuth"))
+ return prd_dat
+
+ def get_RHI_data(self, az, field_name="dBZ"):
+ """
+ 获取RHI剖面数据
+ :param az:
+ :param field_name:
+ :return:
+ """
+ mesh_RHI = []
+ mesh_RANGE = []
+ mesh_Z = []
+ order_dat = self.ordered_az()
+ for isweep in self.scan_info.sweep.values:
+ if (isweep>0) and (self.scan_info.fixed_angle.values[isweep] < self.scan_info.fixed_angle.values[isweep-1]):
+ continue ##remove VCP26 Type data
+ isweep_data = order_dat.fields[isweep].sel(azimuth=az, method='nearest')[field_name]
+ x, y, z = antenna_vectors_to_cartesian_rhi(isweep_data.range, isweep_data.azimuth,\
+ isweep_data.elevation, self.scan_info.altitude.values)
+ mesh_xy = np.sqrt(x**2 + y**2)
+ mesh_RHI.append(isweep_data.values.reshape(1,-1))
+ mesh_RANGE.append(mesh_xy)
+ mesh_Z.append(z)
+ return mesh_RANGE, mesh_Z, mesh_RHI
+
+ def get_vcs_data(self, start_point, end_point, field_name):
+ """
+ :param start_point:
+ :param end_point:
+ :param field_name:
+ :return:
+ """
+ order_dat = self.ordered_az() ##求排序后的数据
  start_x, start_y = start_point
  end_x, end_y = end_point
  bins_res = (self.fields[0].range[1] - self.fields[0].range[0]).values
- start_end_dis = np.sqrt((start_x-end_x)**2 + (start_y-end_y)**2)
- npoints = int(start_end_dis/(bins_res/2.) + 1)
+ start_end_dis = np.sqrt((start_x - end_x) ** 2 + (start_y - end_y) ** 2)
+ npoints = int(start_end_dis/bins_res + 1)
  x_line = np.linspace(start_x, end_x, npoints)
  y_line = np.linspace(start_y, end_y, npoints)
- target = np.c_[x_line, y_line]
  xy_line_1d = np.linspace(0, start_end_dis, npoints)
- z_values = []
- field_values = []
- #先对剖线取最邻近点
- for ifield in self.fields:
- _x, _y, _z = antenna_vectors_to_cartesian(ifield.range.values, \
- ifield.azimuth.values, ifield.elevation.values)
- kdtree = spatial.cKDTree(np.c_[_x.ravel(), _y.ravel()])
- _distance, _idx = kdtree.query(target, k=1, n_jobs=-1)
- _z_value = np.where(_distance > bins_res*2, np.nan, _z.ravel()[_idx])
- _field_value = np.where(_distance > bins_res*2, np.nan, ifield[field_name].values.ravel()[_idx])
- z_values.append(_z_value)
- field_values.append(_field_value)
- z_values = np.asarray(z_values)
- field_values = np.asarray(field_values)
- xy_line_values = np.stack([xy_line_1d,]*len(self.fields), axis=0)
- mask_flag = (np.isnan(z_values.ravel()) | np.isnan(field_values.ravel()))
- z_values_ravel = z_values.ravel()[~mask_flag]
- xy_line_values_ravel = xy_line_values.ravel()[~mask_flag]
- field_values_ravel = field_values.ravel()[~mask_flag]
- mesh_xy, mesh_z = np.mgrid[0:start_end_dis:npoints*1j, 0:np.nanmax(z_values):bins_res/2.] #生成剖面的网格
- grid_field = radar_interp2d_var(np.c_[xy_line_values_ravel, z_values_ravel], field_values_ravel,\
- (mesh_xy, mesh_z), bandwidth=360/self.fields[0].azimuth.size)
- return mesh_xy, mesh_z, grid_field
+ xy = np.stack([xy_line_1d, xy_line_1d], axis=0)
+ mesh_Z = []
+ mesh_vcs = []
+ mesh_xy = []
+ # 先对剖线取最邻近点
+ for isweep, ifield in enumerate(order_dat.fields):
+ if (isweep>0) and (self.scan_info.fixed_angle.values[isweep] < self.scan_info.fixed_angle.values[isweep-1]):
+ continue ##remove VCP26 Type data
+ az, ranges, _ = cartesian_to_antenna_cwr(x_line, y_line, self.scan_info.fixed_angle.values[isweep],\
+ self.scan_info.altitude.values)
+ vcs_data = ifield[field_name].sel(azimuth=xr.DataArray(az, dims="vcs_r"),\
+ range=xr.DataArray(ranges, dims="vcs_r"), method="nearest") ##选取插值的点
+ _, _, z = antenna_vectors_to_cartesian_vcs(vcs_data.range, vcs_data.azimuth, vcs_data.elevation,\
+ order_dat.scan_info.altitude.values, \
+ self.scan_info.beam_width.values[isweep])
+ mesh_xy.append(xy)
+ mesh_Z.append(z)
+ mesh_vcs.append(vcs_data.values.reshape(1,-1))
+ return mesh_xy, mesh_Z, mesh_vcs
+
+class PRD_AZ:
+ """
+ data obj for radar data, AZ as dims!
+ """
+ def __init__(self):
+ self.scan_info = None
+ self.fields = []
+
+
+
+
+

pycwr/core/transforms.py

@@ -153,19 +153,10 @@ def antenna_to_cartesian_cwr(ranges, azimuths, elevations, h):
 
  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)
@@ -264,6 +255,42 @@ def antenna_vectors_to_cartesian(ranges, azimuths, elevations, edges=False):
  rg, eleg = np.meshgrid(ranges, elevations)
  return antenna_to_cartesian(rg, azg, eleg)
 
+def antenna_vectors_to_cartesian_rhi(ranges, azimuths, elevations, h, BeamWidth=1):
+ """
+ 考虑波束宽度,来构建RHI扫描的坐标系
+ :param ranges: 距离
+ :param azimuths: 方位角
+ :param elevations: 获取该仰角的信息
+ :param BeamWidth: 波束宽度
+ :param h: 雷达的高度
+ :return:
+ """
+ assert elevations.size == 1, "not rhi, check!"
+ assert azimuths.size == 1, "not rhi, check!"
+ elevs = np.array([elevations-BeamWidth/2., elevations+BeamWidth/2.])
+ rg, azg = np.meshgrid(ranges, azimuths)
+ rg, eleg = np.meshgrid(ranges, elevs)
+ return antenna_to_cartesian_cwr(rg, azg, eleg, h)
+
+def antenna_vectors_to_cartesian_vcs(ranges, azimuths, elevations, h, BeamWidth=1):
+ """
+ 考虑波束宽度,来构建任意剖面的坐标系
+ :param ranges: 距离
+ :param azimuths: 方位角
+ :param elevations: 获取该仰角的信息
+ :param BeamWidth: 波束宽度
+ :param h: 雷达的高度
+ :return:
+ """
+ assert elevations.ndim == 1, "check ,may input data is not right!"
+ assert azimuths.ndim == 1, "check ,may input data is not right!"
+ assert ranges.ndim == 1, "check ,may input data is not right!"
+
+ eleg = np.stack([elevations-BeamWidth/2., elevations+BeamWidth/2.], axis=0)
+ rg = np.stack([ranges, ranges], axis=0)
+ azg = np.stack([azimuths, azimuths], axis=0)
+ return antenna_to_cartesian_cwr(rg, azg, eleg, h)
+
 
 def _interpolate_range_edges(ranges):
  """ Interpolate the edges of the range gates from their centers. """

pycwr/data/default_opendir.json

@@ -1 +1 @@
-{"lastOpenDir": "E:/RadarBaseData/WSR98D/hangzhou"}
+{"lastOpenDir": "C:/Users/zy/Desktop"}