Matlab实现Holland风场

更新：

日期	内容
2022.09.13	代码中 Vg_Holland2 = r * abs(f) / 2; 中 r 的单位是 km，需要乘以 1000 转化为 m。

---

仅使用Matlab实现Holland圆对称风场，不考虑热带气旋的移动速度和前进方位角。

---

主要的参数设定：

参数	公式
空气密度	$\rho =1.15kg\cdot m^{-3}$
外围环境气压	$P_n=1010.0hPa$
气压差	$\Delta P=P_n-P_c$
科氏力参数	$f=2\ast 7.292\ast 0.00001\ast sin(La{t}_c\ast \pi /180.0)$
最大风速半径	$R_m=-18.18\ast log(\Delta P)+112.2$
Holland B系数	$B=1.881-0.00557\ast R_m-0.01295\ast |La{t}_c|)$
气压场	$P_r=P_c+\Delta P\ast exp(-(R_m/r{)}^B)$
梯度风场	$V_g(r)=\sqrt{\Delta P\frac{B}{\rho }(\frac{R_m}{r}{)}^{B}exp(-(\frac{R_m}{r}{)}^B))+(\frac{r\cdot f}{2}{)}^2}+\frac{r\cdot |f|}{2}$
近地表10m处的平均风速	$V_{10m}(r)=K_m\ast V_g(r)$
$K_m$取值标准	$$K_m=\{\begin{array}{l}0.81,V_g<6\\ 0.81-2.96\ast 10^{-3}\ast V_g,6\le V_g<19.5\\ 0.77-4.31\ast 10^{-3}\ast V_g,19.5\le V_g<45\\ 0.66,V_g\ge 45\end{array}$$

---

生成0.25°分辨率的风场：

---

附Matlab代码：
其中，basetif.tif是一个0.25°的影像
中心经纬度以及中心气压：
$P_c=995hPa$
$La{t}_c=10.125°$
$Lo{n}_c=120.125°$

%%
clc;
clear;

%%
fullpath = mfilename('fullpath');
[path, name] = fileparts(fullpath);

%%
%
basetif = strcat(path, '\basetif.tif');
[A_M, R_M] =  geotiffread(basetif);

info = geotiffinfo(basetif);
A_I = zeros(size(A_M));
A_J = zeros(size(A_M));

for i = 1 : size(A_I, 1)
    A_I(i, :) = i;
end

for j = 1 : size(A_J, 2)
    A_J(:, j) = j;
end

[A_Lat, A_Lon] = pix2latlon(info.RefMatrix, A_I, A_J);

%%
P_c = 995; % hPa
Lat_c = 10.125; % °
Lon_c = 120.125; % °

%%
%
rou_air = 1.15; % kg m-3

%
P_n = 1010.0;
dertP = P_n - P_c;
Rm = -18.18 * log(dertP) + 112.2; % R2 = 0.84
B = 1.881 - 0.00557 * Rm - 0.01295 * abs(Lat_c);
f = 2 * 7.292 * 0.00001 * sin(Lat_c * pi / 180.0);

% great cricle arc
R_earth = 6371;
arcLAT1 = A_Lat * pi / 180.0;
arcLAT2 = Lat_c * pi / 180.0;
arcLON1 = A_Lon * pi / 180.0;
arcLON2 = Lon_c * pi / 180.0;

cos_sita = sin(arcLAT1) .* sin(arcLAT2) + ...
    cos(arcLAT1) .* cos(arcLAT2) .* cos(arcLON1 - arcLON2);
sita = acos(cos_sita);
r = R_earth * sita + 0.001;

P_r = P_c + dertP * exp(-(Rm ./ r) .^ B);

Vg_Holland1 = 100 * dertP * B / rou_air * (Rm ./ r) .^ B .* exp(-(Rm ./ r) .^ B);
Vg_Holland2 = r * abs(f) / 2 * 1000;
Vg_Holland = (Vg_Holland1 + Vg_Holland2 .^ 2) .^ 0.5 - Vg_Holland2;

%%
Vg = Vg_Holland;

Km = zeros(size(Vg_Holland));
Km(Vg < 6) = 0.81;
Km(Vg >= 6 & Vg < 19.5) = 0.81 - 2.96 * 10 ^ (-3) * (Vg(Vg >= 6 & Vg < 19.5) - 6);
Km(Vg >= 19.5 & Vg < 45) = 0.77 - 4.31 * 10 ^ (-3) * (Vg(Vg >= 19.5 & Vg < 45) - 19.5);
Km(Vg >= 45) = 0.66;

V_10m = Km .* Vg;

outfile = strcat(path, '\V_10m.tif');
geotiffwrite(outfile, single(V_10m), R_M);