Python | 绘制影像地图（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支持哪些投影 
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代码的例子 
# 使用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()

# 读取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()

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()

# 绘制全球降水的例子 
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%')
cb.set_label(label = cb_label ,size=12)#weight='bold' 
cb.ax.tick_params(labelsize=10)

plt.show()
fig.savefig(fr'{title_name}_Greenland.png', format='png', dpi=300) #保存图片 
plt.close()