Plotting cylindrical map data over a 3D sphere in

Question

Suppose I have an image with a cylindrical map of a planet, say one of these:

http://www.johnstonsarchive.net/spaceart/cylmaps.html

And I want to plot it over a 3D sphere to recover the original figure of the planet.

Is there a way to do it using a Python package like matplotlib, mayavi, basemap or similar?

Answer 1

the basemap toolkit is ideal for this task.

The Problem seems to be that you don't really have the data, i.e. something like lat,lon,value for each pixel. The problem is also that for a given image, you normally don't know the projection which was used to create the image, therefore you cannot do the back-transformation to get the original data.

If you just like to plot a image of the earth, use the bluemarble() function from the basemap toolkit:

http://wiki.scipy.org/Cookbook/Matplotlib/Maps

Answer 2

Update: This is the new version using Cartopy, as basemap is EOL. Below is the original answer.

import matplotlib.pyplot as plt
import cartopy.crs as ccrs

img = plt.imread("/tmp/venuscyl4.tif")

plt.figure(figsize=(3, 3))

ax = plt.axes(projection=ccrs.Orthographic(-10, 45))
ax.gridlines(color='black', linestyle='dotted')
ax.imshow(img, origin="upper", extent=(-180, 180, -90, 90),
          transform=ccrs.PlateCarree())  # Important

plt.show()

---

Thanks to Raphael Roth answer I finally found what I was looking for: the basemap method warpimage.

Here is a very minimal example. Using this cylindrical map of Venus, and based on the simple example of the cookbook:

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
# set up orthographic map projection with
# perspective of satellite looking down at 50N, 100W.
# use low resolution coastlines.
# don't plot features that are smaller than 1000 square km.
bmap = Basemap(projection='ortho', lat_0 = 50, lon_0 = -100,
              resolution = 'l', area_thresh = 1000.)
# plot surface
bmap.warpimage(image='venuscyl4.jpg')
# draw the edge of the map projection region (the projection limb)
bmap.drawmapboundary()
# draw lat/lon grid lines every 30 degrees.
bmap.drawmeridians(np.arange(0, 360, 30))
bmap.drawparallels(np.arange(-90, 90, 30))
plt.show()

produces the following output: