Plot RAAS PPI using python

ppi_example.py

import struct
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import cwbplot.cwb_colorbar as cwbcolor

##### C O N F I G #####
fname_x = 'PPI_20220629074850_b_x_01_00.5'
fname_y = 'PPI_20220629074850_b_y_01_00.5'
fname_feild = 'PPI_20220629074850_bDZ_01_00.5'

radarlon = 120.9075
radarlat = 24.8191

plot_range = 300 #km
range_index = np.arange(50, 301, 50) #50,100,.....,300
azi_index = np.arange(0, 361, 30) #0,30,60,.....,360

title = fname_feild[4:18] + ' ' + fname_feild[-3:] + ' (DZ)'
#######################

def read_radar_ppi(fname):
    data = []
    with open(fname, 'rb') as f:
        ifactor = struct.unpack('i', f.read(4))[0]
        nx = struct.unpack('i', f.read(4))[0]
        ny = struct.unpack('i', f.read(4))[0]
        #print('ifactor=', ifactor, 'nx=', nx, 'ny=', ny)
        for i in range(nx*ny):
            b = f.read(4)
            data.append(struct.unpack('i', b)[0])
            
    data = np.array(data, dtype=np.float32).reshape((nx, ny))
    data = data / ifactor
    data[data < -999] = np.nan
    return data

def xy2latlon(centerlon, centerlat, xdis, ydis):
    A = 6371.22
    gridlat = (ydis / A / np.pi * 180) + centerlat
    gridlon = (xdis / A / np.pi * 180) / np.cos((gridlat + centerlat) * np.pi / 360.0) + centerlon
    return gridlon, gridlat
    
def radar_distance_ring(rrr, sss, color, width, center_x, center_y):
    for i in rrr:
        y = i * np.sin(np.arange(0, 2*np.pi, np.pi/100))
        x = i * np.cos(np.arange(0, 2*np.pi, np.pi/100))
        x, y = xy2latlon(center_x, center_y, x, y)
        map_x, map_y = map(x, y)
        plt.plot(map_x, map_y, '--', color=color, linewidth=width)
        
    for i in sss:
        x = np.array([0, rrr[-1]*np.cos(i/180*np.pi)])
        y = np.array([0, rrr[-1]*np.sin(i/180*np.pi)])
        x, y = xy2latlon(center_x, center_y, x, y)
        map_x, map_y = map(x, y)
        plt.plot(map_x, map_y, '--', color=color, linewidth=width)

    
#read data from binary files
xdis = read_radar_ppi(fname_x)
ydis = read_radar_ppi(fname_y)
feild = read_radar_ppi(fname_feild)

#convert distance to lat/lon
gridlon, gridlat = xy2latlon(radarlon, radarlat, xdis, ydis)
llcrnrlon, llcrnrlat = xy2latlon(radarlon, radarlat, -plot_range, -plot_range)
urcrnrlon, urcrnrlat = xy2latlon(radarlon, radarlat, plot_range, plot_range)

#draw map
map = Basemap(projection="mill", llcrnrlon=llcrnrlon, urcrnrlon=urcrnrlon, llcrnrlat=llcrnrlat, urcrnrlat=urcrnrlat, resolution="h")
map.drawcoastlines()
map.drawcountries()
map.drawrivers()
map.drawparallels(np.arange(-90, 90, 1), labels=[1,0,0,0])
map.drawmeridians(np.arange(-180, 180, 1), labels=[0,0,0,1])
map.fillcontinents(color='lightgrey')
map.drawmapboundary()

radar_distance_ring(range_index, azi_index, 'b', 1, radarlon, radarlat)

#plot feild
map_x, map_y = map(gridlon, gridlat)
cb = cwbcolor.radar()
plt.contourf(map_x, map_y, feild, vmin=0, levels=cb['levels'][2:-1], norm=cb['norm'], cmap=cb['cmap'])

#colorbar
cbstr = [str(i) for i in cb["levels"][2:-1]]
cbar = plt.colorbar(ticks=cb["levels"][2:-1])
cbar.ax.set_yticklabels(cbstr)

plt.title(title)
plt.show()