Basemap 在2D圖形背景加上海陸地圖

Basemap is a module provided by mpl_toolkits tool box, to overlay the 2D gridded plot with a background of continental outlines.

Basemap is available on study server. Some of you may not have Basemap installed in the python system on your computer. (check this by "import mpl_toolkits.basemap"). To install basemap, please see the steps here:

Reference for basemap package:

To import Basemap:

In [6]:
from mpl_toolkits.basemap import Basemap
# other packages we need 
import netCDF4 as nc   # read nc file
import numpy as np     # processing data array
import matplotlib.pyplot as plt   # plotting contour
import as cm        # color map

As an example, we will plot the global map of temperature from the NetCDF file:

In [2]:
rootgrp = nc.Dataset("")  # see example in Lecture NetCDF

# extract the lat, lon, pressure level, and time information
lat = rootgrp.variables['lat'][:]
lon = rootgrp.variables['lon'][:]
lev = rootgrp.variables['lev'][:]
time = rootgrp.variables['time'][:] # get value of time
time_u=rootgrp.variables['time'].units # get unit 
time_d=nc.num2date(time,units = time_u,calendar='standard') # get date using NetCDF num2date function

lon2,lat2=np.meshgrid(lon,lat) # convert 1D lat/lon to 2D lat/lon arrays

# extract global temperature at first time step at first level
t = rootgrp.variables['tmp'][0,0,:,:] # first time step (0), first level (0), all lat/lon (:)


m = Basemap(projection='...',lon_0=...,lat_0=..., llcrnrlat=..., llcrnrlon=...,urcrnrlat=..., urcrnrlon=...)

  • projection: 投影法
  • lon_0,lat_0: 圖中心經緯度 lon/lat at the center of map
  • llcrnrlat, llcrnrlon: 圖左下角經緯度 lon/lat at the lower-left corner
  • urcrnrlat, urcrnrlon: 圖右上角經緯度 lon/lat at the upper-right corner
In [3]:
# Ex1. Global map, basic projection and layout
m = Basemap(projection='cyl',lon_0=180,lat_0=0)

m.drawcoastlines()   #畫海岸線

cx,cy =m(lon2,lat2)  # convert to map projection coordinate
CS = m.pcolormesh(cx,cy,t,cmap=cm.jet,shading='gouraud')
plt.title('Cylindrical, T at 1000 hPa, GFS'+date)
In [4]:
# Ex2. Global map, orthographic projection
m = Basemap(projection='ortho',lon_0=120,lat_0=25)

m.drawcoastlines()   #畫海岸線

parallels = np.arange(-90.,90,30.)
m.drawparallels(parallels) # draw parallels 畫緯度線
meridians = np.arange(0.,360.,20.)
m.drawmeridians(meridians) # draw meridians 畫經度線

cx2,cy2 =m(lon2,lat2)  # convert to map projection coordinate
CS2 = m.pcolormesh(cx2,cy2,t,cmap=cm.jet,shading='gouraud')
plt.title('Orthographic, T at 1000 hPa, GFS'+date)
In [5]:
# plot Asia Only, maskout continents, using stereographic projection
m = Basemap(projection='stere',lon_0=120,lat_0=20,llcrnrlat=-10, 
            llcrnrlon=50,urcrnrlat=50, urcrnrlon=190)

m.drawcoastlines()   #畫海岸線
m.drawcountries()    #畫國界

m.drawparallels(parallels,labels=[1,1,0,0],fontsize=10)  # 緯度度線、在左右兩邊加緯度標籤
m.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)  # 經度線、在下方加經度標籤

cx3,cy3 =m(lon2,lat2)  # convert to map projection coordinate
clevs = np.arange(260,310,5)
CS3 = m.contourf(cx3,cy3,t,clevs,cmap=cm.jet)

plt.title('Stereographic, T at 1000 hPa, GFS'+date)

Summarize the steps to draw data on a map

  1. Prepare your 2D (lat-lon) data (e.g., from nc file)
  2. Generate lat-lon meshgrid (2D array)
  3. m = Basemap(...) decide projection, range of map
  4. m.drawcoastlines(), m.drawcountries(), latitude/longitude lines, etc
  5. cx,cy = m(lon2,lat2) # convert 2D lat/lon to map projection coordinate (!!!! Important !!!)
  6. m.contour(cx,cy,data), m.pcolormesh(cx,cy,data), plotting contour of data.
  7. add color bar, title, etc.
  8. plt.savefigs(),