Error plotting data with Orthographic() map projection

[joshdorrington]

trafficstars

I am trying to plot a netcdf file, which can be found here: testcube.zip

The issue is that the Orthographic projection incorrectly plots the data. Seemingly, all negative values are set to 0 during plotting. This is not the case for other projections. My code is below:

Code to reproduce

import matplotlib.pyplot as plt
import iris
import cartopy
import cartopy.crs as ccrs
import iris.plot as iplt
print(cartopy.__version__)

cube=iris.load_cube("testcube.nc")
print(cube)
print(cube.coord("latitude"))
print(cube.coord("longitude"))

fig,axes=plt.subplots(4)
clevs=np.linspace(-abs(cube.data).max(),abs(cube.data).max(),21)
projs=[ccrs.PlateCarree(),ccrs.Mollweide(),ccrs.Orthographic(0,90),ccrs.NearsidePerspective(central_latitude=50,central_longitude=-20)]
for i,proj in enumerate(projs):
    ax=plt.subplot(4,1,i+1,projection=proj)
    ax.coastlines()
    plot=iplt.contourf(cube, levels=clevs, cmap=plt.cm.RdBu_r,axes=ax)
    plt.colorbar(mappable=plot)

Which produces the following output:

'0.17.0'

unknown / (unknown)                 (latitude: 73; longitude: 144)
     Dimension coordinates:
          latitude                           x              -
          longitude                          -              x
     Scalar coordinates:
          mean cluster composites: 0
DimCoord(array([ 90. ,  87.5,  85. ,  82.5,  80. ,  77.5,  75. ,  72.5,  70. ,
        67.5,  65. ,  62.5,  60. ,  57.5,  55. ,  52.5,  50. ,  47.5,
        45. ,  42.5,  40. ,  37.5,  35. ,  32.5,  30. ,  27.5,  25. ,
        22.5,  20. ,  17.5,  15. ,  12.5,  10. ,   7.5,   5. ,   2.5,
         0. ,  -2.5,  -5. ,  -7.5, -10. , -12.5, -15. , -17.5, -20. ,
       -22.5, -25. , -27.5, -30. , -32.5, -35. , -37.5, -40. , -42.5,
       -45. , -47.5, -50. , -52.5, -55. , -57.5, -60. , -62.5, -65. ,
       -67.5, -70. , -72.5, -75. , -77.5, -80. , -82.5, -85. , -87.5,
       -90. ], dtype=float32), standard_name='latitude', units=Unit('degrees'), var_name='latitude')
DimCoord(array([  0. ,   2.5,   5. ,   7.5,  10. ,  12.5,  15. ,  17.5,  20. ,
        22.5,  25. ,  27.5,  30. ,  32.5,  35. ,  37.5,  40. ,  42.5,
        45. ,  47.5,  50. ,  52.5,  55. ,  57.5,  60. ,  62.5,  65. ,
        67.5,  70. ,  72.5,  75. ,  77.5,  80. ,  82.5,  85. ,  87.5,
        90. ,  92.5,  95. ,  97.5, 100. , 102.5, 105. , 107.5, 110. ,
       112.5, 115. , 117.5, 120. , 122.5, 125. , 127.5, 130. , 132.5,
       135. , 137.5, 140. , 142.5, 145. , 147.5, 150. , 152.5, 155. ,
       157.5, 160. , 162.5, 165. , 167.5, 170. , 172.5, 175. , 177.5,
       180. , 182.5, 185. , 187.5, 190. , 192.5, 195. , 197.5, 200. ,
       202.5, 205. , 207.5, 210. , 212.5, 215. , 217.5, 220. , 222.5,
       225. , 227.5, 230. , 232.5, 235. , 237.5, 240. , 242.5, 245. ,
       247.5, 250. , 252.5, 255. , 257.5, 260. , 262.5, 265. , 267.5,
       270. , 272.5, 275. , 277.5, 280. , 282.5, 285. , 287.5, 290. ,
       292.5, 295. , 297.5, 300. , 302.5, 305. , 307.5, 310. , 312.5,
       315. , 317.5, 320. , 322.5, 325. , 327.5, 330. , 332.5, 335. ,
       337.5, 340. , 342.5, 345. , 347.5, 350. , 352.5, 355. , 357.5],
      dtype=float32), standard_name='longitude', units=Unit('degrees'), var_name='longitude', circular=True)

[joshdorrington]

I have just noted that this is a problem with the contourf function, but does not occur when using pcolormesh:

projs=[ccrs.PlateCarree(),ccrs.Mollweide(),\
       ccrs.Orthographic(0,90),ccrs.NearsidePerspective(central_latitude=50,central_longitude=-20)]

for i,proj in enumerate(projs):
    ax=plt.subplot(4,1,i+1,projection=proj)
    ax.coastlines()
    plot=iplt.pcolormesh(cube, cmap=plt.cm.RdBu_r,axes=ax)
    plt.colorbar(mappable=plot)

[greglucas]

I'm pretty sure this is the same issue as https://github.com/SciTools/cartopy/issues/1076 where one of the contour levels is essentially covering up the rest of the underlying contours and not being clipped properly. There is a proposed fix linked in a PR at the bottom of that issue.

[rcomer]

Here is a reproducer without Iris:

import matplotlib.pyplot as plt
import cartopy
import cartopy.util
import cartopy.crs as ccrs
import numpy as np
from xarray import open_dataset
print(cartopy.__version__)


ds = open_dataset('testcube.nc')
da = ds.unknown

data, lons, lats = cartopy.util.add_cyclic(da.data, x=da.longitude.data, y=da.latitude.data)

clevs=np.linspace(-abs(da.data).max(),abs(da.data).max(),21)
projs=[ccrs.PlateCarree(),ccrs.Mollweide(),ccrs.Orthographic(0,90),ccrs.NearsidePerspective(central_latitude=50,central_longitude=-20)]
for i,proj in enumerate(projs):
    ax=plt.subplot(4,1,i+1,projection=proj)
    ax.coastlines()
    plot=ax.contourf(lons, lats, data, levels=clevs, cmap=plt.cm.RdBu_r, transform=ccrs.PlateCarree())
    plt.colorbar(mappable=plot)

plt.show()