CODEWITHSUNDEEP
matplotlibmap-projectionsmatplotlib-basemap
alex_reader
alex_reader

Reputation: 711

Plotting sphere grid with matplotlib basemap

I would like to draw a grid covering all the sphere on an orthographic projection. The issue is cells outside the projection are not drawed correctly. This happened with drawgreatcircles as pointed here.

I have also tried to use Polygons as described here, but same problem.

Finally, I have coded a custom check based on Wikipedia. The idea is for each point of each segment, we check cos c (cf Wikipedia) and do not plot it if the cosinus is negative.

My question is : can we do this kind of check with basemap own functions ? This strategy would not work for other projections.

Also, why is this kind of check not included in Basemap ?

Upvotes: 2

Views: 1816

Answers (2)

pelson
pelson

Reputation: 21839

Thanks to your example, I took the data and plotted it with cartopy. The following changes were needed to create the plot:

import cartopy.crs as ccrs
ax =plt.axes(projection=ccrs.Orthographic())
plt.pcolormesh(lons, lats,val, edgecolors='k', 
               linewidths=1, transform=ccrs.PlateCarree())

ax.coastlines()
ax.gridlines()

plt.show()

output using cartopy

This is using pcolormesh so is pretty quick (though your example wasn't that slow on my machine in the first place).

Upvotes: 3

alex_reader
alex_reader

Reputation: 711

Here is a solution using pcolor :

import pylab as plt
from mpl_toolkits.basemap import Basemap
import numpy as np

nb_lat2 = 20
nb_lat = 2*nb_lat2
nb_lon = 3*(2*(nb_lat+1) - 1)
lats = np.zeros((2*nb_lat, nb_lon))
lons = np.zeros((2*nb_lat, nb_lon))
val = np.zeros((2*nb_lat, nb_lon))
dlat = 90./nb_lat2
for i in range(nb_lat):

  nb_lon = 2*(i+1)-1
  if ((i+1) > nb_lat2):
    nb_lon = 2*(nb_lat - i)-1

  dlon = 120./nb_lon
  lats[2*i][:] = 90 - i*dlat
  lats[2*i+1][:] = 90 - (i+1)*dlat
  for j in range(nb_lon):
    lons[2*i][j] = j*dlon
    lons[2*i+1][j] = j*dlon

    for k in range(1,3):
      lons[2*i][j + k*nb_lon] = j*dlon + 120.*k
      lons[2*i+1][j + k*nb_lon] = j*dlon + 120.*k

    lons[2*i][3*nb_lon:] = nb_lon*dlon + 240.
    lons[2*i+1][3*nb_lon:] = nb_lon*dlon + 240.

lons = lons - 180
val = lats + lons
# Crash
##m = Basemap(projection='robin',lon_0=0,resolution=None)

#m = Basemap(projection='mill',lon_0=0)
m = Basemap(projection='ortho', lat_0=0,lon_0=0)
x, y = m(lons, lats)
m.pcolor(x,y,val, edgecolors='k', linewidths=1)

m.drawcoastlines()
m.drawparallels(np.arange(-90.,91.,30.))
m.drawmeridians(np.arange(-180.,181.,60.))

plt.show()

This does exactly what I want : drawing rectangles and filling them with one color. But it is very slow (too slow). A lot of cells are unused : at the end of a latidude line, we set the width of unused cells to 0.

Another issue is some projections crash (Robin for example).

Upvotes: 1

Related Questions