Problem with ortho projection and pcolormesh

[ricitron]

I have trouble with the ortho projection and pcolormesh.

It should plot a mesh of grid points. Instead, in the upper right portion of the sphere it plots strange lines instead of grid points. The mapping of the mesh looks off.

The below example is for a mesh of random noise, the mesh is clearly distorted / misplotted in the upper right portion of the ortho projection sphere.

Example:

from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt

plt.clf()

dpp =1 # degrees per pixel
lons = np.arange(-180,180+dpp,dpp)
lats = -1*np.arange(-90,90+dpp,dpp)

m = Basemap(projection='ortho', lon_0=0, lat_0=-60, resolution='l')
data = np.random.random((np.size(lats), np.size(lons)))
lons, lats = np.meshgrid(lons, lats)
x, y = m(lons, lats)

im = m.pcolormesh(x, y, data, latlon=False, cmap='RdBu')
#im = m.pcolormesh(lons, lats, data, latlon=True, cmap='RdBu')

m.colorbar(im)
plt.show()

[2sn]

i have the same issue, irrespective of first mapping coordinates or using latlon=True. Basemap.pcolor works (correct image):

but replaceing the same call by Basemap.colormesh does not work

[2sn]

After some experimenting, I find that issues comes from ax.pcolormesh not truncating triangles with data points of 1e30 - to which invalid points are mapped by the projection. Setting them to np.nan, which would have seem the easiest solution, does not work

from mpl_toolkits.basemap import Basemap
import numpy as np
m = Basemap(projection='ortho',lat_0=45,lon_0=30)
# define lon lat mesh and data array
x,y = m(lon,lat)
x[x>1e20]=np.nan
y[y>1e20]=np.nan
cs = m.pcolormesh(x, y, data)

the result is

ValueError: x and y arguments to pcolormesh cannot have non-finite values or be of type numpy.ma.core.MaskedArray with masked values

What one can do, however, is to set data values to np.nan where any of the edges is set to invalid value:

from mpl_toolkits.basemap import Basemap
import numpy as np
m = Basemap(projection='ortho',lat_0=45,lon_0=30)
# define lon lat mesh and data array
# assume data.shape == (N,M) and lon.shape == lat.shape == (N+1, M+1)
x,y = m(lon,lat)
ii = (
     (x[:-1,:-1] > 1e20) |
     (x[1:,:-1] > 1e20) |
     (x[:-1,1:] > 1e20) |
     (x[1:,1:] > 1e20) |
     (y[:-1,:-1] > 1e20) |
     (y[1:,:-1] > 1e20) |
     (y[:-1,1:] > 1e20) |
     (y[1:,1:] > 1e20)
     )
data[ii] = np.nan
cs = m.pcolormesh(x, y, data)

This works for my test case.

[2sn]

suggested change

    @_transform
    def pcolormesh(self,x,y,data,**kwargs):
        """
        Make a pseudo-color plot over the map
        (see matplotlib.pyplot.pcolormesh documentation).

        If ``latlon`` keyword is set to True, x,y are intrepreted as
        longitude and latitude in degrees.  Data and longitudes are
        automatically shifted to match map projection region for cylindrical
        and pseudocylindrical projections, and x,y are transformed to map
        projection coordinates. If ``latlon`` is False (default), x and y
        are assumed to be map projection coordinates.

        Extra keyword ``ax`` can be used to override the default axis instance.

        Other \**kwargs passed on to matplotlib.pyplot.pcolormesh.

        Note: (taken from matplotlib.pyplot.pcolor documentation)
        Ideally the dimensions of x and y should be one greater than those of data;
        if the dimensions are the same, then the last row and column of data will be ignored.
        """
        ax, plt = self._ax_plt_from_kw(kwargs)
        self._save_use_hold(ax, kwargs)
        try:
            ret =  ax.pcolormesh(x,y,data,**kwargs)
        finally:
            self._restore_hold(ax)
        # reset current active image (only if pyplot is imported).
        if plt:
            plt.sci(ret)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # clip to map limbs
        ret,c = self._cliplimb(ax,ret)
        if self.round:
            # for some reason, frame gets turned on.
            ax.set_frame_on(False)
        return ret

to

    @_transform
    def pcolormesh(self,x,y,data,**kwargs):
        """
        Make a pseudo-color plot over the map
        (see matplotlib.pyplot.pcolormesh documentation).

        If ``latlon`` keyword is set to True, x,y are intrepreted as
        longitude and latitude in degrees.  Data and longitudes are
        automatically shifted to match map projection region for cylindrical
        and pseudocylindrical projections, and x,y are transformed to map
        projection coordinates. If ``latlon`` is False (default), x and y
        are assumed to be map projection coordinates.

        Extra keyword ``ax`` can be used to override the default axis instance.

        Other \**kwargs passed on to matplotlib.pyplot.pcolormesh.

        Note: (taken from matplotlib.pyplot.pcolor documentation)
        Ideally the dimensions of x and y should be one greater than those of data;
        if the dimensions are the same, then the last row and column of data will be ignored.
        """
        ax, plt = self._ax_plt_from_kw(kwargs)
        # fix for invalid grid points
        if ((np.any(x > 1e20) or np.any (y > 1e20)) and
            len(x.shape) == 2 and len(y.shape) == 2):
            if not x.shape == y.shape:
                raise Exception('pcolormesh: x and y need same dimension')
            nx,ny = x.shape
            if nx < data.shape[0] or ny < data.shape[1]:
                raise Exception('pcolormesh: data dimension needs to be at least that of x and y.')
            mask = (
                (x[:-1,:-1] > 1e20) |
                (x[1:,:-1] > 1e20) |
                (x[:-1,1:] > 1e20) |
                (x[1:,1:] > 1e20) |
                (y[:-1,:-1] > 1e20) |
                (y[1:,:-1] > 1e20) |
                (y[:-1,1:] > 1e20) |
                (y[1:,1:] > 1e20)
                )
            # we do not want to overwrite original array
            data = data[:nx-1,:ny-1].copy()
            data[mask] = np.nan
        self._save_use_hold(ax, kwargs)
        try:
            ret =  ax.pcolormesh(x,y,data,**kwargs)
        finally:
            self._restore_hold(ax)
        # reset current active image (only if pyplot is imported).
        if plt:
            plt.sci(ret)
        # set axes limits to fit map region.
        self.set_axes_limits(ax=ax)
        # clip to map limbs
        ret,c = self._cliplimb(ax,ret)
        if self.round:
            # for some reason, frame gets turned on.
            ax.set_frame_on(False)
        return ret

see #476

[molinav]

After merging pull request #476 this problem should be fixed now in the master branch.

Thanks @2sn!