Basemap is a module provided by mpl_toolkits tool box, to overlay the 2D gridded plot with a background of continental outlines.
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
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¶
- Prepare your 2D (lat-lon) data (e.g., from nc file)
- Generate lat-lon meshgrid (2D array)
- m = Basemap(...) decide projection, range of map
- m.drawcoastlines(), m.drawcountries(), latitude/longitude lines, etc
- cx,cy = m(lon2,lat2) # convert 2D lat/lon to map projection coordinate (!!!! Important !!!)
- m.contour(cx,cy,data), m.pcolormesh(cx,cy,data), plotting contour of data.
- add color bar, title, etc.
- plt.savefigs(), plt.show()
