add CAPPI and CR method

kaomiao · Mar 1, 2020 · 7d9f2be · 7d9f2be
1 parent 3f09bab
commit 7d9f2be
Show file tree

Hide file tree

Showing 5 changed files with 182 additions and 54 deletions.
diff --git a/pycwr/core/NRadar.py b/pycwr/core/NRadar.py
@@ -5,10 +5,12 @@
 """
 import numpy as np
 import xarray as xr
+import pyproj
 from ..configure.default_config import DEFAULT_METADATA, CINRAD_field_mapping
 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
+from .RadarGridC import get_CR_xy, get_CAPPI_xy
 
 class PRD(object):
  """
@@ -126,6 +128,8 @@ def __init__(self, fields, scan_type, time, range, azimuth, elevation,latitude,
  self.nsweeps = nsweeps
  self.nrays = nrays
  self.sitename = sitename
+ self.get_vol_data()
+ self.product = xr.Dataset()
 
  def ordered_az(self, inplace=False):
  """
@@ -143,6 +147,159 @@ def ordered_az(self, inplace=False):
  prd_dat.fields.append(self.fields[isweep].swap_dims({"time": "azimuth"}).sortby("azimuth"))
  return prd_dat
 
+ def add_product_CR_xy(self, XRange, YRange):
+ """
+ 计算给定范围的组合反射率
+ :param XRange: np.ndarray, 1d, units:meters
+ :param YRange: np.ndarray, 1d, units:meters
+ :return:
+ """
+ GridX, GridY = np.meshgrid(XRange, YRange, indexing="ij")
+ vol_azimuth, vol_range, fix_elevation, vol_value, radar_height,\
+ radar_lon_0, radar_lat_0 = self.vol
+ fillvalue = -999.
+ GridV = get_CR_xy(vol_azimuth, vol_range, fix_elevation, vol_value,\
+ radar_height, GridX.astype(np.float64), GridY.astype(np.float64), -999.)
+ self.product.coords["x_cr"] = XRange
+ self.product.coords["y_cr"] = YRange
+ self.product["CR"] = (('x_cr', 'y_cr'), np.where(GridV==fillvalue, np.nan, GridV))
+ self.product.coords["x_cr"].attrs = {'units': 'meters',
+ 'standard_name': 'CR_product_x_axis ',
+ 'long_name': 'east_distance_from_radar',
+ 'axis': 'xy_coordinate',
+ 'comment': 'Distance from radar in east'}
+ self.product.coords["y_cr"].attrs = {'units': 'meters',
+ 'standard_name': 'CR_product_y_axis ',
+ 'long_name': 'north_distance_from_radar',
+ 'axis': 'xy_coordinate',
+ 'comment': 'Distance from radar in north'}
+ self.product["CR"].attrs = {'units': 'dBZ',
+ 'standard_name': 'Composite_reflectivity_factor',
+ 'long_name': 'Composite_reflectivity_factor',
+ 'axis': 'xy_coordinate',
+ 'comment': 'Maximum reflectance of all level',}
+
+ def add_product_CAPPI_xy(self, XRange, YRange, level_height):
+ """
+ 计算给定范围的CAPPI的图像
+ :param XRange: np.ndarray, 1d, units:meters
+ :param YRange: np.ndarray, 1d, units:meters
+ :param level_height: 要插值的高度，常量, units:meters
+ :return:
+ """
+ GridX, GridY = np.meshgrid(XRange, YRange, indexing="ij")
+ vol_azimuth, vol_range, fix_elevation, vol_value, radar_height, \
+ radar_lon_0, radar_lat_0 = self.vol
+ fillvalue = -999.
+ GridV = get_CAPPI_xy(vol_azimuth, vol_range, fix_elevation, vol_value, radar_height,
+ GridX.astype(np.float64), GridY.astype(np.float64), level_height, fillvalue)
+ self.product.coords["x_cappi_%d"%level_height] = XRange
+ self.product.coords["y_cappi_%d"%level_height] = YRange
+ self.product["CAPPI_%d"%level_height] = (("x_cappi_%d"%level_height, "y_cappi_%d"%level_height),
+ np.where(GridV==fillvalue, np.nan, GridV))
+ self.product.coords["x_cappi_%d"%level_height].attrs = {'units': 'meters',
+ 'standard_name': 'CAPPI_product_x_axis ',
+ 'long_name': 'east_distance_from_radar',
+ 'axis': 'xy_coordinate',
+ 'comment': 'Distance from radar in east'}
+ self.product.coords["y_cappi_%d"%level_height].attrs = {'units': 'meters',
+ 'standard_name': 'CAPPI_product_y_axis ',
+ 'long_name': 'north_distance_from_radar',
+ 'axis': 'xy_coordinate',
+ 'comment': 'Distance from radar in north'}
+ self.product["CAPPI_%d"%level_height].attrs = {'units': 'dBZ',
+ 'standard_name': 'Constant_altitude_plan_position_indicator',
+ 'long_name': 'Constant_altitude_plan_position_indicator',
+ 'axis': 'xy_coordinate',
+ 'comment': 'CAPPI of level %d m.'%level_height, }
+
+ def add_product_CR_lonlat(self, XLon, YLat):
+ """
+ 计算给定经纬度范围的组合反射率
+ :param XLon:np.ndarray, 1d, units:degree
+ :param YLat:np.ndarray, 1d, units:degree
+ :return:
+ """
+ fillvalue = -999.
+ GridLon, GridLat = np.meshgrid(XLon, YLat, indexing="ij")
+ projparams = {"proj": "aeqd", "lon_0": self.scan_info["longitude"].values,
+ "lat_0": self.scan_info["latitude"].values}
+ proj = pyproj.Proj(projparams)
+ GridX, GridY = proj(GridLon, GridLat, inverse=False)
+ vol_azimuth, vol_range, fix_elevation, vol_value, radar_height, \
+ radar_lon_0, radar_lat_0 = self.vol
+ GridV = get_CR_xy(vol_azimuth, vol_range, fix_elevation, vol_value, \
+ radar_height, GridX.astype(np.float64), GridY.astype(np.float64), -999.)
+ self.product.coords["lon_cr"] = XLon
+ self.product.coords["lat_cr"] = YLat
+ self.product["CR_geo"] = (('lon_cr', 'lat_cr'), np.where(GridV == fillvalue, np.nan, GridV))
+ self.product.coords["lon_cr"].attrs = {'units': 'degrees',
+ 'standard_name': 'CR_product_lon_axis ',
+ 'long_name': 'longitude_cr',
+ 'axis': 'lonlat_coordinate'}
+ self.product.coords["lat_cr"].attrs = {'units': 'degrees',
+ 'standard_name': 'CR_product_lat_axis ',
+ 'long_name': 'latitude_cr',
+ 'axis': 'lonlat_coordinate'}
+ self.product["CR_geo"].attrs = {'units': 'dBZ',
+ 'standard_name': 'Composite_reflectivity_factor_lonlat_grid',
+ 'long_name': 'Composite_reflectivity_factor',
+ 'axis': 'lonlat_coordinate',
+ 'comment': 'Maximum reflectance of all level', }
+
+ def add_product_CAPPI_lonlat(self, XLon, YLat, level_height):
+ """
+ 计算给定经纬度范围的CAPPI
+ :param XLon:np.ndarray, 1d, units:degrees
+ :param YLat:np.ndarray, 1d, units:degrees
+ :param level_height:常量，要计算的高度
+ :return:
+ """
+ fillvalue = -999.
+ GridLon, GridLat = np.meshgrid(XLon, YLat, indexing="ij")
+ projparams = {"proj": "aeqd", "lon_0": self.scan_info["longitude"].values,
+ "lat_0": self.scan_info["latitude"].values}
+ proj = pyproj.Proj(projparams)
+ GridX, GridY = proj(GridLon, GridLat, inverse=False)
+ vol_azimuth, vol_range, fix_elevation, vol_value, radar_height, \
+ radar_lon_0, radar_lat_0 = self.vol
+ GridV = get_CAPPI_xy(vol_azimuth, vol_range, fix_elevation, vol_value, radar_height,
+ GridX.astype(np.float64), GridY.astype(np.float64), level_height, fillvalue)
+ self.product.coords["lon_cappi_%d" % level_height] = XLon
+ self.product.coords["lat_cappi_%d" % level_height] = YLat
+ self.product["CAPPI_geo_%d" % level_height] = (("lon_cappi_%d" % level_height, "lat_cappi_%d" % level_height),
+ np.where(GridV == fillvalue, np.nan, GridV))
+ self.product.coords["lon_cappi_%d" % level_height].attrs = {'units': 'degrees',
+ 'standard_name': 'CAPPI_product_lon_axis ',
+ 'long_name': 'longitude_CAPPI',
+ 'axis': 'lonlat_coordinate',}
+ self.product.coords["lat_cappi_%d" % level_height].attrs = {'units': 'degrees',
+ 'standard_name': 'CAPPI_product_lat_axis ',
+ 'long_name': 'latitude_CAPPI',
+ 'axis': 'lonlat_coordinate',}
+ self.product["CAPPI_geo_%d" % level_height].attrs = {'units': 'dBZ',
+ 'standard_name': 'Constant_altitude_plan_position_indicator',
+ 'long_name': 'Constant_altitude_plan_position_indicator',
+ 'axis': 'lonlat_coordinate',
+ 'comment': 'CAPPI of level %d m' % level_height, }
+
+ def get_vol_data(self, field_name="dBZ", fillvalue=-999.):
+ """
+ 获取用于插值的雷达体扫数据
+ :return:
+ """
+ order_vol = self.ordered_az()
+ order_vol.fields = [order_vol.fields[i] for i in order_vol.scan_info["fixed_angle"].argsort().values]
+ order_vol.scan_info["fixed_angle"] = order_vol.scan_info["fixed_angle"].sortby("sweep")
+ vol_azimuth = [ppi.azimuth.values for ppi in order_vol.fields]
+ vol_range = [ppi.range.values for ppi in order_vol.fields]
+ fix_elevation = order_vol.scan_info["fixed_angle"].values
+ vol_value = [np.where(np.isnan(ppi[field_name].values), fillvalue, ppi[field_name].values) for ppi in order_vol.fields]
+ radar_height = float(order_vol.scan_info["altitude"].values)
+ radar_lon_0 = float(order_vol.scan_info["longitude"].values)
+ radar_lat_0 = float(order_vol.scan_info["latitude"].values)
+ self.vol = vol_azimuth, vol_range, fix_elevation.astype(np.float64), vol_value, radar_height, radar_lon_0, radar_lat_0
+
  def get_RHI_data(self, az, field_name="dBZ"):
  """
  获取RHI剖面数据
@@ -209,6 +366,7 @@ class PRD_AZ:
  def __init__(self):
  self.scan_info = None
  self.fields = []
+ self.product = xr.Dataset()
 
 
 

diff --git a/pycwr/core/RadarGrid.py b/pycwr/core/RadarGrid.py
@@ -1,5 +1,4 @@
 import numpy as np
-import pyproj
 from .transforms import cartesian_to_antenna_cwr, cartesian_xyz_to_antenna
 
 def interp_ppi(az, r, az_0, az_1, r_0, r_1, mat_00, mat_01, mat_10, mat_11, fillvalue=-999.):
@@ -107,27 +106,6 @@ def get_CR_xy(vol_azimuth, vol_range, fix_elevation, vol_value, radar_height, Gr
  GridValue = np.nanmax(np.where(GridValue == fillvalue, np.nan, GridValue), axis=0)
  return np.where(np.isnan(GridValue), fillvalue, GridValue)
 
-def get_CR_lonlat(vol_azimuth, vol_range, fix_elevation, vol_value, radar_height,
- radar_lon_0, radar_lat_0, GridLon, GridLat, fillvalue=-999.):
- """
- 对应经纬度网格, 插值组合反射率因子
- :param vol_azimuth:存放多个仰角体扫方位角的列表, list, units:degree
- :param vol_range:存放多个仰角体扫距离的列表, list, units:meters
- :param fix_elevation:每个仰角体扫对应的仰角， np.ndarray， 1d
- :param vol_value:存放多个仰角体扫数据的列表, list
- :param radar_height:常量, 雷达距离海平面的高度， units:meters
- :param radar_lon_0:常量，雷达站点所在的经度, units:degree
- :param radar_lat_0:常量，雷达站点所在的纬度, units:degree
- :param GridLon:要插值的二维格点Lon, np.ndarray, 2d, units:meters
- :param GridLat:要插值的二维格点Lat, np.ndarray, 2d, units:meters
- :param fillvalue: 常量，缺测值
- :return:
- """
- projparams = {"proj":"aeqd", "lon_0":radar_lon_0, "lat_0":radar_lat_0}
- proj = pyproj.Proj(projparams)
- GridX, GridY = proj(GridLon, GridLat, inverse=False)
- return get_CR_xy(vol_azimuth, vol_range, fix_elevation, vol_value, radar_height, GridX, GridY, fillvalue)
-
 def get_CAPPI_xy(vol_azimuth, vol_range, fix_elevation, vol_value, radar_height,
  GridX, GridY, level_height, fillvalue=-999.):
  """
@@ -197,33 +175,3 @@ def get_CAPPI_xy(vol_azimuth, vol_range, fix_elevation, vol_value, radar_height,
  IER1 = interp_azimuth(Grid_range[ix, iy], vol_range[ie][range_1 - 1], vol_range[ie][range_1], ER10, ER11, fillvalue)
  GridValue[ix, iy] = interp_azimuth(Grid_el[ix, iy], fix_elevation[ie - 1], fix_elevation[ie], IER0, IER1, fillvalue)
  return GridValue
-
-def get_CAPPI_lonlat(vol_azimuth, vol_range, fix_elevation, vol_value, radar_height,
- radar_lon_0, radar_lat_0, GridLon, GridLat, level_height, fillvalue=-999.):
- """
- 对应经纬度网格，插值CAPPI图像
- :param vol_azimuth:存放多个仰角体扫方位角的列表, list, units:degree
- :param vol_range:存放多个仰角体扫距离的列表, list, units:meters
- :param fix_elevation:每个仰角体扫对应的仰角， np.ndarray， 1d
- :param vol_value:存放多个仰角体扫数据的列表, list
- :param radar_height:常量, 雷达距离海平面的高度， units:meters
- :param radar_lon_0:常量，雷达站点所在的经度, units:degree
- :param radar_lat_0:常量，雷达站点所在的纬度, units:degree
- :param GridLon:要插值的二维格点经度, np.ndarray, 2d, units:degree
- :param GridLat:要插值的二维格点纬度, np.ndarray, 2d, units:degree
- :param level_height:常量,待插值的高度， units:meters
- :param fillvalue:常量, 缺测值
- :return:
- """
- projparams = {"proj": "aeqd", "lon_0": radar_lon_0, "lat_0": radar_lat_0}
- proj = pyproj.Proj(projparams)
- GridX, GridY = proj(GridLon, GridLat, inverse=False)
- return get_CAPPI_xy(vol_azimuth, vol_range, fix_elevation, vol_value, radar_height,
- GridX, GridY, level_height, fillvalue)
-
-
-if __name__ == "__main__":
- GridLon, GridLat = np.mgrid[110:115:0.01, 29:35:0.01]
- projparams = {"proj": "aeqd", "lon_0": 112, "lat_0": 31}
- proj = pyproj.Proj(projparams)
- GridX, GridY = proj(GridLon, GridLat, inverse=False)
diff --git a/pycwr/core/RadarGridC.pyx b/pycwr/core/RadarGridC.pyx
@@ -96,6 +96,7 @@ def interp_ppi(double az, double r, double az_0, double az_1, double r_0, double
  interped = fillvalue
  return interped
 
+@cython.boundscheck(False) # Deactivate bounds checking
 def ppi_to_grid(cnp.ndarray[cnp.float64_t, ndim=1] azimuth, cnp.ndarray[cnp.float64_t, ndim=1] ranges,
  double elevation, cnp.ndarray[cnp.float64_t, ndim=2] mat_ppi, double radar_height,
  cnp.ndarray[cnp.float64_t, ndim=2] GridX, cnp.ndarray[cnp.float64_t, ndim=2] GridY,
@@ -145,6 +146,7 @@ def ppi_to_grid(cnp.ndarray[cnp.float64_t, ndim=1] azimuth, cnp.ndarray[cnp.floa
  GridValue[ix, iy] = interp_ppi(az, r, az_last, azimuth[iaz], ranges[ir-1], ranges[ir], mat_ppi[iaz-1, ir-1], mat_ppi[iaz-1, ir], mat_ppi[iaz, ir-1], mat_ppi[iaz, ir], fillvalue)
  return GridValue
 
+@cython.boundscheck(False) # Deactivate bounds checking
 def get_CAPPI_xy(vol_azimuth, vol_range, cnp.ndarray[cnp.float64_t, ndim=1] fix_elevation, vol_value,
  double radar_height, cnp.ndarray[cnp.float64_t, ndim=2] GridX,
  cnp.ndarray[cnp.float64_t, ndim=2] GridY, double level_height, double fillvalue):
@@ -229,6 +231,7 @@ def interp_azimuth(double az, double az_0, double az_1, double dat_0, double dat
  else:
  return dat_0
 
+@cython.boundscheck(False) # Deactivate bounds checking
 def get_CR_xy(vol_azimuth, vol_range, cnp.ndarray[cnp.float64_t, ndim=1] fix_elevation, vol_value, double radar_height,
  cnp.ndarray[cnp.float64_t, ndim=2] GridX, cnp.ndarray[cnp.float64_t, ndim=2] GridY, double fillvalue):
  """

diff --git a/pycwr/core/transforms.py b/pycwr/core/transforms.py
@@ -153,6 +153,21 @@ def antenna_to_cartesian_cwr(ranges, azimuths, elevations, h):
 
  return x, y, z
 
+def cartesian_xyz_to_antenna(x, y, z, h):
+ """
+ 根据采样点距离雷达的x,y的水平距离,以及高度z, 以及雷达高度h
+ x, units:meters
+ y, units:meters
+ z, units:meters
+ h, units:meters
+ return ranges, azimuth, elevation
+ """
+ R = 8494666.6666666661 #等效地球半径
+ ranges = ((R + h) ** 2 + (R + z) ** 2 - 2 * (R + h) * (R + z) * np.cos((x ** 2 + y ** 2) ** 0.5 / R)) ** 0.5
+ elevation = np.arccos((R + z) * np.sin((x ** 2 + y ** 2) ** 0.5 / R) / ranges) * 180. / np.pi
+ azimuth = _azimuth(x, y)
+ return azimuth, ranges, elevation
+
 def cartesian_to_antenna_cwr(x, y, elevation , h):
  """根据采样点距离雷达的x,y的水平距离,以及雷达仰角
  return x, y, z
@@ -165,7 +180,7 @@ def cartesian_to_antenna_cwr(x, y, elevation , h):
  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)
+ ranges = np.tan(s/R)*(R+h)/np.cosh(El)
  z = (R+h)/np.cos(El + s/R)*np.cos(El) - R ##计算高度
  az = _azimuth(x, y) ##计算方位角
  return az, ranges, z

diff --git a/setup.py b/setup.py
@@ -5,13 +5,15 @@
 
 from setuptools import find_packages, setup
 import sys, os
+from Cython.Build import cythonize
+import numpy
 
 parent_dir = os.path.dirname(os.path.abspath(__file__))
 sys.path.append(parent_dir)
 
 DISTNAME = "pycwr"
 AUTHOR = "pycwr developers"
-AUTHOR_EMAIL = "zhengyunuist@gmail.com"
+AUTHOR_EMAIL = "YuZheng1206@outlook.com"
 URL = "https://github.com/YvZheng/pycwr"
 LICENSE='MIT'
 PYTHON_REQUIRES = ">=3.6"
@@ -50,5 +52,7 @@
  classifiers=CLASSIFIERS,
  include_package_data = True,
  packages=find_packages(parent_dir),
+ ext_modules=cythonize('./pycwr/core/RadarGridC.pyx'),
+ include_dirs=[numpy.get_include()]
 )