Matlab实现Holland风场

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

---

主要的参数设定：

---

生成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;
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);
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74