Python结合Matlab气候数据工具箱CDT

#!/usr/bin/env python
from ecmwfapi import ECMWFDataServer
server = ECMWFDataServer()
server.retrieve({
    "class": "ei",
    "dataset": "interim",
    "date": "20170101/20170201/20170301/20170401/20170501/20170601/20170701/20170801/20170901/20171001/20171101/20171201",
    "expver": "1",
    "grid": "0.75/0.75",
    "levtype": "sfc",
    "param": "134.128/165.128/166.128/167.128/228.128",
    "step": "12",
    "stream": "mnth",
    "time": "00:00:00",
    "type": "fc",
    "area": "90/-180/-90/180",
    "format": "netcdf",
    "target": "F:/Rpython/lp36/ERA_Interim_2017.nc",
})

clc;clear;close all
% 在命令行窗口中以文本形式显示 NetCDF 数据源 source 中的所有组、维度、变量定义，以及所有属性。
ncdisp('F:/Rpython/lp36/ERA_Interim_2017.nc');
data='F:/Rpython/lp36/ERA_Interim_2017.nc';
lat = ncread(data,'latitude');
lon = ncread(data,'longitude');
T = ncread(data,'t2m');
t = ncdateread(data,'time');
% Calculate mean temperature of 2017:
Tm = mean(T,3);
figure(1);
% Plot mean temperature:
imagesc(lon,lat,Tm)
% Rotate and flip T:
T = flipud(rot90(T));
% Recalculate Tm using the rotated, flipped data matrix:
Tm = mean(T,3);
figure(2);
% Re-plot the mean temperature:
imagesc(lon,lat,Tm);
figure(3);
imagesc(lon,lat,Tm);
axis xy
figure(4);
[Lon,Lat] = meshgrid(lon,lat);
pcolor(double(Lon),double(Lat),Tm)
shading interp % eliminates black lines between grid cells
% cmocean thermal % cold-to-hot colormap
xlabel longitude
ylabel latitude
figure(5);
[Lat,Lon,T,Tm] = recenter(Lat,Lon,T,Tm);
pcolor(double(Lon),double(Lat),Tm)
shading interp % eliminates black lines between grid cells
hold on
borders('countries','color',rgb('dark gray'))
cmocean thermal % cold-to-hot colormap
xlabel longitude
ylabel latitude

% Load longitude array:
lon = double(ncread(data,'longitude'));
% determine which indices correspond to dimension 1:
ind1 = find(lon>=-180 & lon<=180);
% Do the same for lat:
lat = double(ncread(data,'latitude'));
ind2 = find(lat>=-90 & lat<=90);
% Clip lat and lon to their specified range:
lat = lat(ind2);
lon = lon(ind1);
% Make a grid:
[Lat,Lon] = meshgrid(lat,lon);
start = [ind1(1) ind2(1) 6];
stride = [length(ind1) length(ind2) 1];
% Load June surface pressure for Mediterranean:
sp = ncread(data,'sp',start,stride);
% Also load temperature and wind:
T = ncread(data,'t2m',start,stride);
u10 = ncread(data,'u10',start,stride);
v10 = ncread(data,'v10',start,stride);
figure(6)
pcolor(Lon,Lat,T-273.15) % temperature in deg C
shading interp
cmocean thermal
caxis([17 44]) % color axis limits
cb = colorbar;
ylabel(cb,'temperature \circC')
hold on
borders('countries','color',rgb('dark gray'))
contour(Lon,Lat,sp,'k') % pressure contours
quiversc(Lon,Lat,u10,v10) % wind vectors
% quiver(Lon,Lat,u10,v10) % wind vectors
xlabel 'longitude'
ylabel 'latitude'

figure(7);
% 作图plot
m_proj('Equidistant Cylindrical','lon',[-180 180],'lat',[-90 90]);
m_coast('patch',[1 .85 .7],'edgecolor',[0 0 0],'linewidth',1);
hold on;
m_pcolor(Lon,Lat,T-273.15) % temperature in deg C
hold on;
m_contour(Lon,Lat,sp,'r') % pressure contours
hold on;
[Lon7,Lat7] = meshgrid(lon,lat);
Lon7=Lon7';
Lat7=Lat7';
u107 = ncread(data,'u10');
v107 = ncread(data,'v10');
T7 = ncread(data,'t2m');
t7 = ncdateread(data,'time');
% Calculate mean temperature of 2017:
Tm7 = mean(T7,3);
U107 = mean(u107,3);
V107 = mean(v107,3);
m_quiver(Lon7(1:10:end,1:10:end),Lat7(1:10:end,1:10:end),U107(1:10:end,1:10:end),V107(1:10:end,1:10:end),0,'k','linewidth',0.5)
m_grid('ytick',[-90:30:90],'xtick',[-180:60:180],'tickdir','out','linest','none','fontname','Times','fontsize',12,'linewidth',1.5);
% 调用ncl色带
load('colorbar-mat/rainbow_r.mat');
colormap(rainbow_r);
hold on
c=colorbar('eastoutside','ticklength',0);
caxis([0,50])
% 调整colorbar
ax = gca;
axpos = ax.Position;
c.Position(3) = 0.5*c.Position(3);
ax.Position = axpos;
cbarrow;