一点总结：Ubuntu+anaconda+Python+basemap(+WRF)

8828 · 发表于 2016-11-29 19:10:30

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

本帖最后由 8828 于 2017-7-14 16:55 编辑

看到我们的Python板块有了关于气象应用的帖子，我也想来分享一些东西。这是将我周末弄得一些东西，系统整理一下。(从零开始，有很多不足，刚开始接触Python的气象应用)
outline：
1、Python科学计算环境的搭建
2、Python气象库的应用例子
3、WRF-Python处理

一、Python科学计算环境的搭建
#推荐两种环境：①canopy
②anaconda
#这里主要应用anaconda
#个人使用平台： Ubuntu16.04(个人科研学习拒绝Windows)

#声明：所有的下载渠道配置到清华镜像上（除个别库外优先选择该channel）

1.下载客户端：客户端传送门
2.开始搭建：
#我开始安装的是miniconda，详细问题百度

$bash *.sh # *.sh为下载的文件

复制代码

$conda #检测是否安装成功
$conda create -n basemap python=2 #创建虚拟环境
###配置清华镜像
$conda config --add channels 'https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/'
$conda config --set show_channel_urls yes
###云大镜像正在编写程序中，欢迎到云大开源中心访问
###气象所需库：basemap （属于matplotlib）
$conda install basemap
###自动配置相关所需库，使用镜像快的飞起。环境算基本搭建完毕，详细内容百度，如遇错误，最好的办法是研究报错信息
###下面介绍basemap的例子，涉及到WRF的内容在第三部分

复制代码

二、应用举例
###由于图片放入代码模式里，代码中只显示图片代码，效果图在最后

#来自相关网站，个人修改，相关网址最后总结
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
import pdb
# 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.
map = Basemap(projection='ortho', lat_0 = 35, lon_0 = 135,
resolution = 'l', area_thresh = 1000.)
# draw coastlines, country boundaries, fill continents.
"""map.drawcoastlines()
map.drawcountries()
map.fillcontinents(color = 'coral')"""
#NASA instead of baove, conda install PIL
map.bluemarble()
# draw the edge of the map projection region (the projection limb)
map.drawmapboundary()
# draw lat/lon grid lines every 30 degrees.
map.drawmeridians(np.arange(0, 360, 30))
map.drawparallels(np.arange(-90, 90, 30))
# lat/lon coordinates of five cities.
"""lats = [40.02, 32.73, 38.55, 48.25, 17.29]
lons = [-105.16, -117.16, -77.00, -114.21, -88.10]
cities=['Boulder, CO','San Diego, CA',
'Washington, DC','Whitefish, MT','Belize City, Belize']"""
#cities
lats = [25.05, 39.91]
lons = [102.73, 116.41]
cities = ['Kunming', 'Tangshan']
# compute the native map projection coordinates for cities.
x,y = map(lons,lats)
# plot filled circles at the locations of the cities.
map.plot(x,y,'bo')
# plot the names of those five cities.
for name,xpt,ypt in zip(cities,x,y):
plt.text(xpt+50000,ypt+50000,name)
# make up some data on a regular lat/lon grid.
nlats = 73; nlons = 145; delta = 2.*np.pi/(nlons-1)
lats = (0.5*np.pi-delta*np.indices((nlats,nlons))[0,:,:])
lons = (delta*np.indices((nlats,nlons))[1,:,:])
wave = 0.75*(np.sin(2.*lats)**8*np.cos(4.*lons))
mean = 0.5*np.cos(2.*lats)*((np.sin(2.*lats))**2 + 2.)
# compute native map projection coordinates of lat/lon grid.
x, y = map(lons*180./np.pi, lats*180./np.pi)
# contour data over the map.
CS = map.contour(x,y,wave+mean,15,linewidths=1.5)
#save figure
plt.savefig('city.png', bbox_inches='tight')
plt.show()
<img src="http://bbs.06climate.com/forum.php?mod=image&aid=59639&size=300x300&key=b5af3a6ca1216f08&nocache=yes&type=fixnone" aid="attachimg_59639" alt="" width="433" border="0" height="388">

复制代码

###Python与NASA数据的结合

###当时保存的图片是前几天的，利用Python的datetime模块可获得此刻时间
import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from datetime import datetime
# miller projection
map = Basemap(projection='mill',lon_0=180)
# plot coastlines, draw label meridians and parallels.
map.drawcoastlines()
map.drawparallels(np.arange(-90,90,30),labels=[1,0,0,0])
map.drawmeridians(np.arange(map.lonmin,map.lonmax+30,60),labels=[0,0,0,1])
# fill continents 'coral' (with zorder=0), color wet areas 'aqua'
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')
# shade the night areas, with alpha transparency so the
# map shows through. Use current time in UTC.
date = datetime.utcnow()
CS=map.nightshade(date)
plt.title('Day/Night Map for %s (UTC)' % date.strftime("%d %b %Y %H:%M:%S"))
plt.savefig('day_night' , bbox_inches='tight')
plt.show()
<img src="http://bbs.06climate.com/forum.php?mod=image&aid=59640&size=300x300&key=7c7d93cc3e23dbab&nocache=yes&type=fixnone" aid="attachimg_59640" alt="" width="416" border="0" height="314">

复制代码

#####Python库很强大##########

#Python的模块可以在网上抓取数据，进而直接提取数据
from mpl_toolkits.basemap import Basemap
from netCDF4 import Dataset, date2index
import numpy as np
import matplotlib.pyplot as plt
from datetime import datetime
date = datetime(2007,01,01,0) # date to plot.
# open dataset.
dataset = \
Dataset('http://www.ncdc.noaa.gov/thredds/dodsC/OISST-V2-AVHRR_agg')
timevar = dataset.variables['time']
timeindex = date2index(date,timevar) # find time index for desired date.
# read sst. Will automatically create a masked array using
# missing_value variable attribute. 'squeeze out' singleton dimensions.
sst = dataset.variables['sst'][timeindex,:].squeeze()
# read ice.
ice = dataset.variables['ice'][timeindex,:].squeeze()
# read lats and lons (representing centers of grid boxes).
lats = dataset.variables['lat'][:]
lons = dataset.variables['lon'][:]
lons, lats = np.meshgrid(lons,lats)
# create figure, axes instances.
fig = plt.figure()
ax = fig.add_axes([0.05,0.05,0.9,0.9])
# create Basemap instance.
# coastlines not used, so resolution set to None to skip
# continent processing (this speeds things up a bit)
m = Basemap(projection='kav7',lon_0=0,resolution=None)
# draw line around map projection limb.
# color background of map projection region.
# missing values over land will show up this color.
m.drawmapboundary(fill_color='0.3')
# plot sst, then ice with pcolor
im1 = m.pcolormesh(lons,lats,sst,shading='flat',cmap=plt.cm.jet,latlon=True)
im2 = m.pcolormesh(lons,lats,ice,shading='flat',cmap=plt.cm.gist_gray,latlon=True)
# draw parallels and meridians, but don't bother labelling them.
m.drawparallels(np.arange(-90.,99.,30.))
m.drawmeridians(np.arange(-180.,180.,60.))
# add colorbar
cb = m.colorbar(im1,"bottom", size="5%", pad="2%")
# add a title.
ax.set_title('SST and ICE analysis for %s'%date)
# save figure
plt.savefig('sst1.png', bbox_inches='tight')
plt.show()
<img src="http://bbs.06climate.com/forum.php?mod=image&aid=59641&size=300x300&key=eab0e1429a062936&nocache=yes&type=fixnone" aid="attachimg_59641" alt="" border="0">

复制代码

############展示一些其他例子，原脚本来自官网，做个人修改##########
SLP and Wind_.png

#########一些介绍############
basemap是matplotlib的一个库，matplotlib是科学计算中的重要库之一，关于matplotlib请百度，
basemap tutorial 在github上均有，请利用好网络资源，百度即可，不再提供网址
##########################

三、WRF-Python
###这里要说明，Python作图是NCL官网提供的，官方手册也有介绍，在最后的几部分####
###WRF应用时需要的库是pynio，用来读取wrfout数据，请Google

>>>import Nio #用于读取variables、global variables等

复制代码

#########################我遇到的问题：############################################
#ncar提供pynio这个库，但版本上与basemap1.0.7冲突，即geoslib3.3.2,与geoslib3.4.2冲突。可以选择其他channel，速度很慢，
#不能保证冲突问题可以解决，也许降低basemap版本可以消除冲突

$conda install -c ncar pynio #如遇问题请百度和参考报错信息
$anaconda search -t conda package #在anaconda客户端搜索package和提供package的channe（可以利用conda安装anaconda client）

复制代码

#########################处理办法################################################
###利用Python的netCDF4模块读取wrfout的变量。目前只读到了variables，global variables还没有读取到，以后会认真找办法

#原脚本来自NCL官网，个人修改。处理basemap与Nio矛盾的问题在脚本中体现出来
#----------------------------------------------------------------------
# This script calculates "TK" from data read off a WRF-ARW output,
# and creates a matplotlib plot over a lambert conformal map. It
# checks first that the WRF data is on a Lambert Conformal projection.
#
# See wrf_plot_TK_shapefiles_mpl.py for an example of adding
# shapefile outlines.
#
# You have to provide this script with a WRF output file:
#
# python wrf_plot_TK_mpl.py -f ../Data/wrfout_d03_2012-04-22_23_00_00
#
#----------------------------------------------------------------------
import os, sys, argparse
import numpy as np
from netCDF4 import Dataset as NetCDFFile #import Nio
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.basemap import Basemap
def wrf_tk(pressure,theta):
"""
Calculates temperature in [K] from ARW WRF model output.
If variable is more than two dimensions, then only the first element
of the leftmost dimensions are plotted.
tk = wrf_tk(pressure,theta)
pressure - full pressure (perturbation + base state pressure) in Pa
theta - potential temperature (perturbation + reference temperature) in K
"""
p1000mb = 100000.
r_d = 287.
cp = 7.*r_d/2.
pi = (pressure/p1000mb)**(r_d/cp)
tk = pi*theta
return(tk)
filename = "wrfout_d02_2016-01-01_00_00_00"
if os.path.exists(filename):
ncfile = NetCDFFile('wrfout_d02_2016-01-01_00_00_00') #Nio.open_file(filename + ".nc")
else:
raise RuntimeError("'%s' not found" % filename)
# Get the dimensions
nx = 100 #ncfile.dimensions["west_east"]
ny = 100 #ncfile.dimensions["south_north"]
nz = 30 #ncfile.dimensions["bottom_top"]
# Get the map projection information
map_proj = 1 #ncfile.attributes["MAP_PROJ"]
cen_lon = 103.49 #ncfile.attributes["STAND_LON"]
cen_lat = 36.03 #ncfile.attributes["MOAD_CEN_LAT"]
truelat1 = 30 #ncfile.attributes["TRUELAT1"]
truelat2 = 60 #ncfile.attributes["TRUELAT2"]
dx = 10000 #ncfile.attributes["DX"]
dy = 10000 #ncfile.attributes["DY"]
if map_proj != 1:
print("Error: the data on %s do not appear to be on a lambert conformal projection." % f)
sys.exit(1)
# Get the latitude and longitude values
xlon = ncfile.variables["XLONG"][0, :, :]
xlat = ncfile.variables["XLAT"][0, :, :]
# Get the lat/lon corner points
lat_ll = xlat[0,0]
lat_ur = xlat[-1,-1]
lon_ll = xlon[0,0]
lon_ur = xlon[-1,-1]
# Create the basemap object
bm = Basemap(projection="lcc",
lon_0=cen_lon,
lat_0=cen_lat,
lat_1=truelat1,
lat_2=truelat2,
llcrnrlat=lat_ll,
urcrnrlat=lat_ur,
llcrnrlon=lon_ll,
urcrnrlon=lon_ur,
resolution='l')
# Transform lon/lat points to x/y points in projected coordinates
x, y = bm(xlon, xlat)
# Get the data for plotting
T = ncfile.variables["T"]
P = ncfile.variables["P"]
PB = ncfile.variables["PB"]
T = T[0] + 300
P = P[0] + PB[0]
TK = wrf_tk(P, T)
rank = len(TK.shape)
vname = "TK"
if(rank == 2):
var = TK[:]
elif(rank == 3):
var = TK[0,:,:]
elif(rank == 4):
var = TK[0,0,:,:]
else:
raise RuntimeError("Only 2, 3, and 4 dimensions are supported")
# Create the figure and add axes
fig = plt.figure(figsize=(8,8))
ax = fig.add_axes([0.1,0.1,0.8,0.8])
# Set contour levels
min_lev = np.amin(var)
max_lev = np.amax(var)
# Make only 10 levels
lev_spacing = (max_lev - min_lev) / 10.0
clevs = np.arange(min_lev, max_lev, lev_spacing)
# Make the contour plot
cplot = bm.contourf(x, y, var, clevs, cmap=cm.jet, extend="both")
# Define and plot the meridians and parallels
min_lat = np.amin(xlat)
max_lat = np.amax(xlat)
min_lon = np.amin(xlon)
max_lon = np.amax(xlon)
# Make only 5 parallels and meridians
parallel_spacing = (max_lat - min_lat) / 5.0
merid_spacing = (max_lon - min_lon) / 5.0
parallels = np.arange(min_lat, max_lat, parallel_spacing)
meridians = np.arange(min_lon, max_lon, merid_spacing)
bm.drawcoastlines(linewidth=1.5)
bm.drawparallels(parallels,labels=[1,0,0,0],fontsize=10)
bm.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10)
# Create the colorbar
cb = bm.colorbar(cplot,"bottom", size="5%", pad="5%")
cb.set_label(vname)
# Set the plot title
ax.set_title(vname)
#save figure
plt.savefig('TK.png', bbox_inches='tight')
# Show the figure (remove this line if only wanting to produce PNG file)
plt.show()
<img src="http://bbs.06climate.com/forum.php?mod=image&aid=59680&size=300x300&key=d9bd9cdee1b7e754&nocache=yes&type=fixnone" aid="attachimg_59680" alt="" width="288" border="0" height="305">

复制代码

############################说明:#######################################
#我在两个虚拟环境里分别安装了basemap和PYnio(来自NCAR)，大家可以尝试换用不同channel提供的PYnio
#

#我对basemap中的参数理解，在非WRF气象数据作图中应用
lon_0=cen_lon, #lon_0 中心经度
lat_0=cen_lat, #lat_0 中心纬度
lat_1=truelat1, #真实纬度1
lat_2=truelat2, #真实纬度2
llcrnrlat=lat_ll, #起始纬度
urcrnrlat=lat_ur, #终止纬度
llcrnrlon=lon_ll, #起始经度
urcrnrlon=lon_ur, #终止经度

复制代码

#
###2016.12.06

| ================ ====================================================
| Variable Name Description
| ================ ====================================================
| projection map projection. Print the module variable
| ``supported_projections`` to see a list of allowed
| values.
| epsg EPSG code defining projection (see
| http://spatialreference.org for a list of
| EPSG codes and their definitions).
| aspect map aspect ratio
| (size of y dimension / size of x dimension).
| llcrnrlon longitude of lower left hand corner of the
| selected map domain.
| llcrnrlat latitude of lower left hand corner of the
| selected map domain.
| urcrnrlon longitude of upper right hand corner of the
| selected map domain.
| urcrnrlat latitude of upper right hand corner of the
| selected map domain.
| llcrnrx x value of lower left hand corner of the
| selected map domain in map projection coordinates.
| llcrnry y value of lower left hand corner of the
| selected map domain in map projection coordinates.
| urcrnrx x value of upper right hand corner of the
| selected map domain in map projection coordinates.
| urcrnry y value of upper right hand corner of the
| selected map domain in map projection coordinates.
| rmajor equatorial radius of ellipsoid used (in meters).
| rminor polar radius of ellipsoid used (in meters).
| resolution resolution of boundary dataset being used (``c``
| for crude, ``l`` for low, etc.).
| If None, no boundary dataset is associated with the
| Basemap instance.
| proj4string the string describing the map projection that is
| used by PROJ.4.
| ================ ====================================================

复制代码

###

############netCDF4暂时代替Nio的办法###################################
#手动添加相应数值（脚本中注释部分，问题主要集中于global variables的读取），在namelist已经设定 #
################################################################

################################################################
canopy提供basemap模块，大家可以尝试一下，我没有用canopy处理以上内容

################################################
先写这么多吧，很多网站没有给出，直接百度即可，fork到github上仔细学习
user guide(其它可直接在github搜索basemap）
#################################################
uvcdat的效果也很棒，已有前辈写出，大二的新人就不乱说话了.
为云南大学开源中心做个广告，云大gitlab
#################################################

2016.11.29
Thanks&Regards

#######################2016.12.08 FOR WRF##############################
#受到Python全自动：下载数据—分析绘图—上传图到邮箱—关机（by 平流层的萝卜）的启发，这两天写了
#一个下载ERA-interim数据并作图，以及做出wrfout相应物理量的图像并发送至某一邮箱的脚本。
#参考资料：飘逸的python - 发送带各种类型附件的邮件
#                详细讲解用Python发送SMTP邮件的教程
#                Brief request syntax(ecmwf)
#                ECMWF官网
#                python 官网
#                NCL官网                                                                                                                                  #
##########################################################################

# -*- coding:utf-8 -*-
#weiyiwang@mail.ynu.edu.cn
#由于一些内存问题，具有科学意义的程序代码被我改了一下，大家可参考官网资料改回来
#Ubuntu16.04-----conda------python
"""
First, we need to download NC data from ECMWF,
Second, we need to plot data,
Third, we need to plot wrfout data,
At last, we need to transport pictures with e-mail
"""
#################################################
##########
#get data#
##########
from ecmwfapi import ECMWFDataServer
server = ECMWFDataServer()
server.retrieve({
'stream' : "oper",
'levtype' : "sfc",
'param' : "165.128/166.128",
'repres' : "sh",
'step' : "0",
'time' : "0000",
'date' : "2016-01-01/to/2016-01-01",
'dataset' : "interim",
'type' : "an",
'class' : "ei",
'grid' : "0.75/0.75",
'area' : "39/96/30/107.25",
'format' : "netcdf",
'target' : "uv.nc"
})
#########################################################################
##############################
#plot uv with data downloaded#
##############################
import numpy as np
from netCDF4 import Dataset as NetCDFFile
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from mpl_toolkits.basemap import Basemap
#load data
data = NetCDFFile('uv.nc')
#extra longitude, latitude, u, v then calculate speed
lon = data.variables['longitude']
lat = data.variables['latitude']
u = data.variables['u10'][0, :, :]
v = data.variables['v10'][0, :, :]
spd = np.sqrt(u**2+v**2)
#draw a map backgroud
m = Basemap(projection='cyl',llcrnrlat=min(lat)-1,urcrnrlat=max(lat)+1,\
resolution='l',llcrnrlon=min(lon)-1,urcrnrlon=max(lon)+1)
#set figure attributes
fig = plt.figure(figsize=(8,10))
ax = fig.add_axes([0.1,0.1,0.8,0.8])
#add lon, lat, u, v as quiver
uproj,vproj,xx,yy = m.rotate_vector(u,v,lon,lat,returnxy=True)
Q = m.quiver(xx,yy,uproj,vproj,spd,cmap=cm.jet,headlength=7)
#draw coastlines, state and country boundaries, edge of map
m.drawcoastlines()
#create and draw meridians and parallels grid lines
m.drawparallels(np.arange( -90., 90.,10.),labels=[1,0,0,0],fontsize=10)
m.drawmeridians(np.arange(-180.,180.,10.),labels=[0,0,0,1],fontsize=10)
#save
plt.savefig('uv_ecmwf.pdf')
#close
plt.close(fig)
###############################################################################
##########################
#plot uv with wrfout data#
##########################
def wrf_unstagger(x,dim):
rank = len(x.shape)
if(rank == 4 and dim == "u"):
xu = 0.5*(x[:,:,:,:-1] + x[:,:,:,1:])
elif(rank == 4 and dim == "v"):
xu = 0.5*(x[:,:,:-1,:] + x[:,:,1:,:])
if(rank == 3 and dim == "u"):
xu = 0.5*(x[:,:,:-1] + x[:,:,1:])
elif(rank == 3 and dim == "v"):
xu = 0.5*(x[:,:-1,:] + x[:,1:,:])
elif(rank == 2 and dim == "u"):
xu = 0.5*(x[:,:-1] + x[:,1:])
elif(rank == 2 and dim == "v"):
xu = 0.5*(x[:-1,:] + x[1:,:])
return xu
#Read data
a = NetCDFFile('/media/wwy/DATA/wrfoutNoah/wrfout_d02_2016-01-01_00_00_00')
u = a.variables["U"]
v = a.variables["V"]
xlatu = a.variables["XLAT_U"]
xlonu = a.variables["XLONG_U"]
# U, V, XLATU, XLONU are on different grids. Unstagger them to the same grid.
xlat = wrf_unstagger(xlatu,"u")
xlon = wrf_unstagger(xlonu,"u")
ua = u
va = v
#First timestep, lowest (bottommost) level, every 5th lat/lon
nl = 0
nt = 0
nstep = 2
u10 = ua[nt,nl,::nstep,::nstep]
v10 = va[nt,nl,::nstep,::nstep]
spd = np.sqrt(u10**2+v10**2)
lat = xlat[nt,::nstep,::nstep]
lon = xlon[nt,::nstep,::nstep]
#draw map background
m = Basemap(projection='cyl',llcrnrlat=lat.min()-1,urcrnrlat=lat.max()+1,\
resolution='l', llcrnrlon=lon.min()-1,urcrnrlon=lon.max()+1)
#create figure, add axes
fig = plt.figure(figsize=(8,10))
ax = fig.add_axes([0.1,0.1,0.8,0.8])
# rotate vectors to projection grid.
uproj,vproj,xx,yy = m.rotate_vector(u10,v10,lon,lat,returnxy=True)
Q = m.quiver(xx,yy,uproj,vproj,spd,cmap=cm.jet,headlength=7)
#draw coastlines, state and country boundaries, edge of map
m.drawcoastlines()
#create and draw meridians and parallels grid lines
m.drawparallels(np.arange( -90., 90.,10.),labels=[1,0,0,0],fontsize=10)
m.drawmeridians(np.arange(-180.,180.,10.),labels=[0,0,0,1],fontsize=10)
#save figure
plt.savefig("wrf_UV.pdf")
#close figure
plt.close(fig)
#####################################################################################
###############
#mail the pdfs#
###############
#e-mail with SMTP
import smtplib
from email import encoders
from email.mime.text import MIMEText
from email.mime.multipart import MIMEMultipart
from email.mime.application import MIMEApplication
from email.header import Header
from email.utils import parseaddr, formataddr
#definite function
def format_addr(s):
name, addr = parseaddr(s)
return formataddr((Header(name, 'utf-8').encode(), addr.encode('utf-8') if isinstance(addr, unicode) else addr))
#info
user = raw_input("user's e-mail: ")
pwd = raw_input("user's password: ")
customer = raw_input("customer's e-mail: ")
sever = raw_input("SMTP sever: ")
#instance
lis1 = user.split('@')
na1 = lis1[0]
lis2 = customer.split('@')
na2 = lis2[0]
msg = MIMEMultipart()
msg["Subject"] = Header("test 附件", "utf-8").encode()
msg["From"] = format_addr(u"%s <%s>" %(na1, user))
msg["To"] = format_addr(u"%s <%s>" %(na2, customer))
#text
txt = MIMEText("附件", 'plain', 'utf-8')
msg.attach(txt)
#pdfs as attachment
atta = MIMEApplication(open('uv_ecmwf.pdf', 'rb').read())
atta.add_header('Content-Disposition', 'attachment', filename='uv_ecmwf.pdf')
msg.attach(atta)
atta = MIMEApplication(open('wrf_UV.pdf', 'rb').read())
atta.add_header('Content-Disposition', 'attachment', filename='wrf_UV.pdf')
msg.attach(atta)
#transport
ser = smtplib.SMTP(sever, 25)
ser.set_debuglevel(1)
ser.login(user, pwd)
ser.sendmail(user, [customer], msg.as_string())
ser.quit()
#####################################################################################
#
#关于底图shape文件的使用，以中国为例(*代表CHN_adm1)
#下载shape文件得到：*.shp, *.shx, *.dbf, *.prj, *.cpg
#可以在上面的代码中加入 m.readshapefile('CHN_adm1', 'CHN_adm1')
#