Problem with ortho projection and pcolormesh
Closed this issue · 4 comments
ricitron commented
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 commented
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 commented
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 commented
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 retto
@_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 retsee #476