python栅格数据出图可视化

import numpy as np
import cartopy.crs as ccrs
import cartopy.mpl.ticker as cticker
import matplotlib.pyplot as plt
import cartopy.feature as cfeature
from mpl_toolkits.basemap import Basemap
import cartopy.io.shapereader as shpreader
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from matplotlib import rcParams
from osgeo import gdal
import matplotlib.patches as mpatches
config = {"font.family":'Times New Roman', "font.size": 18, "mathtext.fontset":'stix'}
rcParams.update(config)
#------------------------------读取栅格------------------------
def Read_Tiif(fileName):
    #输入信息，文件名
    #输出信息，data_array,im_bands,im_geotrans,im_proj
    #分别是，矩阵数据，波段数，放射变换矩阵信息，投影信息
    dataset = gdal.Open(fileName)
    if dataset == None:
        print(fileName+"文件无法打开")
        return
    im_width = dataset.RasterXSize #栅格矩阵的列数
    im_height = dataset.RasterYSize #栅格矩阵的行数
    im_bands = dataset.RasterCount #波段数
    if im_bands < 1:
        print("影像没有可使用数据！")
        return
    im_data = dataset.ReadAsArray(0,0,im_width,im_height)#获取数据
    im_geotrans = dataset.GetGeoTransform()#获取仿射矩阵信息
    im_proj = dataset.GetProjection()#获取投影信息
    if im_bands == 1:
        data_array =  im_data[0:im_height,0:im_width]#获取第Band_num波段
    if im_bands > 1:
        data_array =  im_data[:,0:im_height,0:im_width]#获取第Band_num波段
    return data_array,im_width,im_height,im_bands,im_geotrans,im_proj

#-----------函数：添加指北针--------------
def add_north(ax, labelsize=16, loc_x=0.1, loc_y=0.85, width=0.04, height=0.13, pad=0.14):
    """
    画一个比例尺带'N'文字注释
    主要参数如下
    :param ax: 要画的坐标区域 Axes实例 plt.gca()获取即可
    :param labelsize: 显示'N'文字的大小
    :param loc_x: 以文字下部为中心的占整个ax横向比例
    :param loc_y: 以文字下部为中心的占整个ax纵向比例
    :param width: 指南针占ax比例宽度
    :param height: 指南针占ax比例高度
    :param pad: 文字符号占ax比例间隙
    :return: None
    """
    minx, maxx = ax.get_xlim()
    miny, maxy = ax.get_ylim()
    ylen = maxy - miny
    xlen = maxx - minx
    left = [minx + xlen*(loc_x - width*.5), miny + ylen*(loc_y - pad)]
    right = [minx + xlen*(loc_x + width*.5), miny + ylen*(loc_y - pad)]
    top = [minx + xlen*loc_x, miny + ylen*(loc_y - pad + height)]
    center = [minx + xlen*loc_x, left[1] + (top[1] - left[1])*.4]
    triangle = mpatches.Polygon([left, top, right, center], color='k')
    ax.text(s='N',
            x=minx + xlen*loc_x,
            y=miny + ylen*(loc_y - pad + height),
            fontsize=labelsize,
            horizontalalignment='center',
            verticalalignment='bottom')
    ax.add_patch(triangle)

## 数据读入
input_file = r'D:\exe\anaconda\file\make_map\PT_1981_08_05.tif'
data_array,im_width,im_height,im_bands,im_geotrans,im_proj = Read_Tiif(input_file)
left_lon,width_res,uperlat,heigh_res = im_geotrans[0], im_geotrans[1], im_geotrans[3],im_geotrans[5]
#lon = np.linspace(left_lon,24,25) * im_width
lon = np.arange(im_width) * width_res + left_lon
lat = np.arange(im_height) * heigh_res + uperlat
## 位势高度
data_array = data_array
data_array[np.where(data_array >1)] = np.nan
data_array[np.where(data_array <-1)] = np.nan
data_array = data_array.transpose(1,0)
#hgt = HGT[:,:,5,:]
HHGT = data_array
[LON,LAT] = np.meshgrid(lon,lat)
LON = LON.transpose(1,0)
LAT = LAT.transpose(1,0)
#print(U.variables.keys())
#创建画布
fig = plt.figure(figsize=(12,8))
ax=plt.axes(projection=ccrs.PlateCarree(central_longitude=85))
# 标注坐标轴
leftlon, rightlon, lowerlat, upperlat = (89,127,25,55)
img_extent = [leftlon, rightlon, lowerlat, upperlat]
ax.set_extent([leftlon, rightlon, lowerlat, upperlat], crs=ccrs.PlateCarree())
ax.set_xticks([ 90, 95, 100, 105, 110, 115, 120,125], crs=ccrs.PlateCarree())
ax.set_yticks([ 25, 30, 35,40,45,50,55], crs=ccrs.PlateCarree())
#zero_direction_label用来设置经度的0度加不加E和W
lon_formatter = LongitudeFormatter(zero_direction_label=False)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.set_title('Nei Meng SSM(v)',loc='center',fontsize =20)
# 加载中国地图
china = shpreader.Reader(r'D:\exe\anaconda\file\make_map\shp\shp\neimengg.dbf').geometries()
#绘制中国国界省界九段线等等
ax.add_geometries(china, ccrs.PlateCarree(),facecolor='none', edgecolor='black',zorder = 1)
# 绘制风场和位势图
cf=ax.contourf(LON[:],LAT[:],HHGT[:],100,cmap='gist_rainbow',transform=ccrs.PlateCarree(),extend='both')
#cf=ax.contourf(LON[:],LAT[:],HHGT[:],40,cmap='gist_rainbow',transform=ccrs.PlateCarree(),extend='both')#其中100是颜色的变化区间数目
#Q = ax.quiver(LON[:], LAT[:], UUwnd[:], VVwnd[:], color='black',width=0.0018,scale=100,headwidth=4,transform=ccrs.PlateCarree())     
# position=fig.add_axes([0.2, 0.1, 0.5, 0.03])#位置[左,下,右,上]
value_min = np.min(HHGT)
value_max = np.max(HHGT)
add_north(ax)
cbar=plt.colorbar(cf,orientation='vertical',shrink=0.92,aspect=20,fraction=.03,pad=0.01)   
#cbar=plt.colorbar(cf,orientation='vertical',shrink=0.92,aspect=20,fraction=.03,pad=0.01)   
#cbar = fig.colorbar(cf,orientation='horizontal',fraction=0.07, pad=0.05)
help(plt.colorbar)
plt.savefig(r'D:\exe\anaconda\file\make_map\plot33.png',dpi=800)
plt.show()

一些网站
Cartopy 绘图示例库
https://gnss.help/2018/04/24/cartopy-gallery/index.html