海洋科学方向MATLAB代码快查

内容引自知乎 https://zhuanlan.zhihu.com/p/554465522?utm_id=0

1 数据处理

1.1平均值和标准差
平均值

求平均值一般用mean函数,不推荐用nanmean。对于矩阵A:

miu = mean(A,'all');%求其总的平均值。
miu = mean(A,3);%求其第3维度的平均值。
miu = mean(A,[2,3]);%求其第2、3维度的平均值。
miu = mean(A,'all','omitnan');%求A总的平均值,忽略nan值。
标准差

求标准差一般用std函数。对于矩阵A:

sigma = std(A,0,3);%求A第3维的标准差。0为权重w的默认值。
最大最小值
求最大最小值一般用max和min函数。值得注意的是其函数内要加[]。

M = max(A,[],2,'omitnan');%求A在第2维度的最大值,忽略nan值。
M = min(A,[],[1,3]);%求A在第1、3维度的最小值。
RMSE

误差就是估计值减去实测值(或真值)。误差大就是高估,误差低就是低估,哪个减哪个不能搞混淆。
在这里插入图片描述

e = data-data_insitu;
RMSE = sqrt(sum(e.^2)/length(e));
网格数据插值

interp2 网址:meshgrid 格式的二维网格数据的插值 - MATLAB interp2 - MathWorks 中国

二维插值一般使用interp2函数。

假设需要将4km空间分辨率的数据产品插值至0.25°。
那么产品1为:

Reso1 = 1/24;
lon1 = [0+Reso1/2:Reso1:360-Reso1/2]';
lat1 = flipud([-90+Reso1/2:Reso1:90-Reso1/2]');
[LON1,LAT1] = meshgrid(lon1,lat1);
Data1 = Data1;%4km空间分辨率,8640*4320

产品2为:

Reso2 = 0.25;
lon2 = [0+Reso1/2:Reso1:360-Reso1/2]';
lat2 = flipud([-90+Reso1/2:Reso1:90-Reso1/2]');
[LON2,LAT2] = meshgrid(lon2,lat2);

插值时需要meshgrid以后的网格可能和Data1的矩阵大小不一致,可能需要转置.
插值代码如下:

Data2 = interp2(LON1,LAT1,Data1,LON2,LAT2);

最终得到了网络大小为1440*720的产品2。

但是需要注意的是,低分辨率产品插值到高分辨率,可能存在数据边缘数据缺失的问题,关于这个问题暂时没有比较好的解决方法,或许可行的方法是对低分辨率数据的边缘一个像元先做数据均值填充,再做插值。

1.2 读写文件

一般需要读取的文件格式有txt、csv、nc、hdf、mat。这里读都比较好读,但是nc和hdf的写入比较麻烦,有时间专门详细写一下。

读写csv

读写纯数字csv格式的文件推荐用readmatrix、writematrix、readtable、writetable函数。

这里展示用readmatrix来读csv,因为用readtable读进来是table格式我不是很喜欢用。detectImportOptions用来获取csv文件的信息,用以readmatrix的调用。preview函数可以在命令行打印这个csv文件的前几行内容。我们假设要读的这个csv有10列数据,除了第5列都是数字,那么就可以这么读:将opts的SelectedVariableNames设置为[1:4 6:10],读的时候就会跳过第5列。而读第5列时readmatrix加上输出格式’OutputType’为’string’,就可以读字符格式了。最后如果第5行是’mm/dd/yyyy’格式的时间,就可以用最后两行代码进行转换。

filename='F:\my_csv_file.csv';
%preview data
opts = detectImportOptions(filename);
preview(filename,opts)
%read data
opts.SelectedVariableNames = [1:4 6:10]; 
Data=readmatrix(filename,opts);
opts.SelectedVariableNames = [5]; 
DATE=readmatrix(filename,opts,'OutputType','string');
formatIn='mm/dd/yyyy';
[YEAR MONTH DAY]=datevec(DATE,formatIn);

然后展示writetable写数据,优点是写出来的文件有头部信息。这里columns是表格里参数的名称,建议写英文格式,如果有汉字,需要在writetable函数的’Encoding’项设置为中文excel常用的’GB2312’编码格式(matlab默认输出UTF-8),这样excel打开中文不会乱码。data里的参数需要长度相等。如果其列数和columns里的变量名数量不一致,可能会出现自动编序号的情况,也可能报错,所以尽量保持一致。

%writetable写csv文件
columns = {'var1','var2','var3','var4'};
data = table(var1,var2,var3,var4,'VariableNames', columns);
filename='D:\my_csv_file.csv';
writetable(data,filename,'Encoding','GB2312');%'Encoding'正常情况下不用设置
读写txt
读写nc

ncread很好用

读写hdf

hdf有hdf4和hdf5两种格式,读取方法不一样

读写mat

读写mat很简单了。这些代码基本够用,其中’-v7.3’是保存较大的文件时需要填写的参数,凭感觉加或者不加。

filename = "data.mat";%也可以用绝对路径
load(filename);%读文件中所有参数
load(filename,'var1','var2','var3');%仅读取文件中有限参数

filename = "data.mat";
save(filename,'var1','var2','var3');%保存3个参数
save(filename,'var1','var2','var3','-v7.3');%以-v7.3格式保存mat文件
字符串处理
时间数据处理
曲线拟合

2 画图

2.1 透明度绘图

以用plot绘制Data为例,可用以下代码将红色线的透明度设置为0.1。

a = plot(Data,'color',[1 0 0],'LineWidth',2);
a.Color(4) = 0.1;
2.2 设置多子窗

现在一般用tiledlayout函数。在每个nexttile下面的画图语句将出现在对应的坐标系中。

t = tiledlayout(1,3);%创建1*3的坐标区。
nexttile
%画图1
nexttile
%画图2
nexttile
%画图3
title(t,"Title",'fontsize',16);%调用句柄t以设置总title

当然也可以用subplot函数。但这个函数是R2006a 之前推出,相比R2019b推出的tiledlayout,功能更少。

subplot(2,2,1);%创建2*2的分块中的坐标区,并画第1个图。
%画图1
subplot(2,2,2);%切到第2个图位置。
%画图2
subplot(2,2,[3 4]);%切到第3个图位置,通过设置3和4的位置使图占据2*2的下半空间。
%画图3
2.3 设置ColorBar

ColorBar 网址:显示色阶的颜色栏 - MATLAB colorbar - MathWorks 中国

ColorBar属性 网址:颜色栏的外观和行为 - MATLAB - MathWorks 中国

设置colorbar函数建议带上句柄,方便更进一步单独设置字符内容、字体和字号。代码如下:

h = colorbar('eastoutside');%设置colorbar及其位置,默认为east outside,可以不写。
h.Label.String = "My Colorbar Label";%沿颜色栏添加一个文本标签。
h.Label.FontSize=16;%修改文本标签字体大小。
h.FontSize=16;%修改刻度字体大小。
2.4 设置Colormap

Colormap 网址:查看并设置当前颜色图 - MATLAB colormap - MathWorks 中国

设置colormap实际上很简单。

colormap("jet");%里面填上官方的色标名字就好了
colormap(Colormap);%Colormap可以是和官方色标里的名字,也可以是你自己定义的n×3颜色矩阵,数值在0~1之间

Jet = jet(8);%官方色标都有对应的函数,可以用如jet(8)的形式创建只有8个颜色的色标,画出来的colorbar就是分成8段的
colormap(Jet);

配上各种官方配色的图片:
parula: 在这里插入图片描述
turbo: 在这里插入图片描述
hsv:在这里插入图片描述

hot:

cool:

spring:

summer:

autumn:

winter:

gray:

bone:

copper:

pink:

jet:

lines:

colorcube:

prism:

flag:

white:

当然你也可以用第三方的色标函数包,这里推荐几个:

Custom Colormap 网站:Custom Colormap

PyColormap4Matlab 网站:PyColormap4Matlab

2.5 设置Legend

Legend 网址:在坐标区上添加图例 - MATLAB legend - MathWorks 中国

Legend属性 网址:图例的外观和行为 - MATLAB - MathWorks 中国

2.6 设置字体

当前坐标系字体设置为Times New Roman,字号设置为16。这个函数会把图里所有字符都该了字号和字体,所以用的时候先后顺序要注意。*gca为当前坐标区或图。

而其他涉及到字符打印的代码段,一般也可在括号中进行类似的设置。

set(gca,'Fontname','Times New Roman','Fontsize', 16);

2.7设置图窗大小

一般采用归一化的单位。中括号中内容为[left bottom width height]。即图窗距离左边,底边的距离,以及图窗的宽度和高度。*gcf为当前图窗的句柄。

set(gcf,'unit','normalized','position',[0.03 0.03 0.7 0.8]);

3 m_map相关

3.1 m_map官网:

A Mapping package for Matlab
​www.eoas.ubc.ca/~rich/map.html

3.2 计算地球表面两点距离

这里用到的是m_map中的m_lldist函数,输入两点的经纬度,如上海的121°E,31°N,北京的116°E,40°N,即可计算得到两地的距离dist (km)。

%example
lon = [121,116];
lat = [31,40];
dist = m_lldist(lon,lat);

输出:

dist = 1.0991e+03
3.3 GSHHS高精度海岸线绘图

下载及安装
GSHHS下载链接:Index of /mgg/shorelines/data/gshhs

下载latest的gshhg-bin-2.3.7.zip文件(版本号以最新的为准)。下载后将解压文件中的*.b文件拷贝至m_map文件夹的private子文件夹中。m_map1.4n以后版本须拷贝至data子文件夹中。

3.4 绘图

一般用m_gshhs函数代替m_coast函数即可。

m_gshhs('hc1','patch',[0.7 0.7 0.7],'edgecolor',[0 0 0]);
1:1对比图
Xdata = Xdata;%X参数为列向量形式,数据大小与Y参数一致
Ydata = Ydata;%Y参数为列向量形式,数据大小与X参数一致
scatter(Xdata,Ydata,16,[1 0 0],'filled');%16为散点大小,[1 0 0]是散点颜色,'filled'为填充散点
box on%四边都有边框
axis equal%横纵坐标长度一致,画图画的四四方方的
xlabel("X");%一般为实测值
ylabel("Y")%一般为估计值
xlim([0 100]);
ylim([0 100]);%一般X,Y轴的范围一致
长时间序列图
画地图散点图
画地图伪彩图
%以NASA 4km级产品经纬度为例
Res = 1/24;
lon = [-180+Res/2:Res:180-Res/2]';
lat = fliplr([-90+Res/2:Res:90-Res/2])';%很多产品纬度都是从大到小
%如果经纬度不是矩阵形式,需要用meshgrid函数处理
[LON,LAT] = meshgrid(lon,lat);

Data = nan(8640,4320);
%设置画图的范围
Lonlim = [0 360];
Latlim = [-90 90];
%以太平洋为中心则最好偏移20°,从非洲大陆切开
Bias = 20;
m=max(size(Lon));
n=Bias/360*m;
%按20°对所有经纬度和数据偏移
ind=[n+1:m 1:n];
LON = LON(:,ind);
LON(LON<0) = LON(LON<0)+360;
LAT = LAT(:,ind);
Data = Data(ind,:);
%设置Colormap,伪彩图色值范围,colorbar标签信息,字体大小,标题
Colormap = "jet";
Caxis = [0 1];
LabelString = "LabelString";
Fontsize = 12;
Title = "Title";

%画图,miller指画图的投影方式
m_proj('miller','lon',Lonlim+Bias,'lat',Latlim);
m_pcolor(LON,LAT,Data');%注意Data是否需要转置
m_coast('patch',[0.7 0.7 0.7],'edgecolor','none');%设置陆地颜色及岸线颜色,这里岸线为无色
m_grid('box','fancy','tickdir','in','linewi',2,'fontsize',Fontsize+4);%设置坐标
colormap(Colormap);
h=colorbar('eastoutside');
h.FontSize=Fontsize-2;
caxis(Caxis);
h.Label.String = LabelString;
h.Label.FontSize=Fontsize;
title(Title,'fontsize',Fontsize+4);
画地图散点图带地形
画地图流场矢量图
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值