|
#!/usr/bin/python3
|
|
|
|
import os
|
|
import sys
|
|
import pathlib
|
|
import numpy as np
|
|
from pyproj import Proj
|
|
from netCDF4 import Dataset
|
|
|
|
# Where will the output path be?
|
|
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 loop through each of the tilts to qc everything
|
|
for tiltnum in range(0,len(az_dict)):
|
|
az = az_dict[tiltnum]
|
|
zh_tilt = zh_dict[tiltnum]
|
|
|
|
# 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)
|
|
|
|
# Output a netcdf file
|
|
full_outfile = 'test.nc'
|
|
|
|
# Start with the netcdf filename and all the good stuff
|
|
ncfile = Dataset(full_outfile,mode='w',format='NETCDF4_CLASSIC')
|
|
|
|
z = zh_dict[0]
|
|
# Add the dimensions of our data, lat, and lons.
|
|
ncfile.createDimension('y',((x.shape[0])))
|
|
ncfile.createDimension('x',((x.shape[1])))
|
|
|
|
# Now we want to add the data
|
|
datavar = ncfile.createVariable('Reflectivity',np.float64,(('y','x')))
|
|
datavar[:,:] = z
|
|
datavar.coordinates = "longitude latitude"
|
|
|
|
# And then add the latitude
|
|
datalat = ncfile.createVariable('latitude',np.float64,(('y','x')))
|
|
datalat[:,:] = lats
|
|
datalat.standard_name = 'latitude'
|
|
datalat.units = 'degrees_east'
|
|
|
|
# And then add the longitude
|
|
datalon = ncfile.createVariable('longitude',np.float64,(('y','x')))
|
|
datalon[:,:] = lons
|
|
datalon.standard_name = 'longitude'
|
|
datalon.units = 'degrees_north'
|
|
|
|
# Close the netcdf file
|
|
ncfile.close()
|
