【pcolor数据可视化】Matlab vs. Python

1、Matlab代码及结果

代码

clear;clc
load('.\nclcolormap.mat')

sl = [0,50,100,200,500,0];
el = [50,100,200,500,1000,200];

for i = 1:length(sl)
    
    file = ['..\data\static_result\VIS_Min-',num2str(sl(i)),'to',num2str(el(i)),'_yearly.npy'];
    data = readNPY(file);
    mask=readNPY('.\mask.npy');
    data =data.*mask;
    siz = 30;
    
    lat = 15:0.05:55;
    lon = 70:0.05:140.05;
    [X,Y] = meshgrid(lon,lat);
    
    set(gcf,'color',[1 1 1],'position',[10 45 1000 800*1.2]);%get(0,'screensize')
    axes('position',[0.1 0.2 0.85 0.6]);
    
    % 数据映射: 数据分布差异较大，使用分位数将数据进行映射
    datax = data(:);
    datax(isnan(datax))=[];
    percentiles = prctile(datax(:), [0,50,80,90:2:98 99.9]);
    percentiles = unique(percentiles);
    x_ = linspace(percentiles(end-1), percentiles(end), 4);
    percentiles = [percentiles(1:end-2), x_];
    percentiles = unique(percentiles);
    xx1=percentiles(1:end-1);
    xx2=percentiles(2:end);
    CC1=0;CC2=length(xx1);
    data_map=data;
    for ii=1:length(xx1)
        data_map(data>=xx1(ii)&data<xx2(ii))=ii-0.5;
    end
    
    % 画图
    mypcolor(X,Y,data_map);shading interp
    hcb=colorbar;
    map = nclcolormap.WhBlGrYeRe;
    map_ = map(round(linspace(1, length(map), CC2)),:);
    colormap(map_)
    box on;hold on
    caxis([CC1 CC2]);
    set(hcb,'xtick',[CC1:CC2],'xticklabel',num2str(percentiles','%.1f'),'FontName','Times New Roman','fontsize',siz-15)
    set(hcb,'Location', 'southoutside','position',[0.1000 0.0903 0.8500 0.0278])
    set(get(hcb,'xlabel'),'string','\fontname{Aril}频次','fontsize',siz-10);
    % 地理信息
    path = '.\China\';
    China1=shaperead([path,'省.shp']);
    China2=shaperead([path,'九段线.shp']);
    plot([China1(:).X],[China1(:).Y],'color',[0.8 0.8 0.8],'linewidth',1); hold on
    plot([China2(:).X],[China2(:).Y],'color',[0.8 0.8 0.8],'linewidth',1); hold on
    axis([70 140 15 55]);hold on
    set(gca,'fontsize',siz-15,'fontname','Times New Roman',...
        'tickdir','out','ticklength',[0.01,0.05],'linewidth',1.2);
    xlabel('Lon(\circE)', 'fontsize',siz-12,'fontweight','bold')
    ylabel('Lat(\circN)', 'fontsize',siz-12,'fontweight','bold')
    title(['\fontname{Aril}2019-2023年平均',num2str(sl(i)),'-',num2str(el(i)),'m能见度频次'],'fontsize',siz-12,'fontweight','bold');
    % 小地图
    ax2=axes('position', [0.81 0.21 0.12 0.12]);
    mypcolor(X,Y,data_map);shading interp
    caxis([CC1 CC2]);
    colormap(map_)
    box on;hold on
    plot([China1(:).X],[China1(:).Y],'color',[0.7 0.7 0.7],'linewidth',0.5); hold on
    plot([China2(:).X],[China2(:).Y],'color',[0.7 0.7 0.7],'linewidth',0.5); hold on
    axis([104.5, 124, 0, 26]);
    set(gca,'xtick','','ytick','','layer','top');
    
    save_name=['../map/matlab_map/',num2str(sl(i)),'-',num2str(el(i)),'.png'];
    export_fig(save_name,'-r200');
    clf
end
close all

结果

2、Python代码及结果

代码

import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import matplotlib as mpl

import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
from matplotlib.colors import LinearSegmentedColormap

import regionmask
import geopandas as gpd
import warnings

warnings.filterwarnings("ignore")
import matplotlib.ticker as ticker


def map_frequency(condition = [0, 50], season = None):
    # the mask of Chinese Region: 中国区域的mask
    file = './china2.shp'
    countries = gpd.read_file(file)
    deg = 0.05
    lat = np.arange(15, 55 + deg, deg)
    lon = np.arange(70, 140 + deg, deg)
    X, Y = np.meshgrid(lon, lat)
    mask = regionmask.mask_geopandas(countries, lon, lat).to_numpy()
    mask[~np.isnan(mask)] = 1

    # colormap: 将两个colormap进行拼接
    cmap_jet = plt.cm.jet
    cmap_gray = plt.cm.gray
    gray = cmap_gray(np.linspace(0, 1, 9))
    colors = np.vstack((gray[8], cmap_jet(np.linspace(0, 1, 9))))

    # load data
    if season is None:
        file = f'../data/static_result/VIS_Min-{condition[0]}to{condition[1]}_yearly.npy'
    else:
        file = f'../data/static_result/VIS_Min{condition[0]}to{condition[1]}_yearly_{season}.npy'

    factor_dict = {'spring': 92, 'summer': 92, 'autumn': 91, 'winter': 90}
    factor = factor_dict.get(season, 365)
    data = np.load(file, allow_pickle=True)

    # plot map: 画地图
    plt.rcParams['font.family'] = ['Microsoft YaHei']
    extents = [70, 140, 15, 55]
    crs = ccrs.PlateCarree()
    fig = plt.figure(figsize=(12, 6))
    ax = fig.add_subplot(111, projection=crs)
    ax.set_extent(extents, crs)

    prov_reader = shpreader.Reader(r'.\bou2_4p.shp')
    nineline_reader = shpreader.Reader(r'.\九段线.shp')

    ax.add_geometries(prov_reader.geometries(), crs, edgecolor='gray', facecolor='none', lw=1)
    ax.add_geometries(nineline_reader.geometries(), crs, edgecolor='gray', facecolor='none', lw=1)

    ax.set_xticks(np.arange(extents[0], extents[1] + 5, 10))
    ax.set_yticks(np.arange(extents[2], extents[3] + 5, 10))
    plt.tick_params(axis='both', which='major', labelsize=18)

    # 数据分布差异较大，使用分位数对数据进行映射
    data = data * mask
    Low_VIS = data.copy()
    datax = data[~np.isnan(data)]
    level = np.unique((np.percentile(datax, np.array([0, 50, 80] + list(range(90, 99, 2)) + [99.9]))))
    if len(level) >=2:
        x_ = np.linspace(level[-2], level[-1], 4)
    else:
        x_ = np.linspace(0, 1, 4)
    level = np.unique(np.concatenate((level[:-2],x_)))
    xx1 = level[:-1]
    xx2 = level[1:]
    CC1 = 0
    CC2 = len(xx1)

    for jj in range(CC2):
        data[np.where((Low_VIS >= xx1[jj]) & (Low_VIS < xx2[jj]))] = jj + 0.5
	
	# 叠加数据层
    norm = mcolors.Normalize(vmin=0, vmax=len(level) - 1)
    cmaps = LinearSegmentedColormap.from_list('Combined', colors, N=len(level) - 1)
    map = ax.pcolormesh(X, Y, data * mask, cmap=cmaps, vmin=CC1, vmax=CC2)

    title_dict = {'spring': '春季', 'summer': '夏季', 'autumn': '秋季', 'winter': '冬季'}
    title_str = title_dict.get(season, '')
    plt.title('2019—2023年{}平均 {}-{}m能见度频次'.format(title_str, condition[0], condition[1]), fontsize=20, pad=10)
    plt.xlabel('Lon(°E)', fontsize=18)
    plt.ylabel('Lat(°N)', fontsize=18)

    # 调整colorbar与图像等高
    cb_ax = inset_axes(ax, width="3%", height="100%", loc='lower left', bbox_to_anchor=(1.01, 0., 1, 1),bbox_transform=ax.transAxes, borderpad=0)
    cbar = plt.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmaps), cax=cb_ax, ax=map)
    cbar.ax.set_title('频次', size=18, fontname='SimHei', pad=10)
    cbar.ax.tick_params(labelsize=15)
    tick_locs = list(range(len(level)))
    tick_labels = np.unique(np.round(level,1)).astype(str)
    cbar.locator = ticker.FixedLocator(tick_locs)
    cbar.formatter = ticker.FixedFormatter(tick_labels)
    cbar.update_ticks()
	
	# 南海小地图
    ax2 = fig.add_axes([0.74, 0.13, 0.1, 0.25], projection=crs)
    ax2.set_extent([104.5, 124, 0, 26], crs=ccrs.PlateCarree())
    ax2.add_geometries(prov_reader.geometries(), crs, edgecolor='gray', facecolor='none', lw=1)
    ax2.add_geometries(nineline_reader.geometries(), crs, edgecolor='gray', facecolor='none', lw=1)
    ax2.pcolormesh(X, Y, data * mask, cmap=cmaps, vmin=CC1, vmax=CC2)

    prov_reader.close()
    nineline_reader.close()

    save_path = '../map/interp_map'
    if not os.path.exists(save_path):
        os.makedirs(save_path)

    if season is None:
        save_name = f'{condition[0]}-{condition[1]}.png'
    else:
        save_name = f'{condition[0]}-{condition[1]}-{season}.png'
    save_file = os.path.join(save_path, save_name)

    plt.savefig(save_file, bbox_inches='tight', dpi=300)


if __name__ == '__main__':
    conditions = [[0, 50], [50, 100], [100, 200], [200, 500], [500, 1000], [0, 200]]
    for condition in conditions:
        # 年平均频次
        map_frequency(condition)
        # 分季节频次
        for season in ['spring', 'summer', 'autumn', 'winter']:
            map_frequency(condition, season)