Matlab实现Holland风场

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

A-Chin

4027人浏览 · 2022-08-13 08:00:00

A-Chin · 2022-08-13 08:00:00 发布

更新：

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

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

主要的参数设定：

参数	公式
空气密度	$ρ=1.15kg⋅m−3\rho=1.15 kg·m^{-3}$
外围环境气压	$P_n = 1010.0 hPa$
气压差	$ΔP=Pn−Pc\Delta P = P_n - P_c$
科氏力参数	$sin(Lat_c * \pi / 180.0)$
最大风速半径	$Rm=−18.18∗log(ΔP)+112.2R_m = -18.18 * log(\Delta P) + 112.2$
Holland B系数	$B = 1.881 - 0.00557 * R_m - 0.01295 * \|Lat_c\|)$
气压场	$Pr=Pc+ΔP∗exp(−(Rm/r)B)P_r = P_c + \Delta P * exp(-(R_m / r) ^ B)$
梯度风场	$Vg(r)=ΔPBρ(Rmr)Bexp(−(Rmr)B))+(r⋅f2)2+r⋅∣f∣2V_g(r)=\sqrt{\Delta P \frac{B}{\rho}(\frac{R_m}{r})^Bexp(-(\frac{R_m}{r})^B))+(\frac{r·f}{2})^2} + \frac{r·\|f\|}{2}$
近地表10m处的平均风速	$V_{10m}(r) = K_m * V_g(r)$
$K_m$ 取值标准	$Km={0.81,Vg<60.81−2.96∗10−3∗Vg,6≤Vg<19.50.77−4.31∗10−3∗Vg,19.5≤Vg<450.66,Vg≥45K_m=\begin{cases}0.81, V_g < 6 \\0.81 - 2.96 * 10 ^ {-3} * V_g, 6\leq V_g<19.5 \\0.77 - 4.31 * 10 ^ {-3} * V_g, 19.5\leq V_g<45\\0.66, V_g \ge 45\end{cases}$

生成0.25°分辨率的风场：
在这里插入图片描述

附Matlab代码：
其中，basetif.tif是一个0.25°的影像
中心经纬度以及中心气压：
$P_c = 995 hPa$
$Lat_c = 10.125 °$
$Lon_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);

2048 AI社区

有“AI”的1024 = 2048，欢迎大家加入2048 AI社区

更多推荐

MCP服务的两种部署方式

是一个基于 Model Context Protocol (MCP) 的服务器，可以让 AI 助手通过标准化接口操作 Microsoft Word 文档。文档创建与管理内容添加（标题、段落、表格、图片、列表等）文本格式化（字体、颜色、样式）表格格式化（边框、单元格合并、对齐）文档保护与注释提取Word 转 PDFstdio 模式：适合本地应用（Claude Desktop、Lingma 编辑器），