spurious horizontal contour lines depending on number of contour levels

[jabaars]

I'm seeing an issue where spurious horizontal contour lines appear on my map depending on the number of contour levels I choose. I can also make the horizontal lines go away by 'zooming out' or broadening my plot domain.

Here's an example that has more contour levels and shows the problem (note the two horizontal lines, one near middle of plot, one near bottom):

And here's an example that has less contour levels and the problem goes away:

[jabaars]

This is the main mapping part of my code:

   ur_lat = maxlat
    ur_lon = maxlon
    ll_lat = minlat
    ll_lon = minlon
    lat_ctr = ll_lat + ((ur_lat - ll_lat) * 0.5)
    lon_ctr = ll_lon + ((ur_lon - ll_lon) * 0.5)

    fig = plt.figure(figsize=(width,height))
    # left, bottom, width, height:
    ax = fig.add_axes([0.00,0.05,0.99,0.91])
   map = Basemap(resolution=res,projection='lcc',\
                  llcrnrlon= ll_lon, llcrnrlat=ll_lat,\
                  urcrnrlon= ur_lon, urcrnrlat= ur_lat,\
                  lat_0=lat_ctr,lon_0=lon_ctr,lat_1=(ur_lat - ll_lat))

    #--- Get lat and lon data in map's x/y coordinates.
    lon,lat  = np.meshgrid(lon,lat) 
    x,y = map(lon, lat)

#    test = np.where(lon >= minlon and lon <= maxlon and \
#                    lat >= minlat and lat <= maxlat)
#    print 'test = ', test

    #--- Draw coastlines, country boundaries, fill continents.
    map.drawcoastlines(linewidth = maplw)
    map.drawstates(linewidth = maplw)
    map.drawcountries(linewidth = maplw)
    map.fillcontinents(color='lightgray')

    #--- Draw the edge of the map projection region (the projection limb)
    map.drawmapboundary(linewidth = maplw)

    (irma_dt, irma_time, irma_lat, irma_lon, irma_wind, \
     irma_pres, irma_type) = read_irma_file(irma_file)

    #-- Add Irma track.
    for n in range(len(irma_dt)):
        xn, yn = map(irma_lon[n], irma_lat[n])
        plt.plot(xn, yn, 'ko', markersize = 5)
        id_tmp = irma_dt[n][0:6]
        idt = re.sub('/', '', irma_dt[n][0:6])
        dt_plot = idt + irma_time[n][0:3] + 'Z'        
        plt.text(xn, yn, dt_plot, fontsize=6, fontweight = 'bold', \
                 va = 'top', ha = 'right')
    xns, yns = map(irma_lon, irma_lat)        
    plt.plot(xns, yns, '-', color = 'black')

#    #--- Draw lat/lon grid lines every 30 degrees.
#    map.drawmeridians(np.arange(0, 360, 30), linewidth = maplw)
#    map.drawparallels(np.arange(-90, 90, 30), linewidth = maplw)

    cs = plt.contour(x, y, grid, levs_pres, nchunk = 1)

#    cbar = map.colorbar(cs, location='bottom', pad="3%", size=0.1, ticks=levs)
#    cbar.set_label(colorbar_lab, fontsize=fs, size=fs-1)
#    cbar.ax.tick_params(labelsize=fs-1)

    plt.title(title, fontsize=titlefs, fontweight='bold')

    #--- Save plot.
    print 'creating:\nxli ', plotfname, ' &'
    plt.savefig(plotfname)

    plt.close()

[QuLogic]

PS, use triple backticks for code blocks.

[jklymak]

Can you make a minimal self-contained example that triggers this error? W/o access to grid or lev_pres, its going to be pretty hard to know why this error happens. Also, when writing such an example 95% of the time shows what I did wrong ;-)

[jabaars]

Sorry for the newbie mistakes-- my first post here! Here's a stripped down version of the code:

#!/usr/bin/python
import  os, os.path, sys, string
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap, cm
import pickle

(lat, lon, grid) = pickle.load(open('github_example.pickle', 'rb'))

ur_lat = 32.0
ll_lat =  15.0
ur_lon = -70.0
ll_lon = -88.0 

levs_pres = np.arange(90000, 103000, 1000).tolist()
#levs_pres = np.arange(90000, 103000, 2000).tolist()

lat_ctr = ll_lat + ((ur_lat - ll_lat) * 0.5)
lon_ctr = ll_lon + ((ur_lon - ll_lon) * 0.5)

fig = plt.figure(figsize=(6.75, 6.75))
ax = fig.add_axes([0.00,0.05,0.99,0.91])

map = Basemap(resolution = 'i', projection='lcc',\
              llcrnrlon= ll_lon, llcrnrlat=ll_lat,\
              urcrnrlon= ur_lon, urcrnrlat= ur_lat,\
              lat_0=lat_ctr,lon_0=lon_ctr,lat_1=(ur_lat - ll_lat))

#--- Get lat and lon data in map's x/y coordinates.
lon,lat  = np.meshgrid(lon,lat) 
x,y = map(lon, lat)

#--- Draw coastlines, country boundaries, fill continents.
maplw = 0.6
map.drawcoastlines(linewidth = maplw)
map.drawstates(linewidth = maplw)
map.drawcountries(linewidth = maplw)
map.fillcontinents(color='lightgray')
map.drawmapboundary(linewidth = maplw)

cs = plt.contour(x, y, grid, levs_pres)

plt.savefig('test.png')
plt.close()

It does require a pickle file containing my grid which can be found here: https://atmos.washington.edu/~jbaars/github_issue/github_example.pickle

Setting levs_pres to increment every 1000 causes the horizontal lines. Every 2000 (commented out) does not.

[jklymak]

No doubt this is a bug, but contouring is very difficult, even on a Cartesian grid.

A workaround is to subset your data somewhat so you aren't using the whole globe's pressure field.

    lon[lon>180.] = lon[lon>180.] - 360.
    inx = (lon>ll_lon-2.) & (lon<=ur_lon+2.)
    iny = (lat>ll_lat-2.) & (lat<=ur_lat+2.)
    lat = lat[iny]
    lon = lon[inx]
    grid = grid[iny,:][:,inx]

[jabaars]

Hey thanks that did the trick! I had to turn lat and lon in np arrays but otherwise that fixed it. Much appreciated!

[jklymak]

No doubt because lcc has a big discontinuity at the dateline (or 180 deg from where you center it). You may have fixed your problem by using mercator as well.

[ThunderRoad75]

Is there a reason why this would be happening on a non-global grid as well?

I download RAP grb2 files from NCDC's NOMADS archive (https://nomads.ncdc.noaa.gov/data/rucanl/) and converted to netcdf via wgrib2.

The 'latitude' and 'longitude' arrays are shape (337, 451). The lats span 16.28 to 58.37, and the lons span 220.14 to 302.62. So clearly this is not a global grid. I've tried significant sub-setting of the grid (taking as much as 50 grid points off in each direction) but the problem persists. Even if I reduce my contour interval to something so high the plot is unusable, like 30 dam at 300 mb, the problem persists.

Attached is an output image of what these spurious horizontal contours look like. (Also note the background contourf cutting off on the edges of the domain because this is an example of the grid being subsetted by 50 grid points in each direction.)