For now I am just trying to plot areas where visibility is less than 1mile and give it a vis value of 2 and then where visibility is less than 10 miles give it a vis value of 5. This idea works perfectly for my wind gust (gst) but when I plot my vis's they plot incorrectly and I know my array (nvis) is correct because when I plot it I get the correct numbers but then when I try to plot my vis it shows the image showed below. I can tell something is wrong with either my np.shape or np.where because when I change the values like (nvis < 10 or nvis < 5) it plots the exact same thing. I tried switching in order my nvis < 10 and nvis < 1 but nothing changes the plot.
Does anyone know what I am doing wrong because it should work perfectly fine since I am following almost the exact same process as my gst.
############# Imports
from datetime import datetime, timedelta
import requests
import metpy.calc as mpcalc
from metpy.units import units
import urllib.request
import xarray as xr
import matplotlib
import matplotlib.pyplot as plt
from shapely import geometry
import numpy as np
###########Time Stamps
run_dt = datetime.utcnow()-timedelta(hours=2) # previous hour
fhr = '00'
fhr1 = '01'
run = run_dt.strftime('%H')
ymd = run_dt.strftime('%Y%m%d')
##########Wind Data Import
# Data download
url = 'https://nomads.ncep.noaa.gov/cgi-bin/filter_rap.pl?file=rap.t{run}z.awp130\
pgrbf{fhr}.grib2&lev_0-6000_m_above_ground=on&lev_surface=on&var_GUST=on&var_VIS=on&var_VUCSH=on&var_VVCSH=on&leftlon=0&rightlon=360&toplat=90&bottomlat=-90&dir=%2Frap.{ymd}'
URL = url.format(run=run, fhr=fhr, ymd=ymd)
r = requests.get(URL)
with open('/home/jupyter-kfa5169/FDI/wind.grib00', 'wb') as fh:
fh.write(r.content)
ds = xr.open_dataset('/home/jupyter-kfa5169/FDI/wind.grib00', engine='cfgrib')
lon, lat = ds.longitude.values-360., ds.latitude.values
gust = mpcalc.smooth_gaussian((ds.variables['gust'].values*units('m/s')).to('mph').magnitude, 4)
visab = mpcalc.smooth_gaussian((ds.variables['vis'].values*units('m')).to('mile').magnitude, 4)
#vshear = ds.variables['vvcsh']
#ushear = ds.variables['vucsh']
################## Makes The Levels for Visability
nvis = np.copy(visab)
vis = np.ones(nvis.shape)
vis = np.where(nvis < 5, 5, vis)
vis = np.where(nvis < 10, 2, vis)
################## Makes The levels Gusts
ngust = np.copy(gust)
gst = np.ones(ngust.shape)
gst = np.where(ngust > 40, 3, gst)
gst = np.where(ngust > 50, 4, gst)
gst = np.where(ngust > 60, 7, gst)
#################################
url = 'https://nomads.ncep.noaa.gov/cgi-bin/filter_rap.pl?file=rap.t{run}z.awp130pgrbf{fhr1}.grib2&lev_surface=on&var_ASNOW=on&var_FROZR=on&var_FRZR=on&leftlon=0&rightlon=360&toplat=90&bottomlat=-90&dir=%2Frap.{ymd}'
URL = url.format(run=run, fhr1=fhr1, ymd=ymd)
r = requests.get(URL)
with open('/home/jupyter-kfa5169/FDI/wind.grib01', 'wb') as fh:
fh.write(r.content)
ds1 = xr.open_dataset('/home/jupyter-kfa5169/FDI/wind.grib01', engine='cfgrib')
vis = mpcalc.smooth_gaussian((ds1.variables['asnow'].values*units('m')).to('inch').magnitude, 4)
frz = ds1.variables['frzr']
############################# PLOT ####################
fig = plt.figure(figsize=(12,10))
ax = fig.add_subplot(111)
cs = ax.contourf(lon, lat, vis, levels=[0,1,2,3,4,5,6,7,10,15]) ## EDIT CONTOUR LEVELS HERE
plt.colorbar
plt.show()
Image of my output:

Image of output when I replace the vis with nvis for contourf and keep the same contour levels:

In addition, I guess you'll want to fix like this.
I recommend you to use a Python debugger to monitor data processing and catch bugs like these.