利用matlab筛选登陆我国的台风

 程序所用的台风数据来源于IBTrACS

省界国界资料来源于

matlab绘制世界地图（含国界）、中国地图（含省界），可下载m_map和shp文件_国界省界shp文件_ClimLady的博客-CSDN博客

 程序如下

clear;
clc;clf;
cd('D:\毕业设计\毕设数值模拟');
ChinaL=shaperead('D:\matlab扩充包\全国省级、地市级、县市级行政区划shp\全国省级、地市级、县市级行政区划shp\bou2_4l.shp');
ChinaP=shaperead('D:\matlab扩充包\全国省级、地市级、县市级行政区划shp\全国省级、地市级、县市级行政区划shp\bou1_4p.shp');
bou2_4lx=[ChinaL(:).X];
bou2_4ly=[ChinaL(:).Y];
bou2_4px=[ChinaP(:).X];
bou2_4py=[ChinaP(:).Y];
World = shaperead('D:\matlab扩充包\world map china line\world map china line.shp');%国界数据
wx = [World(:).X];wy = [World(:).Y];
cd('D:\毕业设计\毕设数值模拟');
f1='IBTrACS.ALL.v04r00.nc';
time=ncread(f1,'time');
landfall=ncread(f1,'landfall');%距离陆地的距离
isotime=ncread(f1,'iso_time');%字符式时间
isotimea=isotime(:,:,13539);
basin=ncread(f1,'basin');%读取台风所在海盆
numobs=ncread(f1,'numobs');
lon=ncread(f1,'lon');
lat=ncread(f1,'lat');
number=ncread(f1,'number');
name=ncread(f1,'name');%台风名称
%%
%筛选2000年以后西北太平洋的低压热带气旋
y1=find(isotime(1,1,:)=='2');y1=y1(1:2386);%筛选2000年以后的台风
basin1=basin(:,:,y1);land1=landfall(:,y1);
lon1=lon(:,y1);
lat1=lat(:,y1);
x=find(basin1(1,1,:)=='W'& basin1(2,1,:)=='P');%筛选台风位置（西太平洋）
land2=land1(:,x);lon2=lon1(:,x);lat2=lat1(:,x);
[x2,y2,v]=find(land2==0);[y3,ia,ic]=unique(y2);%筛选登陆的台风
for i=1:length(y2)
    lon3(i)=lon2(x2(i),y2(i));lat3(i)=lat2(x2(i),y2(i));
    [in(i) on(i)]=inpolygon(lon3(i),lat3(i),bou2_4lx,bou2_4ly);1
end%筛选登陆中国的台风
y4=unique(y2((find(in==1))));%登陆中国的台风序号
lat4=lat2(:,y4);lon4=lon2(:,y4);
%%
figure(1);
m_proj('mercator','lon',[100 150],'lat',[0 45]);
m_coast('patch',[0.92 .92 .92],'edgecolor','none');
m_grid('linestyle','-','gridcolor','w','backcolor',[.5294 .80784 0.92156],'ytick',[0 10 20 30 40 45],'xtick',[100 110 120 130 140 150]);
hold on
m_plot(wx,wy,'k','linewidth',0.4,'color','k');%绘制国界
hold on
m_plot(bou2_4lx,bou2_4ly,'k','linewidth',0.1);%绘制中国省界
hold on;
for i=1:length(y4)
    m_plot(lon4(:,i),lat4(:,i));
    hold on
end   
hold off
xlabel('经度','fontsize',12);
ylabel('纬度','fontsize',12);
title('2000年至2022年西北太平洋登录台风轨迹图','fontsize',12);
saveas(gca,'2000年至2022年西北太平洋登录的台风轨迹图.jpg');

所作图 

 其实关键的就那两部分

[x2,y2,v]=find(land2==0);[y3,ia,ic]=unique(y2);%筛选登陆的台风（landfall为0表示已经登陆）
for i=1:length(y2)
    lon3(i)=lon2(x2(i),y2(i));lat3(i)=lat2(x2(i),y2(i));
    [in(i) on(i)]=inpolygon(lon3(i),lat3(i),bou2_4lx,bou2_4ly);1
end%利用省界将所有登陆的台风轨迹点进行筛选，找出登陆中国的台风
y4=unique(y2((find(in==1))));%登陆中国的台风序号
lat4=lat2(:,y4);lon4=lon2(:,y4);

一、利用IBTrCS资料中的landfall向量（即台风眼距离大陆的距离）筛选出所有台风所有登陆的时间点；（所用函数find（二维查找台风登陆时间点）、unique（除去重复的台风编号））

二、判断所有登陆台风的经纬度是否在国界之内，这里主要用到inpolygon函数（判断点是否在不规则闭合区域内）用法可以参照其他博主分享，链接：MATLAB 判断点（坐标）是否在.shp边界（不规则多边形）范围内_matlab怎么判断点在不在封闭图形中_天外来客Jonty的博客-CSDN博客

对于第二部分，我用的时中国省界进行筛选的，因为没找到单独的中国国界，如果用世界国界资料，会将所有登陆台风都筛选进去，

还有一个容易犯的错就是只筛选每个编号台风第一个登陆时间点的经纬度，这也不对，因为有的台风会先经过其他国家再来到中国，所以要筛选所有landfall为零的时间点，我就犯了这个错，结果筛选出来，22年就25个，巨离谱。