使用Matlab从.shp文件中提取浙江省内的气象站点

1.要求

以浙江省为例，从全国气象站中获取浙江省内的气象站。

2.数据格式

2.1 stations.txt

2014年全国45330个气象站数据

第一列是站点编号，第二列是站点经度，第三列是站点纬度，第四列的站点海拔。

data = importdata('stations.txt');

在工作区可以发现data是一个45330*4的矩阵，说明全国有45330个站点（本文使用的气象站点数据是2014年的）。

2.2 浙江省.shp

浙江省shapefile文件包括浙江省shp、shx、dbf文件

map = shaperead('浙江省.shp');

双击工作区中的map变量，可以查看其内部信息如下：

map是一个结构体数组，其中每个数组元素map(i)均包含4个成员变量Geometry、BoundingBox、X、Y。例如，可以使用map(2).BoundingBox获取得到

这里简单介绍一下，map结构体数组将浙江省划分成了90个多边形（Polygon），BoundingBox表示每个多边形的经纬度上下限，即能完全包裹该多边形的最小矩形。

通过搜索90个BoudingBox，可以初步确定能框选浙江省的最小矩形[最大经度 最小经度 最大纬度 最小纬度]。

箭头标注的点也属于浙江省范围（不排除上面有气象站）

3.提取站点

3.1 初步筛选

因此，我们可以由这个矩形进行初步过滤，得到在矩形框内的所有站点，之后再对每个站点进行判断，判断其是否在浙江省内。

station = importdata('stations.txt');
map = shaperead('浙江省.shp');
bound = [map.BoundingBox];% 将90个矩阵合并为一个矩阵
lon = bound(:,1:2:end);   % 取出矩阵的奇数列，即包含经度的列
lat = bound(:,2:2:end);   % 取出矩阵的偶数列，即包含纬度的列
% [min, max] = bounds(a) 返回一维向量a中的最小值和最大值
% b(:)将一个n行*m列的矩阵b转换为n*m行的一维列向量
[lon_min, lon_max] = bounds(lon(:));
[lat_min, lat_max] = bounds(lat(:));
% index1 经度不在矩形范围内的站点索引
index1 = station(:,2) > lon_max | station(:,2) < lon_min;
% index2 纬度不在矩形范围内的站点索引
index2 = station(:,3) > lat_max | station(:,3) < lat_min;
% 将这些站点从矩阵中去除
station(index1|index2,:)=[];

3.2 精确查找

这里介绍一下inpolygon这个函数：

in = inpolygon(xq,yq,xv,yv) 返回 in，以指明 xq 和 yq 所指定的查询点是否在 xv 和 yv 定义的多边形区域的边缘内部或边缘上

在边缘内部或边缘上则 in = 1, 在多边形外部则 in = 0

map结构体中的X和Y用来描述一个多边形区域，相当于上述的xv和yv。因此，我们只需要把站点的经纬度（xq,yq）带入到函数中，逐个判断其是否在多边形内即可。

index = [];                % 记录不在浙江省内的站点索引
for i = 1:length(station)  % 遍历所有站点
    flag = 0;              % 初始化标记为0，表示当前站点不在浙江省内
    for j = 1:length(map)  % 遍历所有多边形
        % 判断该站点是否在该多边形内
        if inpolygon(station(i,2),station(i,3),map(j).X,map(j).Y)
            flag = 1;      % 在多边形内，更新flag=1，表示当前站点在浙江省内
            break;         % 结束循环
        end
    end
    if ~flag               % 站点不在浙江省内
        index = [index; i];% 记录下站点的索引
    end
end
station(index,:)=[];       % 通过站点索引，删除不在浙江省内的站点

最后得到2014年的2205个浙江省内的站点信息数据。

4. 推广

function station = getStationsOfArea(fileMap, fileStation)
%% 作者微信：qczsbwjzjn
%% 时间：2021.04.22
map = shaperead(fileMap);          % 读取shp文件
station = importdata(fileStation); % 读取站点文件
%% 初步缩小范围
bound = [map.BoundingBox];% 将所有BoundingBox矩阵合并为一个矩阵
lon = bound(:,1:2:end);   % 取出矩阵的奇数列，即包含经度的列
lat = bound(:,2:2:end);   % 取出矩阵的偶数列，即包含纬度的列
% [min, max] = bounds(a) 返回一维向量a中的最小值和最大值
% a(:)将一个n行*m列的矩阵a转换为n*m行的一维列向量
[lon_min, lon_max] = bounds(lon(:));
[lat_min, lat_max] = bounds(lat(:));
% index1 经度不在矩形范围内的站点索引
index1 = station(:,2) > lon_max | station(:,2) < lon_min;
% index2 纬度不在矩形范围内的站点索引
index2 = station(:,3) > lat_max | station(:,3) < lat_min;
% 将这些站点从矩阵中去除
station(index1|index2,:) = [];
%% 精确查找
index = [];                % 记录不在map内的站点索引
for i = 1:length(station)  % 遍历所有站点
    flag = 0;              % 初始化标记为0，表示当前站点不在map内
    for j = 1:length(map)  % 遍历所有多边形
        % 判断该站点是否在该多边形内
        if inpolygon(station(i,2),station(i,3),map(j).X,map(j).Y)
            flag = 1;      % 在多边形内，更新flag=1，表示当前站点在map内
            break;         % 结束循环
        end
    end
    if ~flag               % 站点不在map内
        index = [index; i];% 记录下站点的索引
    end
end
station(index,:) = [];     % 通过站点索引，删除不在map内的站点
end

通过调用函数station = getStationsOfArea(fileMap, fileStation)获取任意.shp区域内的站点

% 作者微信：qczsbwjzjn
% 调用样例
station = getStationsOfArea('浙江省.shp', 'stations.txt');

提供文件demo（含浙江省.shp文件，stations.txt文件，getStationsOfArea.m文件， demo.m文件）下载链接