MATLAB 提取shapefile范围内的雨量站点

问题：有大范围地区的雨量站数据，如何批量提取某个指定地区的雨量站点数据？

这是收集到的逐小时雨量数据格式，分别包含了站码、经纬度、高度、降雨量。

现在需要提取出shapefile范围内的站点。

1、首先，使用shaperead函数读取shapefile文件，识别出shp文件的经纬度信息。

%读取shape经纬度，需要注意的是shp文件坐标为WGS1984，这是用的是广州的shp文件
GZ=shaperead('GZ1984.shp');

下图的X、Y代表了shp的经度、纬度

2、然后提取逐小时数据站点的经纬度。

%读取站点经纬度，经纬度分别存储在第4行以下，第2、3列范围内。
Sx=shuju(4:length(shuju),2);     %经度
Sy=shuju(4:length(shuju),3);     %纬度

3、使用inpolygon函数判断站点经纬度是否在shp范围内。

%判断是否在shape范围内，其中GZ.X和GZ.Y为shp经度、纬度
GZs=inpolygon(Sx,Sy,GZ.X,GZ.Y);

此时会返回一个逻辑变量，1代表站点在shp内，0则代表在范围外，如下图

4、使用find函数确定shp范围内站点的位置，并将位置标记为P

%取出广州站点数据位置
P=find(GZs(:)==1);

%读取shuju中的雨量，其中P为shp范围内的雨量站所在位置

Data=shuju(P,1:5);

%Data就是所需的站点信息，包括站码、经纬度、高程、雨量值

5、可视化。使用plot函数画出提取得到站点的分布情况

%其中 ‘~’为取非，  Sx(GZs),Sy(GZs)为shp范围内的站点，标记为,'r+'，Sx(~GZs),Sy(~GZs)为范围外的站点分布，标记为蓝色的o，'bo'
plot(GZ.X,GZ.Y,Sx(GZs),Sy(GZs),'r+',Sx(~GZs),Sy(~GZs),'bo');