Admin 1 features incompletely drawn in some cases

[mhearne-usgs]

Description

Admin level 1 Natural Earth features are incompletely drawn in some cases - the script below creates a map showing the border between Illinois and Indiana, but that ends before the bottom of the map where Illinois, Indiana, and Kentucky all meet. Attached is a map of the same data from Natural Earth plotted in QGIS, with the junction of the three states being shown.

Code to reproduce

#!/usr/bin/env python


import cartopy.crs as ccrs  # projections
import cartopy.feature as cfeature
import matplotlib.pyplot as plt

BORDER_ZORDER = 1110

if __name__ == '__main__':
    bounds = (-88.76666666666667, -86.93333333333334,
              37.73333333333333, 39.18333333333333)
    figsize = (10.0, 9.947416426343834)
    xmin, xmax, ymin, ymax = bounds
    clon = (xmin + xmax) / 2
    if xmax < 0 and xmax < xmin:
        clon = (xmin + (xmax+360))/2

    clat = (ymin + ymax) / 2
    proj = ccrs.Mercator(central_longitude=clon,
                         min_latitude=ymin,
                         max_latitude=ymax,
                         globe=None)
    geoproj = ccrs.PlateCarree()
    figure = plt.figure(figsize=figsize)

    aspect = 0.9941078696774519

    # Note: dimensions are: [left, bottom, width, height]
    dim_left = 0.1
    dim_bottom = 0.19
    dim_width = 0.8
    dim_height = dim_width/aspect

    dimensions = [dim_left, dim_bottom, dim_width, dim_height]

    ax = figure.add_axes(dimensions, projection=proj)
    ax.set_extent([xmin, xmax, ymin, ymax], crs=geoproj)
    states_provinces = cfeature.NaturalEarthFeature(
        category='cultural',
        name='admin_1_states_provinces_lines',
        scale='10m',
        facecolor='none')

    ax.add_feature(states_provinces, edgecolor='black', zorder=BORDER_ZORDER)
    plt.savefig('maptest.png')

Traceback

Full environment definition

Operating system

Cartopy version

conda list

pip list

[QuLogic]

This seems to have something to do with the thresholding; you can get it to work by patching it with proj._threshold /= 6.

The threshold is determined by the x limits, but should probably check the y limits as well. I don't really know how the threshold calculation was chosen though, so I'm not sure what the best calculation should be there (just doing the same thing would produce a threshold much much smaller than 1/6).

[QuLogic]

I think it also works fine if you take out min_latitude because then it's not clipping that close to the map border.