利用python按照不规则区域处理 利用RGB绘制填色图

xunmengfei8

# coding=utf-8
from typing import Any

import netCDF4 as nc
import numpy as np
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import cartopy.mpl.ticker as cticker
from numpy import random, ndarray, dtype, floating
from cartopy.io import shapereader
from wrf import to_np, getvar, interplevel, smooth2d, get_cartopy, cartopy_xlim, cartopy_ylim, latlon_coords, vertcross, smooth2d, CoordPair, GeoBounds,interpline
from matplotlib import colors
from matplotlib.colors import LinearSegmentedColormap
import xarray as xr
import matplotlib.colors as col
import cmaps

def rgb_to_hex(rgb):
    return '#{:02x}{:02x}{:02x}'.format(rgb[0], rgb[1], rgb[2])
from matplotlib.ticker import FuncFormatter
def format_func(value, tick_number):
        # 自定义格式化函数
        if value < 1:
                return f'{value:.1f}'  # 小于0.1时使用科学计数法
        else:
                return f'{value:.0f}'  # 其他情况保留两位小数
# cbar = plt.colorbar(format=FuncFormatter(format_func))
# NCL precip3_16lev
NCL1 = rgb_to_hex((254, 255, 255))
NCL2 = rgb_to_hex((255, 245, 204))
NCL3 = rgb_to_hex((255, 230, 112))
NCL4 = rgb_to_hex((255, 204, 51))
NCL5 = rgb_to_hex((255, 175, 51))
NCL6 = rgb_to_hex((255, 153, 51))
NCL7 = rgb_to_hex((255, 111, 51))
NCL8 = rgb_to_hex((255, 85, 0))
NCL9 = rgb_to_hex((230, 40, 30))
NCL10 = rgb_to_hex((200, 20, 30))
NCL11 = rgb_to_hex((166, 15, 20))
NCL12 = rgb_to_hex((120, 10, 15))
NCL13 = rgb_to_hex((120, 0, 90))
NCL14 = rgb_to_hex((148, 0, 211))
NCL15 = rgb_to_hex((186, 85, 211))
NCL16 = rgb_to_hex((160, 160, 160))
NCL17 = rgb_to_hex((192, 192, 192))
color = ( NCL1, NCL2, NCL3, NCL4, NCL5, NCL6, NCL7, NCL8, NCL9,
              NCL10, NCL11, NCL12, NCL13, NCL14, NCL15,NCL16,NCL17)
# cpool = [ NCL1, NCL2, NCL3, NCL4, NCL5, NCL6, NCL7, NCL8, NCL9,
#              NCL10, NCL11, NCL12, NCL13, NCL14, NCL15,NCL16,NCL17]
# cmap = col.ListedColormap(cpool, 'indexed')
# print(cpool)  # 输出: #ffa500
#mycolor=['white', 'lightgreen' ,'cyan','lightyellow','yellow','orange' , 'red']
#cmap = colors.LinearSegmentedColormap.from_list('my_list', mycolor)

shp_path = r'/mnt/d/LPCC/ncl/chinamap/Chinap_Areas.shp'
shpn_path = r'/mnt/d/LPCC/ncl/chinamap/nine.shp'
wrf_file = '/mnt/d/BaiduNetdiskDownload/wrfoutnco985.nc'

wrf_data = nc.Dataset(wrf_file, 'r')
#print(wrf_data.variables)
lat = wrf_data.variables['XLAT'][0, :, :]
lon = wrf_data.variables['XLONG'][0, :, :]
print(lat,lon.shape)

for i in range(0,145):#左闭右开
        dtim = wrf_data.variables['Times'][i,:]
        dstr =  ''.join([p.decode() for p in dtim])
        stim = dstr[0:4] + dstr[5:7] + dstr[8:10] + dstr[11:13] + '00'
        print(stim)
        plt.close()
        ua = wrf_data.variables['U'][i,:,:]
        va = wrf_data.variables['V'][i,:,:]
        p = wrf_data.variables['PB'][i,:,:]/100
        z = wrf_data.variables['Z'][i,:,:]
        u = ua[:,:,0:479]
        v = va[:,0:319,:]
        pm1 = wrf_data.variables['DUST_1'][i,:,:]
        pm2 = wrf_data.variables['DUST_2'][i, :, :]
        pm3 = wrf_data.variables['DUST_3'][i, :, :]
        pm4 = wrf_data.variables['DUST_4'][i, :, :]
        pm5 = wrf_data.variables['DUST_5'][i, :, :]
        pm = pm1+pm2+pm3+pm4+pm5
        #print(p.shape,z.shape,u.shape,v.shape)
        pp=700
        h_p = interplevel(z, p, pp)
        u_p = interplevel(u, p, pp)
        v_p = interplevel(v, p, pp)
        w_p = np.sqrt(u_p**2 + v_p**2)
        pm_p = interplevel(pm, p, pp)
        dsat = pm_p
        # print(h_p.shape,u_p.shape,v_p.shape,pm_p.shape,w_p.shape)

        proj = ccrs.PlateCarree(central_longitude=245)     #选择投影
        fig = plt.figure(figsize=(6, 4.9), dpi=200)  # 创建图形
        fig_ax1 = fig.add_axes([0.1, 0.1, 0.87, 0.87], projection=proj)
        proj = ccrs.PlateCarree(central_longitude=90)  # 设置边界
        leftlon, rightlon, lowerlat, upperlat = (70, 141, 15, 61)
        img_extent = [leftlon, rightlon, lowerlat, upperlat]
        fig_ax1.set_extent(img_extent, crs=ccrs.PlateCarree())
        # 绘制海岸线和湖泊等地理特征
        reader = shapereader.Reader(shp_path)
        for record in reader.records():
                fig_ax1.add_geometries([record.geometry], ccrs.PlateCarree(), facecolor='none')
        reader = shapereader.Reader(shpn_path)
        for record in reader.records():
                fig_ax1.add_geometries([record.geometry], ccrs.PlateCarree(), facecolor='none')
        #设置刻度及刻度标签格式
        fig_ax1.set_xticks(np.arange(leftlon,rightlon,10), crs=ccrs.PlateCarree())
        fig_ax1.set_yticks(np.arange(lowerlat,upperlat,5), crs=ccrs.PlateCarree())
        lon_formatter = cticker.LongitudeFormatter()
        lat_formatter = cticker.LatitudeFormatter()
        fig_ax1.xaxis.set_major_formatter(lon_formatter)
        fig_ax1.yaxis.set_major_formatter(lat_formatter)
        #绘制相关系数填色
        #print(r)
        #根据数据序列调整色标levels =np.arange(-.05,.05,0.005)
        # print(lon.shape,lat.shape,dsat.shape)
        #clevs = [0, 0.1, 0.2, 0.3, 0.5, 0.7, 1, 1.5, 2, 2.5, 3, 3.5, 4, 5, 6, 7, 8]  # 200
        #clevs = [50,60,80,100,120,150,180, 210, 250,290, 330, 370, 410,450,500,550,600] #400
        #clevs = [50,60,80,100,120,150,180, 210, 250,290, 330, 370, 410,450,500,550,600] #500
        clevs = [50,60,80,100,130,160,190, 220, 250,300, 350, 400, 450,500,550,600,700] #700
        #clevs = [50,60,80,100,130,160,190, 220, 250,300, 350, 400, 450,500,550,600,700] #850
        cf1 = fig_ax1.contourf(lon, lat, dsat, zorder=0, levels=clevs, extend='both', transform=ccrs.PlateCarree(),
                               colors=color)  # 绘制显著性打点。思路为将0-0.05范围内的区域用点的标记来填充，来表示显著性95%水平。
        # 色标
        position = fig.add_axes([0.2, 0.1, 0.6, 0.025])
        cbar = fig.colorbar(cf1, cax=position, orientation='horizontal', format=FuncFormatter(format_func))
        cbar.ax.xaxis.set_tick_params(labelrotation=90)
        cbar.set_ticks(clevs)  # 设置刻度
        # plt.show()
        lonc = to_np(lon[::20,::20])
        latc = to_np(lat[::20,::20])
        uc = to_np(u_p[::20,::20])
        vc = to_np(v_p[::20,::20])
        hm = to_np(h_p[::20,::20])/9.8
        hp = smooth2d(smooth2d(hm, 3, cenweight=4), 3, cenweight=4)
        # print(lonc.dtype,latc.dtype,uc.dtype,vc.dtype)
        m = fig_ax1.contour(lonc, latc, hp, linewidths=1.3, colors='b', levels=np.arange(100, 1600,4),transform=ccrs.PlateCarree())  # 200_700
        #m = fig_ax1.contour(lonc,latc, hp, linewidths=1.3, colors='b',  levels=np.arange(200,1600,4),transform=ccrs.PlateCarree()) #200_700
        #m = fig_ax1.contour(lonc, latc, hp, linewidths=1.3, colors='b', levels=np.arange(100, 240, 2),transform=ccrs.PlateCarree())  #850
        fig_ax1.clabel(m, inline=True, fmt='%.f', fontsize=10)

        # c = fig_ax1.quiver(lonc,latc, uc, vc, color='b', angles='xy', scale=1000, width=0.002, transform=ccrs.PlateCarree())
        # plt.quiverkey(c, X=0.98, Y=1.02, U=20, angle=90, labelpos='W', label='Wind:20m/s', color='r', labelcolor='r')
        # c = fig_ax1.streamplot(lon, lat, u, v, color = 'navy', density = 3.0, linewidth = 1.2, transform = ccrs.PlateCarree())

        s = fig_ax1.barbs(lonc,latc, uc, vc, linewidth=0.4,
                          barb_increments={'half': 2, 'full': 4, 'flag': 20},
                          sizes=dict(spacing=0.13, width=0.18, emptybarb=0.26, height=0.3),
                          length=5.4, zorder=9,transform=ccrs.PlateCarree())

        plt.title("Dust(ug/kg) & UV (m/s) & High(dgpm)    "+stim, x=0.5, y=31)

        # plt.show()
        plt.savefig('/mnt/d/LPCC/ncl/WRF_LES/HIGH/COMBIN/7001/'+stim+'.png')
        #plt.close()