|
#!/usr/bin/python3
|
|
|
|
import os
|
|
import sys
|
|
import pathlib
|
|
import numpy as np
|
|
from pyproj import Proj
|
|
from netCDF4 import Dataset
|
|
import cartopy
|
|
import matplotlib
|
|
import cartopy.crs as ccrs
|
|
import matplotlib.pyplot as plt
|
|
from metpy.plots import USSTATES, USCOUNTIES
|
|
import cartopy.feature as cfeature
|
|
|
|
# what is the top path where all of your folders are?
|
|
sys.path.append('/mnt/c/Users/MichealSimpson/Documents/Functions')
|
|
|
|
# Read in the data
|
|
radarfile = 'test.netcdf'
|
|
data = netcdf4.Dataset(radarfile)
|
|
|
|
def extract_single_variable(data,variable_name):
|
|
var = data[variable_name][:].data
|
|
return var
|
|
|
|
# Extract the metadata from the radar
|
|
zh = extract_single_variable(data,'DBZH'); zh_dict = {};
|
|
rangee = extract_single_variable(data,'range')
|
|
az_all = extract_single_variable(data,'azimuth'); az_dict = {};
|
|
sweep_start = extract_single_variable(data,'sweep_start_ray_index')
|
|
sweep_end = extract_single_variable(data,'sweep_end_ray_index')
|
|
sweep_num = extract_single_variable(data,'sweep_number')
|
|
|
|
# Loop through the number of sweeps
|
|
for i in sweep_num:
|
|
zh_dict[i] = zh[sweep_start[i]:sweep_end[i]+1]
|
|
az_dict[i] = az_all[sweep_start[i]:sweep_end[i]+1]
|
|
|
|
# Now we will try to convert from polar to cartesian so we can plot easily
|
|
x = rangee * np.sin(np.deg2rad(az_dict[0]))[:,None]
|
|
y = rangee * np.cos(np.deg2rad(az_dict[0]))[:,None]
|
|
|
|
# Setup a projection
|
|
dataproj = Proj(f"+proj=eqc +lat_0={latitude} +lat_ts={latitude} +lon_0={longitude} +ellps=WGS84 +units=m")
|
|
lons,lats = dataproj(x,y,inverse=True)
|
|
|
|
# Plot
|
|
# Mask data
|
|
data = zh_dict[0]
|
|
data = np.ma.masked_where(data == -999,data,copy=True)
|
|
|
|
# Color stuff
|
|
colormap = [[220,255,255],[210,180,210],[160,130,200],[115,70,120],[210,210,200],[200,200,120],[120,120,120],[0,240,240],[1,160,240],[0,0,240],[0,255,0],[0,200,0],[0,140,0],[255,255,0],[230,190,0],[255,140,0],[255,0,0],[220,0,0],[190,0,0],[255,0,255],[150,85,200],[255,255,255]]
|
|
boundaries = [-30,-25,-20,-15,-10,-5,0,5,10,15,20,25,30,35,40,45,50,55,60,65,70,75]
|
|
cmap_data = matplotlib.colors.ListedColormap(colormap)
|
|
norm_cmap = matplotlib.colors.BoundaryNorm(boundaries,cmap_data.N,clip=True)
|
|
|
|
# Plot the background
|
|
fig,ax = plt.subplots(1,1, figsize=(10,9), subplot_kw=dict(projection=ccrs.PlateCarree()))
|
|
im = ax.pcolormesh(lons,lats,data,cmap=cmap_data,norm=norm_cmap)
|
|
ax.set_ylim([-81,-79])
|
|
ax.set_xlim([35,37])
|
|
ax.add_feature(USCOUNTIES.with_scale('5m'),edgecolor='black')
|
|
ax.add_feature(USSTATES.with_scale('5m'),edgecolor='black')
|
|
ax.set_title(varname + '\n' + str(ymd) + '.' + str(hms))
|
|
gl = ax.gridlines(draw_labels=True)
|
|
gl.top_labels = False
|
|
gl.right_labels = False
|
|
cax = fig.add_axes([ax.get_position().x1+0.01,ax.get_position().y0,0.02,ax.get_position().height])
|
|
cbar = plt.colorbar(im,cax=cax); cbar.set_ticks(boundaries); cbar.set_ticklabels(boundaries)
|
|
plt.scatter(longitude,latitude,c='r',s=30)
|
|
plt.savefig('test.png',dpi=100)
|
|
plt.close()
|
