MATLAB处理shpfile地理信息数据

数据简介：ESRI Shapefile（shp），或简称shapefile，是美国环境系统研究所公司（ESRI）开发的一种空间数据开放格式。该文件格式已经成为了地理信息软件界的一个开放标准，通常可用于描述几何体对象：点、折线与多边形。例如，Shapefile文件可以存储井、河流、湖泊等空间对象的几何位置。除了几何位置，shp文件也可以存储这些空间对象的属性，例如一条河流的名字，一个城市的温度等等。

编程工具：MATLAB2018b

实现功能：通过shaperead函数读取*.shp，使用geoshow绘制行政边界。

访问数据：

江苏省shp数据

截取徐州shp数据

绘图效果：

徐州市县（区）地理信息图

clc;
clear;
clf;

tic;  

% 使用worldmap([],[])先选定范围，再使用textm()
Lat_min = 33.625;
Lat_max = 35;
Lon_min = 116.25;
Lon_max = 118.75;
worldmap([Lat_min Lat_max],[Lon_min Lon_max]);                             %纬度经度范围显示
setm(gca,'MLineLocation',0.125);                                           %设置经度间隔为0.125°
setm(gca,'PLineLocation',0.125);                                           %设置纬度度间隔为0.125°
setm(gca,'MLabelLocation',0.5)                                             %设置经度标签间隔0.5°
setm(gca,'PLabelLocation',0.5)                                             %设置纬度标签间隔0.5°

% 读取江苏省县（区）边界数据
jiangsu = shaperead('D:\TEST\Matlab程序\绘制地图\江苏省精确到乡镇的shp矢量图\市辖县界.shp','UseGeoCoords',true);   
% 从江苏省县（区）边界数据提取徐州市行政区域边界数据
xuzhou = [jiangsu(9),jiangsu(10),jiangsu(11),jiangsu(16),jiangsu(17),jiangsu(94),jiangsu(1),jiangsu(5)];  

% 画出所有的县（区）
for i = 1:length(xuzhou)
    geoshow(xuzhou(i).Lat,xuzhou(i).Lon);
    hold on;                                                               %保留已经绘制的图像不被覆盖
end
grid on;
% 求解各县（区）中心点，并画在地图上，添加文字标注
center_Lon = zeros(8,1);                                                   %初始化中心点位置矩阵
center_Lat = zeros(8,1);   
name = cell(8,1);                                                          %初始化区（县）名称矩阵

% 执行循环操作
for i = 1:length(xuzhou)
    % 求中心点经度
    a = xuzhou(i).Lon;
    b = a(~isnan(a));                                                      %xuzhou(i).Lon中含有NaN，影响平均值计算，故删除掉NaN
    center_Lon(i) = mean(b);
    % 求中心点纬度
    c = xuzhou(i).Lat;
    d = c(~isnan(c));                                                      %xuzhou(i).Lat中含有NaN，影响平均值计算，故删除掉NaN
    center_Lat(i) = mean(d);
    % 求中心点名称
    name(i) = cellstr(xuzhou(i).NAME); 
    % 绘制中心点
    g(i) = geoshow(center_Lat(i),center_Lon(i),'marker','.','MarkerSize',14,'MarkerEdgeColor','r');  
    % 添加文字标记
    h(i) = textm(center_Lat(i),center_Lon(i)+0.025,name(i),'fontsize',9);  %注意：经纬度顺序要和x,y轴对应，否则文字会到显示范围之外
end

%% 对标记位置进行调整
% 1、修改“丰县”
new_center_Lon(7) = center_Lon(7)-0.1;
new_center_Lat(7) = center_Lat(7);
set(g(7),'visible','off');                                                 %设置原有的中心点不可见
geoshow(new_center_Lat(7),new_center_Lon(7),'marker','.','MarkerSize',14,'MarkerEdgeColor','r');  % 重新绘制适合的中心点
set(h(7),'visible','off');                                                 %设置原有的text不可见
textm(new_center_Lat(7),new_center_Lon(7)+0.025,name(7),'fontsize',9);     %重新写入一遍name
% 2、修改“铜山县”
new_center_Lon(1) = center_Lon(1)-0.1;
new_center_Lat(1) = center_Lat(1)-0.15;
set(g(1),'visible','off');                                                 %设置原有的中心点不可见
geoshow(new_center_Lat(1),new_center_Lon(1),'marker','.','MarkerSize',14,'MarkerEdgeColor','r');  % 重新绘制适合的中心点
set(h(1),'visible','off');                                                 %设置原有的text不可见
textm(new_center_Lat(1),new_center_Lon(1)+0.025,'铜山区','fontsize',9);    %重新写入新的text，将“铜山县”改为“铜山区”

%% 将观测点显示在地图上
[Station_data,Station_text] = xlsread('D:\TEST\Matlab程序\绘制地图\江苏省精确到乡镇的shp矢量图\区域气象自动站信息表.xls');
Station_No_and_Name = Station_text(2:length(Station_text),1);              %站名和站号
% 将站号和站名分离开
Station_No_and_Name_temp = char(Station_No_and_Name);                      %将cell转化为字符串
Station_No = Station_No_and_Name_temp(:,1:5);                              %获取自动站号
Station_Name = Station_No_and_Name_temp(1:length(Station_No_and_Name_temp),7:11);  %获取自动站名

% 将度分秒转化为带小数的度
Station_Lon_char = char(Station_text(2:length(Station_text),8));           %经度
Station_Lat_char = char(Station_text(2:length(Station_text),9));           %纬度
Station_Lon_degree = str2num(Station_Lon_char(:,1:3));                     %获取经度的  度 °     
Station_Lon_min = str2num(Station_Lon_char(:,5:6));                        %获取经度的 分 ′
Station_Lon_sec = str2num(Station_Lon_char(:,8:9));                        %获取经度的 秒 ″
Station_Lat_degree = str2num(Station_Lat_char(:,1:2));                     %获取纬度的  度 °     
Station_Lat_min = str2num(Station_Lat_char(:,4:5));                        %获取纬度的 分 ′
Station_Lat_sec = str2num(Station_Lat_char(:,7:8));                        %获取纬度的 秒 ″
Station_Lon = dms2degrees([Station_Lon_degree Station_Lon_min Station_Lon_sec]); %区域自动气象站坐标的小数形式的经度
Station_Lat = dms2degrees([Station_Lat_degree Station_Lat_min Station_Lat_sec]); %区域自动气象站坐标的小数形式的纬度

% 绘制自动站散点
for i = 1:length(Station_data)
    geoshow(Station_Lat(i),Station_Lon(i),'marker','.','MarkerSize',5,'MarkerEdgeColor','k');
end
% set(Station,'visible','off');                                            %设置区域气象自动站点不可见

toc;