使用Pycinrad读取标准格式雷达基数据生成透明png图片

mrsoberli

# 脚本功能：读取标准格式雷达基数据，生成透明png图片，并获取png图片的四个角的经纬度位置
# 可供WEBGIS调用显示，例如leaflet、高德地图等；
# 要求：透明图没有任何边框，确保经纬度位置完成正确。

自己乱写了个，但不知道经纬度位置到底正确不？图片有没有空白边框？请教大佬们帮忙看看。

1. # -*- encoding: utf8 -*-
   
2. 
   
3. # 脚本功能：读取标准格式雷达基数据，生成透明png图片，并获取png图片的四个角的经纬度位置
   
4. # 可供WEBGIS调用显示，例如leaflet、高德地图等；
   
5. # 要求：透明图没有任何边框，确保经纬度位置完成正确。
   
6. import os
   
7. import datetime
   
8. import sys
   
9. import time
   
10. import cinrad
    
11. import cartopy.crs as ccrs
    
12. import matplotlib.pyplot as plt
    
13. from cartopy.mpl.geoaxes import GeoAxes
    
14. from cinrad.visualize.utils import *
    
15. from cinrad.common import get_dtype, is_radial
    
16. 
    
17. def plotTransparent(data,savePngFile,**kwargs):
    
18.     fig = plt.figure(figsize=(10,10),dpi=160)
    
19.     plt.axis("off")
    
20.     plt.gca().xaxis.set_major_locator(plt.NullLocator())
    
21.     plt.gca().yaxis.set_major_locator(plt.NullLocator())
    
22.     plt.subplots_adjust(top = 1, bottom = 0, right = 1, left = 0, hspace = 0, wspace = 0)
    
23.     plt.margins(0,0)
    
24.     proj = ccrs.AzimuthalEquidistant(
    
25.                 central_longitude=data.site_longitude,
    
26.                 central_latitude=data.site_latitude,
    
27.             )
    
28.     ax = fig.add_axes([0,0,1,1], projection=proj)
    
29.     # ax.background_patch.set_visible(False)
    
30.     # ax.outline_patch.set_visible(False)
    
31.     ax.set_aspect("equal")
    
32.     lon = data["longitude"].values
    
33.     lat = data["latitude"].values
    
34.     extent= [lon.min(), lon.max(), lat.min(), lat.max()]
    
35.     ax.set_extent(extent, crs=ccrs.PlateCarree())
    
36.     geoax: GeoAxes = ax
    
37.     dtype = get_dtype(data)
    
38.     pnorm = norm_plot[dtype]
    
39.     pcmap = cmap_plot[dtype]
    
40.     var = data[dtype].values
    
41.     data_crs = ccrs.PlateCarree()
    
42.     geoax.pcolormesh(lon, lat, var, norm=pnorm, cmap=pcmap, transform=data_crs, **kwargs)
    
43.     autoscale(geoax)
    
44.     plt.savefig(savePngFile,transparent=True, pad_inches=0)
    
45.    
    
46. def autoscale(geoax):
    
47.     llon, ulon, llat, ulat =geoax.get_extent()
    
48.     lon_delta = ulon - llon
    
49.     lat_delta = ulat - llat
    
50.     if lon_delta == lat_delta:
    
51.         return
    
52.     if lon_delta > lat_delta:
    
53.         # The long axis is x-axis
    
54.         lat_center = (ulat + llat) / 2
    
55.         lat_extend = lon_delta / 2
    
56.         llat = lat_center - lat_extend
    
57.         ulat = lat_center + lat_extend
    
58.     elif lon_delta < lat_delta:
    
59.         # The long axis is y-axis
    
60.         lon_center = (ulon + llon) / 2
    
61.         lon_extend = lat_delta / 2
    
62.         llon = lon_center - lon_extend
    
63.         ulon = lon_center + lon_extend
    
64.     # print([llon, ulon, llat, ulat])
    
65.     geoax.set_extent([llon, ulon, llat, ulat], geoax.projection)
    
66. 
    
67. bz2file="Z_RADR_I_Z9999_20210911190300_O_DOR_CAD_CAP_FMT.bin.bz2"
    
68. savePngFile="Z_RADR_I_Z9999_20210911190300_O_DOR_CAD_CAP_FMT.png"
    
69. f = cinrad.io.StandardData(bz2file)
    
70. rl = list(f.iter_tilt(230,'REF'))
    
71. cr = cinrad.calc.quick_cr(rl)
    
72. 
    
73. #获取png图片四个角的经纬度---start
    
74. lon = cr["longitude"].values
    
75. lat = cr["latitude"].values
    
76. extent= [lon.min(), lon.max(), lat.min(), lat.max()]
    
77. ex= [ str(x) for x in extent ]
    
78. #获取png图片四个角的经纬度---end
    
79. 
    
80. plotTransparent(cr,savePngFile)
    
81. plt.close("all")

复制代码

豆豆豆

您好想请问楼主现在怎么获取雷达基数据呀

btsldzz

同问，找雷达基数据找的好头疼！

小黄鹿

厉害！！！！！！！

VV09

您好想请问楼主现在怎么获取雷达基数据呀？头大

mrsoberli

大哥们，基数据是保密数据，除了内部门能获取，其他肯定不行啊，只能和气象局开展合作科研项目。

咱能看一下问题吗？

miracle

雷达基数据是保密的，我们是没有雷达站。虽然是气象局的，用基数据都很难要到

回想之后

网上有开源项目有示例数据。知乎搜搜，不谢。

mrsoberli

回想之后 发表于 2022-2-17 16:21
网上有开源项目有示例数据。知乎搜搜，不谢。

大神明示一下。。。哪个库有示例啊？

回想之后

你这个用的是三方包处理的。按你说的，生成的图片是没有边界框的。你说的四个角的位置是需要另一份文件提供的。你直接拿到一个基数据你是没办法拿到四个图片位置的。除非你能拿到气象局雷达的位置文件。啊哈哈哈哈哈哈哈哈哈哈哈哈。