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:http://homepage.ntu.edu.tw/~weitingc/fortran_lecture/Lecture_P_4_Installing_Basemap_NetCDF4.slides.html

Reference for basemap package: https://matplotlib.org/basemap/users/index.html

To import Basemap:

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 matplotlib.cm as cm        # color map

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

rootgrp = nc.Dataset("GFS20171217_tmp.nc")  # 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 (:)

Basemap

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

• projection: 投影法
  • cyl: direct latitude/longitude (equidistance cylindrical)
  • ortho: orthographic
  • stereo: stereographic
  • many other projection options (each may have specific syntax to control the layout):
    • https://matplotlib.org/basemap/users/mapsetup.html
• 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

# 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.colorbar(CS,orientation='horizontal')
date=time_d[0].strftime('%Y/%m/%d')  
plt.title('Cylindrical, T at 1000 hPa, GFS'+date)
plt.show()

# 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.colorbar(CS,orientation='vertical')
plt.title('Orthographic, T at 1000 hPa, GFS'+date)
plt.show()

# 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.fillcontinents(color='grey')

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.colorbar(CS3,orientation='horizontal').set_label('K')
plt.title('Stereographic, T at 1000 hPa, GFS'+date)
plt.show()

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(), plt.show()