# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # CAPPI插值示例
# 要求pycwr版本>=0.2.15

# %matplotlib inline
from pycwr.io import read_auto
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from pycwr.draw.RadarPlot import plot_xy, Graph, plot_lonlat_map, GraphMap

# ## 读取数据部分
#
#

PRD = read_auto("/Users/zhengyu/OneDrive/Work/13_天气雷达库_pycwr/test_data/Z9040.20190905.175751.AR2.bz2")

PRD.fields[0]

# ## 根据距离雷达中心的x和y生成cappi数据

x1d = np.arange(-200000, 200001, 1000) ##x方向1km等间距， -200km～200km范围
y1d = np.arange(-200000, 200001, 1000) ##y方向1km等间距， -200km～200km范围
PRD.add_product_CAPPI_xy(XRange=x1d, YRange=y1d, level_height=1500) ##插值1500m高度的
#XRange: np.ndarray, 1d, units:meters
# YRange: np.ndarray, 1d, units:meters
# level_height: 要插值的高度，常量, units:meters

PRD.product ##可以查看1500m的cappi的产品

grid_x, grid_y = np.meshgrid(x1d, y1d, indexing="ij")
fig, ax = plt.subplots()
plot_xy(ax, grid_x, grid_y, PRD.product.CAPPI_1500) ##画图显示

# +
#显示 0.5度仰角进行对比

fig, ax = plt.subplots()
graph = Graph(PRD)
graph.plot_ppi(ax, 0, "dBZ", cmap="CN_ref") ## 0代表第一层, dBZ代表反射率产品，cmap
plt.show()
# -

# ## 根据经纬度信息生成cappi数据

lon1d = np.arange(117, 121.0001, 0.01) ##lon方向0.01等间距，117-121 范围
lat1d = np.arange(28, 32.0001, 0.01) ##lat方向0.01等间距， 28～32度范围
PRD.add_product_CAPPI_lonlat(XLon=lon1d, YLat=lat1d, level_height=1500) ##插值1500m高度的
# XLon:np.ndarray, 1d, units:degrees
# YLat:np.ndarray, 1d, units:degrees
# level_height:常量，要计算的高度 units:meters

PRD.product ##可查看lat lon坐标的cappi

grid_lon, grid_lat = np.meshgrid(lon1d, lat1d, indexing="ij")
ax = plt.axes(projection=ccrs.PlateCarree())
plot_lonlat_map(ax, grid_lon, grid_lat, PRD.product.CAPPI_geo_1500, transform=ccrs.PlateCarree())
plt.show()##画图显示

ax = plt.axes(projection=ccrs.PlateCarree())
graph = GraphMap(PRD, ccrs.PlateCarree()) ##叠加地图
graph.plot_ppi_map(ax, 0, "dBZ", cmap="CN_ref") ## 0代表第一层, dBZ代表反射率产品
ax.set_title("Using pycwr for ploting data with map", fontsize=16)
plt.show()



StarcoderdataJupyterScriptsDedupFiltered