%matplotlib inline
Pier Fogli asked:
I'm working with data on the cubed sphere grid:
.
Even though the data is interpolated on a regular lat-lon grid for plotting, I would like to plot the edges of the 6 cube faces for visual reference.
Moreover I would like this to be independent from the map projection type and map center lon/lat
But all I can get using Basemap is:
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import numpy as np
def plot_cube(mp):
[y,x]=mp.gcpoints(45,-45,45,45,100)
mp.plot(y,x,'b',linewidth=1.5)
[y,x]=mp.gcpoints(135,-45,135,45,100)
mp.plot(y,x,'b',linewidth=1.5)
[y,x]=mp.gcpoints(225,-45,225,45,100)
mp.plot(y,x,'b',linewidth=1.5)
[y,x]=mp.gcpoints(315,-45,315,45,100)
mp.plot(y,x,'b',linewidth=1.5)
[y,x]=mp.gcpoints(45,-45,135,-45,100)
mp.plot(y,x,'b',linewidth=1.5)
[y,x]=mp.gcpoints(135,-45,225,-45,100)
mp.plot(y,x,'b',linewidth=1.5)
[y,x]=mp.gcpoints(225,-45,315,-45,100)
mp.plot(y,x,'b',linewidth=1.5)
[y,x]=mp.gcpoints(-45,-45,45,-45,100)
mp.plot(y,x,'b',linewidth=1.5)
#
[y,x]=mp.gcpoints(45,45,135,45,100)
mp.plot(y,x,'b',linewidth=1.5)
[y,x]=mp.gcpoints(135,45,225,45,100)
mp.plot(y,x,'b',linewidth=1.5)
[y,x]=mp.gcpoints(225,45,315,45,100)
mp.plot(y,x,'b',linewidth=1.5)
[y,x]=mp.gcpoints(-45,45,45,45,100)
mp.plot(y,x,'b',linewidth=1.5)
fh=plt.figure(figsize=(8, 8))
mp1=Basemap(projection='cyl',lon_0=180.0)
mp1.drawcoastlines()
mp1.drawparallels(np.arange(-90.,90.01,30.0))
mp1.drawmeridians(np.arange(0.,360.01,60.0))
plot_cube(mp1)
plt.title('Equidistant Cylindrical [0 360]')
plt.show()
Which is definitely not what I would like to do.
Does anyone have any advice on how to solve this issue?
Answer¶
The fundamental problem is that you're converting line end-points into the coordinates of the map, and in doing so are disregarding the topology of the actual line. This is the fundamental problem with your code, and indeed with the approach encouraged with basemap.
It might be worth taking a look at cartopy, which has been designed from the ground-up to avoid the assumed topology problem.
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import numpy as np
def cubedsphere_coords():
"""
Return the ``lons, lats`` for the 6 faces of the cubed-sphere.
"""
lons = np.array([45, 135, 135, 45, 45])
lats = np.array([-45, -45, 45, 45, -45])
# Add the 4 easy faces.
coords = [[lons + face * 90, lats] for face in range(4)]
# The north pole face.
coords.append([[-135, -45, 45, 135, -135], [45] * 5])
# And the south pole face, in reverse order to the north pole patch.
coords.append([[135, 45, -45, -135, 135], [-45] * 5])
return coords
The key with cartopy is that we can define coordinates in any coordinate reference system, and any points, lines and polygons will be transformed appropriately (splitting geometries where necessary).
def plot_cubed_sphere_lines(ax):
for lons, lats in cubedsphere_coords():
ax.plot(lons, lats, transform=ccrs.Geodetic(), color='black')
plt.figure(figsize=(10, 10))
target_projection = ccrs.Orthographic(central_longitude=75, central_latitude=20)
ax = plt.axes(projection=target_projection)
ax.coastlines()
plot_cubed_sphere_lines(ax)
ax.gridlines()
ax.set_global()
ax.stock_img()
plt.show()
Drawing filled polygons is slightly more involved, and requires the construction of a Polygon geometry, but the principle is the same - define the coordinates in the coordinate system that they are in, and let Cartopy do the transformation for you:
import shapely.geometry as sgeom
def plot_cube_sphere_filled(ax):
colors = ['red', 'purple', 'green', 'orange', 'blue', 'yellow']
coords = cubedsphere_coords()
for color, (lons, lats) in zip(colors, cubedsphere_coords()):
poly = sgeom.Polygon(zip(lons, lats))
ax.add_geometries([poly], ccrs.Geodetic(), facecolor=color, zorder=0)
plt.figure(figsize=(10, 10))
target_projection = ccrs.Orthographic(central_longitude=75, central_latitude=20)
ax = plt.axes(projection=target_projection)
ax.coastlines()
plot_cube_sphere_filled(ax)
ax.gridlines()
ax.set_global()
plt.show()
