小白学习Basemap气象画地图的第三天(中国温度分布图,mask外部) |
您所在的位置:网站首页 › 中国地图空白填充图 › 小白学习Basemap气象画地图的第三天(中国温度分布图,mask外部) |
小白学习Basemap气象画地图的第三天(中国温度分布图,mask外部)
首先还是感谢公众号(气象学家),代码和测试数据来自与他,不过这次有长进了,自己学会修改了。还是逐条向大家解释。 (和大家分享一个经验,在代码中查看某个函数或者变量的定义,可以利用快捷键快速寻找,pycharm里面是ctrl + 鼠标左键点击,这样的意义在于可以定位到最原始的构造函数,定义等位置,从而了解其用法,参数设置。要学会看__init__()里面的东西) 这次是用Basemap库画的,这个库在线安装好像已经停止了,只能下载离线包进行安装了,大家注意一下,首先效果图奉上: 导入库 import cartopy.crs as ccrs import matplotlib.pyplot as plt from matplotlib.path import Path from matplotlib.patches import PathPatch import numpy as np import shapefile import xarray as xr from mpl_toolkits.basemap import Basemap 绘图区设置 #设置画图字体的大小 plt.rcParams.update({'font.size':20}) #解决中文乱码问题 plt.rcParams['font.sans-serif'] = ['SimHei'] #解决负号乱码问题 plt.rcParams['axes.unicode_minus'] = False #设置画布和绘图区 fig = plt.figure(figsize=[12,18]) ax = fig.add_subplot(111) 读取地图,设置mask路径(重点)这里是重点部分,小白容易看不懂,注释比较多,大家不要耐烦,文末作者会给出注释少的代码版本 #读取shp格式的地图,有三个文件分别为.dbf,.shp,.shx缺一不可 #这里的shp文件其实就是一堆点,把各个国家圈出来 #sf得到的是一个类,存放了全球的地理信息 sf = shapefile.Reader("country1") #sf.shapeRecords()里面放了各个国家地区的信息 #循环的目的是单个拿出来,找到中国 shape_rec是单个 for shape_rec in sf.shapeRecords(): #shape_rec.record内部存放了单一国家的很多信息,比如名字,货币等等 #其中shape_rec.record[2]放的是国家名 #可以print(shape_rec.record)看看 if shape_rec.record[2] == 'China': codes = []#用来存放移动路径(画图动作) #shape_rec.shape是一个类,图形类 #里面三个属性shapeType:图形类型 points图形边界坐标 parts图形起始索引 #解释一下parts属性,一个国家的边界可能不是全连在一起的,会分为一块一块,那么就相当于多个图形,存在shp文件内就不连续,parts里面就放了每个区域的起始索引(下标) pts = shape_rec.shape.points #上文已经说过了parts的意义,设想一下,两块区域的的起始坐标中间夹的不就是一块区域的所有坐标吗,但是最后一块区域没有结束值,所有加了[len(pts)],这就是最后一个点的索引。 prt = list(shape_rec.shape.parts) + [len(pts)] #下面的循环主要作用是:建立一个绘图路径,利用区块起始点的索引生成一个画图动作 for i in range(len(prt) - 1): codes += [Path.MOVETO]#点移动 codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)#画线 codes += [Path.CLOSEPOLY]#这块画完,循环结束,下一块 clip = Path(pts, codes)#利用数据和路径生成一个画图动作 clip = PathPatch(clip, transform=ax.transData)#再加入ax的变换大致原理是读取地图文件,大家俗称的shp文件,然后利用里面的点,得到一个或者多个画图路径(一般都是多个,国家边界很容易是很多段的),利用这个路径来把路径外部区域全部画成白的。细节可以看代码块的注释。 读取NC数据这一部分比较简单 #这里是读取NC数据部分,上一个帖子已经介绍了,不再赘述 ds = xr.open_dataset('EC-Interim_monthly_2018.nc') lat = ds.latitude lon = ds.longitude data = (ds['t2m'][0,::-1,:] - 273.15) # 把温度转换为℃ [0,::-1,:]表示第一个时次、纬度反向 画温度分部图这一部分和之前也大同小异 cbar_kwargs = { 'orientation': 'horizontal', 'label': 'Temperature (℃)', 'shrink': 0.02, 'ticks': np.arange(-25, 25 + 1, 5), 'pad': -0.28, 'shrink': 0.95 } levels = np.arange(-25, 25 + 1, 1) cs = data.plot.contourf(ax=ax,levels=levels,cbar_kwargs=cbar_kwargs, cmap='Spectral_r')到这里,咱们画出来的全球的温度分补(别忘了plt.show()) 利用掩膜路径画中国区域其实就一条代码,简单的很,拿来用就好了 #添加掩膜路径,白化外部的分部 for contour in cs.collections: contour.set_clip_path(clip)但是结果出现的问题,因为坐标变换的问题,生成的结果为 利用Basemap类调整图形 #生成一个Basemap类,用来变换坐标,画出合适尺寸的图 m = Basemap(llcrnrlon=72.0,#地图左边经度 llcrnrlat=10.0,#地图下方纬度 urcrnrlon=137.0,#地图右边经度 urcrnrlat=55.0,#地图上方纬度 resolution = None,#分辨率,这里设置为无 projection = 'cyl')#投影方式:默认,圆柱投影 #读取图形文件,画中国行政边界 m.readshapefile('bou2_4l','China Map',color='k',linewidth=1.2)readshapefile()函数的作用就是读取bou2_4l文件画图,结果为 图像修饰在进行一些修饰,让图片更好看 #画纬度 parallels = np.arange(10,55,10) m.drawparallels(parallels,labels=[True,True,True,True],color='dimgrey',dashes=[1, 3]) #画经度 meridians = np.arange(70,135,10) m.drawmeridians(meridians,labels=[True,True,False,True],color='dimgrey',dashes=[1, 3]) #移除原来的坐标轴标签 plt.ylabel('') plt.xlabel('') #设置标题 plt.rcParams.update({'font.size':30}) ax.set_title(u' 中国区域2018年1月平均气温',color='blue',fontsize= 25) #画南海 with open('CN-border-La.dat') as src: context = src.read() blocks = [cnt for cnt in context.split('>') if len(cnt) > 0] borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks] sub_ax = fig.add_axes([0.758, 0.249, 0.14, 0.155],projection=ccrs.LambertConformal(central_latitude=90, central_longitude=115)) for line in borders: sub_ax.plot(line[0::2], line[1::2], '-', lw=1, color='k',transform=ccrs.Geodetic()) sub_ax.set_extent([106, 127, 0, 25],crs=ccrs.PlateCarree()) 完整代码,注释简洁版最后有注释详细版,怕大家烦,就放最后了 import cartopy.crs as ccrs import matplotlib.pyplot as plt from matplotlib.path import Path from matplotlib.patches import PathPatch import numpy as np import shapefile import xarray as xr from mpl_toolkits.basemap import Basemap #设置画图字体的大小 plt.rcParams.update({'font.size':20}) #解决中文乱码问题 plt.rcParams['font.sans-serif'] = ['SimHei'] #解决负号乱码问题 plt.rcParams['axes.unicode_minus'] = False #设置画布和绘图区 fig = plt.figure(figsize=[12,18]) ax = fig.add_subplot(111) #读取shp格式的地图 sf = shapefile.Reader("country1") #得到mask白化路径 for shape_rec in sf.shapeRecords(): if shape_rec.record[2] == 'China': codes = [] pts = shape_rec.shape.points#边界点 prt = list(shape_rec.shape.parts) + [len(pts)]#区块起始索引 for i in range(len(prt) - 1): codes += [Path.MOVETO]#点移动 codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)#画线 codes += [Path.CLOSEPOLY]#这块画完,循环结束,下一块 clip = Path(pts, codes)#利用数据和路径生成一个画图动作 clip = PathPatch(clip, transform=ax.transData)#再加入ax的变换 #########################至此,所需国界图形读取完毕#################### #这里是读取NC数据部分,上一个帖子已经介绍了,不再赘述 ds = xr.open_dataset('EC-Interim_monthly_2018.nc') lat = ds.latitude lon = ds.longitude data = (ds['t2m'][0,::-1,:] - 273.15) # 把温度转换为℃ [0,::-1,:]表示第一个时次、纬度反向 ############################至此读取温度数据结束###################### #下面画 cbar_kwargs = { 'orientation': 'horizontal', 'label': 'Temperature (℃)', 'shrink': 0.02, 'ticks': np.arange(-25, 25 + 1, 5), 'pad': -0.28, 'shrink': 0.95 } levels = np.arange(-25, 25 + 1, 1) cs = data.plot.contourf(ax=ax,levels=levels,cbar_kwargs=cbar_kwargs, cmap='Spectral_r') #但是画出的全球温度分布 ############################至此温度画完了###################### #添加掩膜路径,白化外部的分部 for contour in cs.collections: contour.set_clip_path(clip) #生成一个Basemap类,用来变换坐标,画出合适尺寸的图 m = Basemap(llcrnrlon=72.0,#地图左边经度 llcrnrlat=10.0,#地图下方纬度 urcrnrlon=137.0,#地图右边经度 urcrnrlat=55.0,#地图上方纬度 resolution = None,#分辨率,这里设置为无 projection = 'cyl')#投影方式:默认,圆柱投影 #读取图形文件,画中国行政边界 m.readshapefile('bou2_4l','China Map',color='k',linewidth=1.2) ######################至此图基本画完了,下面是一些修饰过程###################### #画纬度 parallels = np.arange(10,55,10) m.drawparallels(parallels,labels=[True,True,True,True],color='dimgrey',dashes=[1, 3]) #画经度 meridians = np.arange(70,135,10) m.drawmeridians(meridians,labels=[True,True,False,True],color='dimgrey',dashes=[1, 3]) #移除原来的坐标轴标签 plt.ylabel('') plt.xlabel('') #设置标题 plt.rcParams.update({'font.size':30}) ax.set_title(u' 中国区域2018年1月平均气温',color='blue',fontsize= 25) #画南海 with open('CN-border-La.dat') as src: context = src.read() blocks = [cnt for cnt in context.split('>') if len(cnt) > 0] borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks] sub_ax = fig.add_axes([0.758, 0.249, 0.14, 0.155],projection=ccrs.LambertConformal(central_latitude=90, central_longitude=115)) for line in borders: sub_ax.plot(line[0::2], line[1::2], '-', lw=1, color='k',transform=ccrs.Geodetic()) sub_ax.set_extent([106, 127, 0, 25],crs=ccrs.PlateCarree()) ######################修饰完毕,出图###################### plt.savefig("China_mask_T2m.png", dpi=300, bbox_inches='tight') plt.show() 完整代码,注释详细版 import cartopy.crs as ccrs import matplotlib.pyplot as plt from matplotlib.path import Path from matplotlib.patches import PathPatch import numpy as np import shapefile import xarray as xr from mpl_toolkits.basemap import Basemap #设置画图字体的大小 plt.rcParams.update({'font.size':20}) #解决中文乱码问题 plt.rcParams['font.sans-serif'] = ['SimHei'] #解决负号乱码问题 plt.rcParams['axes.unicode_minus'] = False #设置画布和绘图区 fig = plt.figure(figsize=[12,18]) ax = fig.add_subplot(111) #读取shp格式的地图,有三个文件分别为.dbf,.shp,.shx缺一不可 #这里的shp文件其实就是一堆点,把各个国家圈出来 #sf得到的是一个类,存放了全球的地理信息 sf = shapefile.Reader("country1") #sf.shapeRecords()里面放了各个国家地区的信息 #循环的目的是单个拿出来,找到中国 shape_rec是单个 for shape_rec in sf.shapeRecords(): #shape_rec.record内部存放了单一国家的很多信息,比如名字,货币等等 #其中shape_rec.record[2]放的是国家名 #可以print(shape_rec.record)看看 if shape_rec.record[2] == 'China': codes = []#用来存放移动路径(画图动作) #shape_rec.shape是一个类,图形类 #里面三个属性shapeType:图形类型 points图形边界坐标 parts图形起始索引 #解释一下parts属性,一个国家的边界可能不是全连在一起的,会分为一块一块,那么就相当于多个图形,存在shp文件内就不连续,parts里面就放了每个区域的起始索引(下标) pts = shape_rec.shape.points #上文已经说过了parts的意义,设想一下,两块区域的的起始坐标中间夹的不就是一块区域的所有坐标吗,但是最后一块区域没有结束值,所有加了[len(pts)],这就是最后一个点的索引。 prt = list(shape_rec.shape.parts) + [len(pts)] #下面的循环主要作用是:建立一个绘图路径,利用区块起始点的索引生成一个画图动作 for i in range(len(prt) - 1): codes += [Path.MOVETO]#点移动 codes += [Path.LINETO] * (prt[i+1] - prt[i] -2)#画线 codes += [Path.CLOSEPOLY]#这块画完,循环结束,下一块 clip = Path(pts, codes)#利用数据和路径生成一个画图动作 clip = PathPatch(clip, transform=ax.transData)#再加入ax的变换 #########################至此,所需国界图形读取完毕#################### #这里是读取NC数据部分,上一个帖子已经介绍了,不再赘述 ds = xr.open_dataset('EC-Interim_monthly_2018.nc') lat = ds.latitude lon = ds.longitude data = (ds['t2m'][0,::-1,:] - 273.15) # 把温度转换为℃ [0,::-1,:]表示第一个时次、纬度反向 ############################至此读取温度数据结束###################### #下面画 cbar_kwargs = { 'orientation': 'horizontal', 'label': 'Temperature (℃)', 'shrink': 0.02, 'ticks': np.arange(-25, 25 + 1, 5), 'pad': -0.28, 'shrink': 0.95 } levels = np.arange(-25, 25 + 1, 1) cs = data.plot.contourf(ax=ax,levels=levels,cbar_kwargs=cbar_kwargs, cmap='Spectral_r') #但是画出的全球温度分布 ############################至此温度画完了###################### #添加掩膜路径,白化外部的分部 for contour in cs.collections: contour.set_clip_path(clip) #生成一个Basemap类,用来变换坐标,画出合适尺寸的图 m = Basemap(llcrnrlon=72.0,#地图左边经度 llcrnrlat=10.0,#地图下方纬度 urcrnrlon=137.0,#地图右边经度 urcrnrlat=55.0,#地图上方纬度 resolution = None,#分辨率,这里设置为无 projection = 'cyl')#投影方式:默认,圆柱投影 #读取图形文件,画中国行政边界 m.readshapefile('bou2_4l','China Map',color='k',linewidth=1.2) ######################至此图基本画完了,下面是一些修饰过程###################### #画纬度 parallels = np.arange(10,55,10) m.drawparallels(parallels,labels=[True,True,True,True],color='dimgrey',dashes=[1, 3]) #画经度 meridians = np.arange(70,135,10) m.drawmeridians(meridians,labels=[True,True,False,True],color='dimgrey',dashes=[1, 3]) #移除原来的坐标轴标签 plt.ylabel('') plt.xlabel('') #设置标题 plt.rcParams.update({'font.size':30}) ax.set_title(u' 中国区域2018年1月平均气温',color='blue',fontsize= 25) #画南海 with open('CN-border-La.dat') as src: context = src.read() blocks = [cnt for cnt in context.split('>') if len(cnt) > 0] borders = [np.fromstring(block, dtype=float, sep=' ') for block in blocks] sub_ax = fig.add_axes([0.758, 0.249, 0.14, 0.155],projection=ccrs.LambertConformal(central_latitude=90, central_longitude=115)) for line in borders: sub_ax.plot(line[0::2], line[1::2], '-', lw=1, color='k',transform=ccrs.Geodetic()) sub_ax.set_extent([106, 127, 0, 25],crs=ccrs.PlateCarree()) ######################修饰完毕,出图###################### plt.savefig("China_mask_T2m.png", dpi=300, bbox_inches='tight') plt.show()测试数据和文件需要的请下面留言,感谢支持 |
今日新闻 |
推荐新闻 |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |