I'm trying to use Level2File to plot data from February of 2002 and I'm getting a number of inconsistent errors.
Sometimes I get errors on the call to Level2File depending on which .Z file I'm attempting to plot. Sometimes it's TypeError: unsupported operand type(s) for /: 'float' and 'NoneType'. Sometimes it's OSError: [Errno 22] Invalid argument.
Sometimes it starts by throwing warnings
Unable to read volume header. Attempting to read messages.
followed by a ton of Unknown message: and a number, before ending with
Remaining buffered message segments for message type(s): 15 (1) and then erroring on the call to Level2File
Occasionally the error is on the the line elv_angle = raddata.sweeps[sweep][0][0].el_angle instead, in which case it is an IndexError: list index out of range.
I just upgraded to metpy 1.1.0 and am running Python 3.9 on Win10.
Note: Right away I noticed the extension for the desired Level II files from Feb 2002 is .Z (whereas I think recent data has an extension of .gz or .ar2v?). Is it possible that these files are just too old to be processed by Level2File?
Below is a simplified version of the code I'm running (though it's basically just copied & pasted from the example code on Github).
Thanks!
import matplotlib as mpl
mpl.use('Agg')
from metpy.io import Level2File
from metpy.plots import colortables
from matplotlib.cm import get_cmap
import numpy as np
import os
import sys
import matplotlib.pylab as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.io.shapereader import Reader
import datetime
import matplotlib.patheffects as PathEffects
from time import time
from datetime import datetime, timedelta
def get_lonlat(x, y, nex_lon, nex_lat):
from pyproj import Proj
p = Proj(proj='aeqd', ellps='sphere',
lon_0=nex_lon,
lat_0=nex_lat)
return p(x,y,inverse=True)
def main(argv):
sweep = 0
varname = 'REF'
raddir = "C:\\Users\\rconn\\Documents\\WORK\\Cases\\20020203_Lake_Effect\\"
filename = "KMQT20020203_200100.Z"
raddata = Level2File(raddir+filename)
elv_angle = raddata.sweeps[sweep][0][0].el_angle
# Pull data out of the file
radtime = raddata.dt
timestr = str(radtime)
timeint = int(timestr[0:4]+timestr[5:7]+timestr[8:10]+timestr[11:13]+timestr[14:16])
# First item in ray is header, which has azimuth angle
az = np.array([ray[0].az_angle for ray in raddata.sweeps[sweep]])
# 5th item is a dict mapping a var name (byte string) to a tuple
# Plotting options that change based on desired variable
if varname=='REF':
var_hdr = raddata.sweeps[sweep][0][4][b'REF'][0]
var_range = np.arange(var_hdr.num_gates) * var_hdr.gate_width + var_hdr.first_gate
var = np.array([ray[4][b'REF'][1] for ray in raddata.sweeps[sweep]])
# Get center lat, lon for plotting
lat_0 = raddata.sweeps[0][0][1].lat
lon_0 = raddata.sweeps[0][0][1].lon
lat_c = lat_0 # For plot extent
lon_c = lon_0
deltalat = 1.8
latsouth = lat_c - deltalat/2.
latnorth = lat_c + deltalat/2.
deltalon = ((5./4.)*deltalat)/np.cos(np.deg2rad(lat_0))
lonwest = lon_c - deltalon/2.
loneast = lon_c + deltalon/2.
# Turn into an array, then mask
data = np.ma.array(var)
data[np.isnan(data)] = np.ma.masked
# Convert az,range to x,y
xlocs = var_range * np.sin(np.deg2rad(az[:, np.newaxis])) * 1000
ylocs = var_range * np.cos(np.deg2rad(az[:, np.newaxis])) * 1000
lon, lat = get_lonlat(xlocs, ylocs, lon_0, lat_0)
# Plot the data
fig = plt.figure(5, figsize=(10,8))
ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
ax.outline_patch.set_linewidth(1.)
ax.outline_patch.set_zorder(6)
# Plot the data
p = ax.pcolormesh(lon, lat, data, cmap=get_cmap('viridis'), vmin=-20, vmax=80, zorder=3)
# Set custom extent for map plot
ax.set_extent([lonwest,loneast,latsouth,latnorth],ccrs.PlateCarree())
if __name__ == '__main__':
main(sys.argv[1:])
