Python | 绘制影像地图（Basemap）
发布
于
2022-04-29 09:12:19
Basemap绘图
- Basemap简介
- Basemap与地图投影
- Basemap绘图实例
文中代码与数据请点击https://pan.bnu.edu.cn/v/link/view/0cd746194a1e42858583e84ac7fc4e40直接下载，不需要转存。
Basemap简介
- Basemap是Matplotlib的扩展，是具有专业标准的地图绘制工具。
- Basemap适合地球科学家，特别是海洋学家和气象学家的需求。
-
最初编写Basemap是用来帮助和研究气候和天气预报的。
- 类似的工具还有CDAT，Python第一种用于绘制地图投影数据。CDAT后来进化为UVCDAT。
- NCL是UNIX/Linux系统中最常用的气象绘图工具，它在python中有PyNIO和PyNGL两个包。目前NCL已经停止支持，被GeoCAT代替，但python的两个包还在。
- Basemap作为轻量的绘图软件包，一直很受欢迎。多年来，Basemap 的功能随着各个学科(如生物学，地质学和地球物理学)的科学家的要求和贡献的新功能而演变。
- Basemap可以与 matplotlib 的一般绘图功能相结合，并在地图上绘制数据。
- Basemap本身不会进行任何绘图，但提供了坐标转换功能(using the PROJ.4 C library)。
- 在Basemap底层使用了GEOS库，用来将海岸线和边界特征剪切到所需的地图投影区域。
- 随着python2.7的寿终正寝，Basemap的支持也结束了，后续会被Cartopy替代。
- 然而Cartopy目前仍不完善，且中文资料匮乏，不便学习。Basemap资源丰富，易学易用。建议仍然以Basemap为主。
Basemap与地图投影
# 首先导入basemap和matplotlib两个包，两者都是必要的。
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
# 新建第一个地图
map = Basemap()
# 在使用 Basemap 类创建地图时具有许多选项。
# 在没有传递任何选项的 情况下，地图具有以经度 =0 和纬度 = 0 为中心的 Plate Carrée 投影(等距圆柱投影)。# 绘制海岸线
map.drawcoastlines()
# 如果使用单独的python程序（.py文件），需要下面这句话才能看到图
plt.show()<matplotlib.collections.LineCollection at 0x11541a978>
- 为了在二维地图上表示地球的曲面，则需要进行地图投影。
- 地图投影的方法有许多种，每种方法都有自己的优点和缺点。
- Basemap提供了34种地图投影方法。有些是全球性的，有些只能代表区域。
-
在创建Basemap类实例时，必须指定所需的地图投影。有两种方法
- 提供矩形映射投影区域的四个角的每一个的纬度和经度值。
- 提供地图投影区域中心的lat/lon值以及地图投影坐标中的 区域的宽度和高度。
- 类变量supported_projections是一个字符串，包含有关 Basemap支持的所有投影的信息。
- 投影需要在创建basemap的实例的时候就指定好。有的投影不需要参数，有的还需要用户指定参数。关于投影的详细信息请参考Basemap官方文档
# 看看Basemap支持哪些投影
import mpl_toolkits.basemap
print(mpl_toolkits.basemap.supported_projections)
# 使用help函数可查看Basemap构造函数的详细信息
# help(Basemap) cyl Cylindrical Equidistant
merc Mercator
tmerc Transverse Mercator
omerc Oblique Mercator
mill Miller Cylindrical
gall Gall Stereographic Cylindrical
cea Cylindrical Equal Area
lcc Lambert Conformal
laea Lambert Azimuthal Equal Area
nplaea North-Polar Lambert Azimuthal
splaea South-Polar Lambert Azimuthal
eqdc Equidistant Conic
aeqd Azimuthal Equidistant
npaeqd North-Polar Azimuthal Equidistant
spaeqd South-Polar Azimuthal Equidistant
aea Albers Equal Area
stere Stereographic
npstere North-Polar Stereographic
spstere South-Polar Stereographic
cass Cassini-Soldner
poly Polyconic
ortho Orthographic
geos Geostationary
nsper Near-Sided Perspective
sinu Sinusoidal
moll Mollweide
hammer Hammer
robin Robinson
kav7 Kavrayskiy VII
eck4 Eckert IV
vandg van der Grinten
mbtfpq McBryde-Thomas Flat-Polar Quartic
gnom Gnomonic
rotpole Rotated Pole
# 绘制海岸线+经纬度
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
map = Basemap(projection='cyl')
map.drawcoastlines()
map.drawmeridians(np.arange(0, 360, 30))
map.drawparallels(np.arange(-90, 90, 30))
plt.show()
# 换个投影试试，绘制南极洲，需要boudninglat和lon_0两个额外的参数
map = Basemap(boundinglat=-60, lon_0=0, projection='spstere')
map.drawcoastlines()
map.drawmeridians(np.arange(0, 360, 30))
map.drawparallels(np.arange(-90, 90, 30))
plt.show()
# 下面展示如何绘制Robinson投影
from mpl_toolkits.basemap import Basemap
import numpy as np
import matplotlib.pyplot as plt
# lon_0表示图的中央经线，一般设为0表示将格林尼治设为地图的中央 。
# 设为-180（注意负值）表示把太平洋放在中央，常用于绘制海表温度 。
# resolution='c' 表示海岸线精度为“粗糙”
m = Basemap(projection='robin',lon_0=-180,resolution='c')
m.drawcoastlines()
m.fillcontinents(color='coral',lake_color='aqua')
m.drawparallels(np.arange(-90.,120.,30.))
m.drawmeridians(np.arange(0.,360.,60.))
m.drawmapboundary(fill_color='aqua')
# 添加标题的方法和matplotlib一样
plt.title("Robinson Projection")
plt.show()
# 下面展示如何进行地图的坐标变换
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
# 兰伯特正轴等角割圆锥投影，设定投影参数
m = Basemap(width=12000000,height=9000000,projection='lcc',
resolution='c',lat_1=45.,lat_2=55,lat_0=50,lon_0=-107.)
# 设定海洋的颜色
m.drawmapboundary(fill_color='aqua')
# 设定陆地和湖泊的颜色
m.fillcontinents(color='coral',lake_color='aqua')
# 设定经纬度坐标，图的左侧和右侧显示维度，上侧和下侧显示经度
parallels = np.arange(0.,81,10.)
m.drawparallels(parallels,labels=[False,True,True,False])
meridians = np.arange(10.,351.,20.)
m.drawmeridians(meridians,labels=[True,False,False,True])
# 对Boulder的经纬度进行坐标变换，将地理坐标系转变为投影坐标系
lon, lat = -104.237, 40.125 # Location of Boulder
# 通过basemap对象，将地理坐标直接变换为投影坐标。basemap对象m的实质是一个坐标转换器。
xpt,ypt = m(lon,lat)
# 设定inverse参数为True，可以把投影坐标转换为地理坐标，再转回来。
lonpt, latpt = m(xpt,ypt,inverse=True)
# 使用m.plot函数，将转化后的坐标绘制在地图上
m.plot(xpt,ypt,'bo') # plot a blue dot there
# 在指定的位置上标注文字
plt.text(xpt+100000,ypt+100000,'Boulder (%5.1fW,%3.1fN)' % (lonpt,latpt))
plt.show()
使用EPSG设置投影
EPSG代码是使用数字代码来命名投影的标准的。Basemap允许使用此表示法来创建地图，但仅在某些情况下可以。若要使用它，需要将epsg参数 传递给带有代码的Basemap构造函数。Basemap支持的epsg代码在文件
/mpl_toolkit/basemap/data/epsg
中。即使所需的epsg出现 在文件中，但有时库并不能使用投影。
# 使用EPSG代码的例子
# 使用UTM 31N绘制西班牙梅诺卡岛
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
map = Basemap(llcrnrlon=3.75,llcrnrlat=39.75,urcrnrlon=4.35,urcrnrlat=40.15, resolution = 'h', epsg=5520)
map.drawmapboundary(fill_color='aqua')
map.fillcontinents(color='coral',lake_color='aqua')
map.drawcoastlines()
plt.show()
Basemap绘图实例：
下面为大家介绍三个绘图实战案例，帮助大家绘制发表论文所需的图表，达到出版水平。主要介绍三个例子：
- 全球地震数据绘图：单点矢量数据
- 格陵兰冰盖数据绘图：坐标变换与colorbar
- 陆面过程数据处理：nc文件读取，用shp文件裁剪数据
USGS全球地震数据集
- 美国政府维护着一系列来自最近地震事件的地震实时数据的接口。
- 可以选择获取过去一小时到过去三十天的数据，可以选择检查具有各种大小的事件的数据。
- 对于这个案例项目，将使用包含过去七天内的所有地震事件的数据集，其具有1.0或更大的量值。
- 这个数据接口有多种形式。在第一个示例中，将以csv格式(逗号分隔值)解析文件。
- 如果查看数据集的文本文件的前几行，可以识别出最相关的信息。
# 读取csv文件，打印前3行
with open('all_week.csv','r') as f:
lines = f.readlines()
for line in lines[:3]:
print(line)time,latitude,longitude,depth,mag,magType,nst,gap,dmin,rms,net,id,updated,place,type,horizontalError,depthError,magError,magNst,status,locationSource,magSource
2017-03-10T15:07:12.330Z,34.154,-117.4248333,12.03,1.93,ml,10,161,0.2114,0.31,ci,ci37824064,2017-03-10T15:10:59.882Z,"7km NNE of Fontana, CA",earthquake,1.37,7.46,0.271,34,automatic,ci,ci
2017-03-10T15:06:28.638Z,62.9005,-149.6816,79.1,1.6,ml,,,,0.88,ak,ak15458857,2017-03-10T15:08:59.369Z,"65km SW of Cantwell, Alaska",earthquake,,0.8,,,automatic,ak,ak
# 目前我们只需要提取地震的纬度和经度。看第一行，我们只需要第二列、第三列、第五列。
# 使用Python的csv模块模块处理数据，这简化了使用csv文件的过程。
# 以下代码生成两个列表，包含文件中每个地震的纬度和经度:
import os, csv
fig_index = 0
filename = 'all_week.csv'
lats, lons, magnitudes = [], [], []
with open(filename) as f:
reader = csv.reader(f)
next(reader)
for row in reader:
lats.append(float(row[1]))
lons.append(float(row[2]))
magnitudes.append(float(row[4]))# 绘制地震的空间位置
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
eq_map = Basemap(projection='robin', resolution='l',
area_thresh=1000.0, lat_0=0, lon_0=-130)
eq_map.drawcoastlines()
eq_map.drawmapboundary()
eq_map.fillcontinents(color='grey')
# 将经纬度转换为投影坐标
x,y = eq_map(lons, lats)
eq_map.plot(x,y,'ro',markersize=6)
plt.show()
# 绘制地震的震级
# 刚刚绘制的地图有个缺陷，大地震和小地震都是一样大小，一样颜色的点，并没有区分。
# 下面使用循环遍历的方式绘制每个点，根据震级调整每个点的大小。
eq_map.drawcoastlines()
eq_map.fillcontinents(color='grey')
min_maker_size = 1.2
# lon,lat,mag分别代表经度、纬度、震级
for lon,lat,mag in zip(lons,lats,magnitudes):
x,y = eq_map(lon,lat)
msize = mag*min_maker_size
eq_map.plot(x,y,'ro',markersize=msize)
plt.show()
# 使用不同的颜色来表示震级，小地震为绿色，中等地震为黄色，大地震为红色
# 首先按照震级设定颜色
def get_maker_color(mag):
if mag < 3.0:
return 'go'
elif mag < 5.0:
return 'yo'
else:
return 'ro'
eq_map.drawcoastlines()
eq_map.fillcontinents(color='grey')
min_maker_size = 1.2
# lon,lat,mag分别代表经度、纬度、震级
for lon,lat,mag in zip(lons,lats,magnitudes):
x,y = eq_map(lon,lat)
msize = mag*min_maker_size
mstr = get_maker_color(mag)
eq_map.plot(x,y,mstr,markersize=msize)
plt.show()
格陵兰冰盖表面质量平衡
- 表面质量平衡（SMB，surface mass balance），表示冰盖表面的质量累积/亏损。
- SMB是陆面过程模式的输出，也是冰盖模式的上边界条件。
- 这个例子比较综合，把矢量数据、栅格数据、colorbar、numpy保存文件等内容都集成在一起了。
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import matplotlib.colors as colors
# 载入数据，事先裁剪好的表面物质平衡（surface mass balance）数据，按经纬度坐标存储为numpy数组
iyr,imm = 1850,7
npzf = np.load('Greenland_%04d-%02d.npz' % (iyr, imm))
lat = npzf['lat']
lon = npzf['lon']
smb = npzf['smb']
# 设定掩膜，排除不合理的值
mask = (smb > 999)
smb = np.ma.array(smb, mask = mask)*1e6
# 设定数据上下限
vmin,vmax = -150,150
# 新建绘图，指定尺寸
fig = plt.figure(figsize=(6,10))
ax = fig.add_subplot(111)
# 绘制北极，注意设定的投影参数
m = Basemap(width=2000000, height=3000000,
resolution = 'l', projection = 'stere',\
lat_ts = 72.5, lat_0 = 72.5, lon_0 = -45, ax =ax)
# 设定标题
title_name = 'Greenland_SMB_%04d-%02d' % (iyr, imm)
ax.set_title(title_name)
# 用一个好看的地形图作为底图
m.etopo(scale=1.5)
# 绘制经纬度线
m.drawmeridians(np.arange(-75, -15, 5),\
linewidth=0.5, fontsize=10, labels=[0,0,0,1],color='silver')
m.drawparallels(np.arange(60, 85, 5),rotation=90,\
linewidth=0.5, fontsize=10, labels=[1,0,0,0],color='silver')
# 绘制格陵兰边界
m.readshapefile('./Greenland_Boundary/GRL_adm0', 'GRL_adm0')
# 对栅格数据，使用meshgrid配合basemap对象进行坐标变换
mlon, mlat = np.meshgrid(lon, lat)
x, y = m(mlon, mlat)
# 设定colormap，把原来的jet转个方向
my_cmap = colors.LinearSegmentedColormap.reversed(plt.cm.jet)
# 设定colorbar的标题
cb_label = 'SMB: Surface Mass Balance (10$^{-6}$mm/s)'
# 绘图，使用指定的colormap
cm = m.pcolormesh(x, y, smb, vmin = vmin, vmax = vmax, cmap = my_cmap)
# 设定colorbar作为图例
cb = m.colorbar(cm, location='bottom',extend='max',size='3%',pad='10%')
cb.set_label(label = cb_label ,size=12)#weight='bold'
cb.ax.tick_params(labelsize=10)
plt.show()
fig.savefig(fr'{title_name}.png', format='png', dpi=300) #保存图片
plt.close()
陆面过程模式数据处理
- 绘制陆面过程模式CLM的输出结果，这里只绘制了降水。
- 这个例子更加综合，还涉及到数据裁剪与netcdf文件读取。
- 建议使用panoply对netcdf文件进行快捷查看，很方便。
- 如果你有linux系统，并且安装好了netcdf的库，可以用ncdump将netcdf文件整个转换为文本文件，可以用ncdump -h快速查看文件头。在anaconda的环境中安装了netcdf的包以后也可以使用ncdump。
# 绘制全球降水的例子
import os, shutil
import numpy as np
import netCDF4 as nc
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import matplotlib.colors as colors
# 加载CLM的surfdata文件，包含经纬度坐标、海陆边界
nc_surfdata_ds = nc.Dataset('../clm_data/surdata_0.9x1.25_simyr2000_c110921.nc', mode='r')
lat = nc_surfdata_ds.variables['LATIXY'][:,0]
lon = nc_surfdata_ds.variables['LONGXY'][0,:]
nlat = lat.shape[0]
nlon = lon.shape[0]
# CLM文件中经度按照0-360度存储，绘图的时候要改成[-180,180]
half = int(nlon/2)
lonidx = [i for i in range(half,nlon)] + [i for i in range(half)]
lon = lon[lonidx]
lon[:half] -= 360
# 提取CLM中的海陆边界，作为数组的mask
landfrac = nc_surfdata_ds.variables['LANDFRAC_PFT'][:]
landfrac = landfrac[:,lonidx]
mask = (landfrac < 0.5)
# 读取nc数据
iyr, imm = 1850, 7
nc_ds = nc.Dataset('../clm_data/BGCN_PI_1deg.clm2.h0.%04d-%02d.nc' % (iyr, imm), mode = 'r')
rain = nc_ds.variables['RAIN'][0,:,:]*1e6
rain = rain[:,lonidx]
# 设定数据上下限
vmin,vmax = 0,100
# 下面准备绘图
fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111)
# 将数据加上掩膜
data = np.ma.array(rain, mask = mask, fill_value =np.nan).filled()
# 建立basemap对象，设定坐标系
m = Basemap(projection = 'cyl', resolution = 'i', ax = ax)
# 设定标题
title_name = 'CLM_RAIN_%04d-%02d' % (iyr, imm)
ax.set_title(title_name)
# 设定经纬度
m.etopo(scale=1.5)
m.drawmeridians(np.arange(-180, 180, 30),\
linewidth=0.5, fontsize=10, labels=[0,0,0,1],color='silver')
m.drawparallels(np.arange(-90, 90, 30),rotation=90,\
linewidth=0.5, fontsize=10, labels=[1,0,0,0],color='silver')
# 设定colormap，把原来的jet转个方向
my_cmap = colors.LinearSegmentedColormap.reversed(plt.cm.jet)
# 设定colorbar的标题
cb_label = 'RAIN: Atmospheric Rain (10$^{-6}$mm/s)'
# 不裁剪
cm = m.pcolormesh(lon, lat, data, vmin = vmin, vmax = vmax,cmap = my_cmap)
cb = m.colorbar(cm, location='bottom',extend='max',size='3%',pad='10%')
cb.set_label(label = cb_label ,size=12)# weight='bold'
cb.ax.tick_params(labelsize=10)
plt.show()
fig.savefig(fr'{title_name}_global.png', format='png', dpi=300) #保存图片
plt.close()
# 只绘制格陵兰岛上的降水，演示如何使用shp文件进行裁剪
import os, shutil
import numpy as np
import netCDF4 as nc
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.path import Path
from matplotlib.patches import PathPatch
import matplotlib.colors as colors
# 加载CLM的surfdata文件，包含经纬度坐标、海陆边界
nc_surfdata_ds = nc.Dataset('../clm_data/surfdata_0.9x1.25_simyr2000_c110921.nc', mode='r')
lat = nc_surfdata_ds.variables['LATIXY'][:,0]
lon = nc_surfdata_ds.variables['LONGXY'][0,:]
nlat = lat.shape[0]
nlon = lon.shape[0]
# CLM文件中经度按照0-360度存储，绘图的时候要改成[-180,180]
half = int(nlon/2)
lonidx = [i for i in range(half,nlon)] + [i for i in range(half)]
lon = lon[lonidx]
lon[:half] -= 360
# 提取CLM中的海陆边界，作为数组的mask
landfrac = nc_surfdata_ds.variables['LANDFRAC_PFT'][:]
landfrac = landfrac[:,lonidx]
mask = (landfrac < 0.5)
# 读取nc数据
iyr, imm = 1850, 7
nc_ds = nc.Dataset('../clm_data/BGCN_PI_1deg.clm2.h0.%04d-%02d.nc' % (iyr, imm), mode = 'r')
rain = nc_ds.variables['RAIN'][0,:,:]*1e6
rain = rain[:,lonidx]
# 设定数据上下限
vmin,vmax = 0,5
# 下面准备绘图
fig = plt.figure(figsize=(10,6))
ax = fig.add_subplot(111)
# 将数据加上掩膜
data = np.ma.array(rain, mask = mask, fill_value =np.nan).filled()
# 建立basemap对象，设定坐标系
m = Basemap(projection = 'cyl', resolution = 'i',\
llcrnrlon = -75, llcrnrlat = 60,\
urcrnrlon = -15, urcrnrlat = 85, ax = ax)
# 设定标题
title_name = 'CLM_RAIN_%04d-%02d' % (iyr, imm)
ax.set_title(title_name)
# 设定经纬度
m.etopo(scale=1.5)
m.drawmeridians(np.arange(-75, -15, 5),\
linewidth=0.5, fontsize=10, labels=[0,0,0,1],color='silver')
m.drawparallels(np.arange(60, 85, 5),rotation=90,\
linewidth=0.5, fontsize=10, labels=[1,0,0,0],color='silver')
# draw Greenland boundary
m.readshapefile('./Greenland_Boundary/GRL_adm0', 'GRL_adm0')
# clip
import shapefile
sf = shapefile.Reader("./Greenland_Boundary/GRL_adm0")
for shape_rec in sf.shapeRecords():
vertices = []
codes = []
pts = shape_rec.shape.points
prt = list(shape_rec.shape.parts) + [len(pts)]
for i in range(len(prt) - 1):
for j in range(prt[i], prt[i+1]):
vertices.append((pts[j][0], pts[j][1]))
codes += [Path.MOVETO]
codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)
codes += [Path.CLOSEPOLY]
clip = Path(vertices, codes)
clip = PathPatch(clip, transform=ax.transData)
# 设定colormap，把原来的jet转个方向
my_cmap = colors.LinearSegmentedColormap.reversed(plt.cm.jet)
# 设定colorbar的标题
cb_label = 'RAIN: Atmospheric Rain (10$^{-6}$mm/s)'
# 按格陵兰的范围裁剪
cm = m.pcolormesh(lon, lat, data, vmin = vmin, vmax = vmax, cmap = my_cmap, clip_path = clip)
cb = m.colorbar(cm, location='bottom',extend='max',size='3%',pad='10%')