I am converting an application that uses matplotlib's toolkit Basemap to using Cartopy in preparation for moving from Python 2 to Python 3. I have found similar functions in Cartopy for Basemap's 'addcyclic()' and 'maskoceans()', However I cannot find something similar in either numpy or Cartopy for Basemap's shiftgrid() function.
This is the code using Basemap: '''
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import numpy as np
from mpl_toolkits.basemap import shiftgrid
bmap = Basemap(projection='ortho', lat_0=0, lon_0=0)
lons = np.arange(30, 410, 30)
lons[1] = 70
lats = np.arange(0, 100, 10)
data = np.indices((lats.shape[0], lons.shape[0]))
data = data[0] + data[1]
data, lons = shiftgrid(180., data, lons, start=False)
llons, llats = np.meshgrid(lons, lats)
x, y = bmap(llons, llats)
bmap.contourf(x, y, data)
bmap.drawcoastlines()
'''
The initial data: data '''
[[ 0 1 2 3 4 5 6 7 8 9 10 11 12]
[ 1 2 3 4 5 6 7 8 9 10 11 12 13]
[ 2 3 4 5 6 7 8 9 10 11 12 13 14]
[ 3 4 5 6 7 8 9 10 11 12 13 14 15]
[ 4 5 6 7 8 9 10 11 12 13 14 15 16]
[ 5 6 7 8 9 10 11 12 13 14 15 16 17]
[ 6 7 8 9 10 11 12 13 14 15 16 17 18]
[ 7 8 9 10 11 12 13 14 15 16 17 18 19]
[ 8 9 10 11 12 13 14 15 16 17 18 19 20]
[ 9 10 11 12 13 14 15 16 17 18 19 20 21]]
lons
[ 30 70 90 120 150 180 210 240 270 300 330 360 390]
After the 'data, lons = shiftgrid(180., data, lons, start=False)':
data
[[ 5 6 7 8 9 10 11 12 1 2 3 4 5]
[ 6 7 8 9 10 11 12 13 2 3 4 5 6]
[ 7 8 9 10 11 12 13 14 3 4 5 6 7]
[ 8 9 10 11 12 13 14 15 4 5 6 7 8]
[ 9 10 11 12 13 14 15 16 5 6 7 8 9]
[10 11 12 13 14 15 16 17 6 7 8 9 10]
[11 12 13 14 15 16 17 18 7 8 9 10 11]
[12 13 14 15 16 17 18 19 8 9 10 11 12]
[13 14 15 16 17 18 19 20 9 10 11 12 13]
[14 15 16 17 18 19 20 21 10 11 12 13 14]]
lons
[-180 -150 -120 -90 -60 -30 0 30 70 90 120 150 180]
''' I have tried the following cartopy code to recreate what the Basemap shiftgrid did. This is the Cartopy code, some things are commented out as I tried them at one time: '''
DATA_CRS = ccrs.PlateCarree()
lons = np.arange(30, 410, 30)
lons[1] = 70
lats = np.arange(0, 100, 10)
data = np.indices((lats.shape[0], lons.shape[0]))
data = data[0] + data[1]
# data2 = np.roll(data, -5)
# lons2 = np.mod(lons2 - 180.0, 360.0) - 180.0
cm_lon = 0
#llons, llats = np.meshgrid(lons2, lats)
llons, llats = np.meshgrid(lons, lats)
PROJECTION = ccrs.Orthographic(central_longitude=cm_lon)
fig1 = plt.figure(num=1, figsize=(11, 8.5), dpi=150)
ax = plt.axes(projection=PROJECTION)
ax.add_feature(cfeature.COASTLINE, linewidths=0.7)
ax.add_feature(cfeature.BORDERS, edgecolor='black', linewidths=0.7)
ax.contourf(llons, llats, data, transform=ccrs.PlateCarree())
'''
The data and the longitudes as original and I just used the 'central_longitude' in the projection. The Basemap image shows the entire globe but the Cartopy image only shows from the equator up. The color of the data seems similar except for the far right side, so I'm concerned the data didn't map the same in Cartopy as it did in Basemap.
So, the question is... Is there anything equivalent to Basemap's shiftgrid() or do I need to figure out something similar to Basemap's shiftgrid() or just use the 'central_longitude' in the projection? I don't seem to be able to paste the .png files. Any help is really appreciated. I have searched the web looking for equivalent functions but haven't found one for the shiftgrid(). Thank you.
This must be the most inelegant solution, but what I have been doing for several of Basemap's useful features that are not (yet?) in cartopy, is just to copy the function definitions from Basemap's source code. It works fine. For example, shiftgrid:
def shiftgrid(lon0,datain,lonsin,start=True,cyclic=360.0):
"""
Shift global lat/lon grid east or west.
.. tabularcolumns:: |l|L|
============== ====================================================
Arguments Description
============== ====================================================
lon0 starting longitude for shifted grid
(ending longitude if start=False). lon0 must be on
input grid (within the range of lonsin).
datain original data with longitude the right-most
dimension.
lonsin original longitudes.
============== ====================================================
.. tabularcolumns:: |l|L|
============== ====================================================
Keywords Description
============== ====================================================
start if True, lon0 represents the starting longitude
of the new grid. if False, lon0 is the ending
longitude. Default True.
cyclic width of periodic domain (default 360)
============== ====================================================
returns ``dataout,lonsout`` (data and longitudes on shifted grid).
"""
if np.fabs(lonsin[-1]-lonsin[0]-cyclic) > 1.e-4:
# Use all data instead of raise ValueError, 'cyclic point not included'
start_idx = 0
else:
# If cyclic, remove the duplicate point
start_idx = 1
if lon0 < lonsin[0] or lon0 > lonsin[-1]:
raise ValueError('lon0 outside of range of lonsin')
i0 = np.argmin(np.fabs(lonsin-lon0))
i0_shift = len(lonsin)-i0
if ma.isMA(datain):
dataout = ma.zeros(datain.shape,datain.dtype)
else:
dataout = np.zeros(datain.shape,datain.dtype)
if ma.isMA(lonsin):
lonsout = ma.zeros(lonsin.shape,lonsin.dtype)
else:
lonsout = np.zeros(lonsin.shape,lonsin.dtype)
if start:
lonsout[0:i0_shift] = lonsin[i0:]
else:
lonsout[0:i0_shift] = lonsin[i0:]-cyclic
dataout[...,0:i0_shift] = datain[...,i0:]
if start:
lonsout[i0_shift:] = lonsin[start_idx:i0+start_idx]+cyclic
else:
lonsout[i0_shift:] = lonsin[start_idx:i0+start_idx]
dataout[...,i0_shift:] = datain[...,start_idx:i0+start_idx]
return dataout,lonsout