M_map(五)

一、圆形区域的画图

1、

clear all

LATLIMS=[14 22];
LONLIMS=[108 118];%南海边界范围
m_proj('miller','lon',LONLIMS,'lat',LATLIMS);%投影模式
m_coast('patch',[0.1 0.7 0.4]);%绿色填充
m_grid('box','fancy','tickdir','in');%没有网格，边框相间，%m_line(lon,lat,'linewi',2,'color','r','linestyle',':');控制线条格式，点画线还是直线

lon=112:1:116;
lat=16:1:20;
m_line(lon,lat,'linewi',2,'color','r');%线宽，2；颜色

[X,Y]=m_ll2xy(117,21);
line(X,Y,'marker','.','markersize',24','color','r')%画点

2、距离

clear all
%distance用法
arclen=distance([37,0],[35,0])%返回两点间的相对球心的角度，以度为单位
d=distance([37,0],[35,0],6371)% [纬度，经度] [纬度，经度] [半径]
D=(arclen/180)*pi*6371
%m_map中函数
dist=m_lldist([0 0],[35 37])%[经度 经度] [纬度 纬度]

3、

%画出矩形散点图和圆形图
clear all

LATLIMS=[14 22];
LONLIMS=[108 118];%南海边界范围
m_proj('lambert','lon',LONLIMS,'lat',LATLIMS);%投影模式
m_coast('patch',[0.1 0.7 0.4]);%绿色填充
m_grid('box','fancy','tickdir','in');%没有网格，边框相间，%m_line(lon,lat,'linewi',2,'color','r','linestyle',':');控制线条格式，点画线还是直线

load EDH_south_sea_2008
load coordi_south_sea_2008 

m_range_ring(114.0623,17.9532,[1e2:1e2:3e2],'linewi',2,'color','b');%红色300km范围圆圈

% 矩形点阵
range_lat=4:24;%21N和15N对应的位置下标
range_lon=20:40;%111.5E和116.5E对应的下标
for i=1:length(range_lon)
    for j=1:length(range_lat)
        [X,Y]=m_ll2xy(lon_south_sea(range_lon(i)),lat_south_sea(range_lat(j)));%化为x,y坐标
        line(X,Y,'marker','.','markersize',10,'color','k')%画点
        hold on
    end
end

%离散圆
[X0,Y0]=m_ll2xy(114.0623,17.9532);%圆心化为x,y坐标
line(X0,Y0,'marker','.','markersize',15,'color','r');%画圆心
DIST=m_lldist([114 115],[18 18])%经度加1度，增加的距离
R=300;%300km
[X1,Y1]=m_ll2xy(114.0623-R/DIST,17.9532);%找到300km的一个点

r=sqrt((X0-X1)^2+(Y0-Y1)^2);%地图距离到图上距离转换
r=linspace(0,r,10);%等分成9份
for i=1:10
    theta=0:pi/45:2*pi;%360度，91个点。0到2pi
    x=X0+r(i)*cos(theta); %（X0,Y0）圆心
    y=Y0+r(i)*sin(theta);
    plot(x,y,'.','markersize',6,'color','r')
end

　　

4、pcolor

clear all
n =18;
r = (0:n)'/n;
theta = pi*(-n:n)/n;
X = r*cos(theta);
Y = r*sin(theta);
C = r*cos(2*theta);
pcolor(X,Y,C)
axis equal tight
colorbar

figure
load PCOLOR %南海坐标和波导高度数据
colormap('jet');
shading flat;%平滑方式
gca=pcolor(Plg,Plt,EDH_south_sea)
set(gca, 'LineStyle','none');%去除网格
axis equal tight %按比例展示
colorbar %颜色条

clear all

LATLIMS=[2 22];
LONLIMS=[105 121];%南海边界范围
m_proj('miller','lon',LONLIMS,'lat',LATLIMS);%投影模式


load coordi_south_sea_2008 %载入南海经纬度信息
load EDH_south_sea_2008
for i=1:length(lon_south_sea) %经度
    [X,Y]=m_ll2xy(lon_south_sea(i),12);%化为x,y坐标,12是纬度的平均值
    lon(i)=X;
end
for i=1:length(lat_south_sea) %经度
    [X,Y]=m_ll2xy(113,lat_south_sea(i));%化为x,y坐标，113是经度的平均值
    lat(i)=Y;
end

[Plg,Plt]=meshgrid(lon,lat);%经纬度 形成网格

colormap('jet');
shading flat;%平滑方式
gca=pcolor(Plg,Plt,EDH_south_sea)
set(gca, 'LineStyle','none');%去除网格
axis equal tight %按比例展示
colorbar %颜色条

m_coast('patch',[0.1 0.7 0.4]);%绿色填充
m_grid('box','fancy','tickdir','in');　

左右两幅图对比，左边是将经纬度换为xy值后画的图，更改绘图模式，利用pcolor画的，看出上边边缘有的吻合不是太好；右图是直接利用m_pcolor画的，吻合度更好

相当于载入了一幅地图，先将其经纬度化为正常的坐标，然后再在上面画图。

5、

clear all
n =6;
r = (0:n)'/n;%0到6，半径上均分的数
theta = pi*(-n:n)/n;%将整个圆分成了13分。
X = r*cos(theta);
Y = r*sin(theta);
C = r*cos(2*theta);
pcolor(X,Y,C)
axis equal tight
colorbar

角度逆时针转

半径分成了6份。

 6、

%温度绘图，只读取南海数据绘图
clear all
load EDH_south_sea_2008
load coordi_south_sea_2008 
LATLIMS=[14 22];
LONLIMS=[108 118];%南海边界范围
m_proj('lambert','lon',LONLIMS,'lat',LATLIMS);%投影模式

hold on %一定要有，否则地图会被覆盖

%% 
%将经纬度换成坐标轴
for i=1:length(lon_south_sea) %经度
    [X,Y]=m_ll2xy(lon_south_sea(i),18);%化为x,y坐标,12是纬度的平均值
    lon(i)=X;
end
for i=1:length(lat_south_sea) %经度
    [X,Y]=m_ll2xy(114,lat_south_sea(i));%化为x,y坐标，113是经度的平均值
    lat(i)=Y;
end

%% 
%圆形点图数据
%m_range_ring(114.0623,17.9532,[3e2],'linewi',2,'color','b');%红色300km范围圆圈
%离散圆
[X0,Y0]=m_ll2xy(114.0623,17.9532);%圆心化为x,y坐标
%[X0,Y0]=m_ll2xy(114.0623,5);%圆心化为x,y坐标
%line(X0,Y0,'marker','.','markersize',15,'color','r');%画圆心
DIST=m_lldist([114 115],[18 18]);%经度加1度，增加的距离
R=300;%300km
[X1,Y1]=m_ll2xy(114.0623-R/DIST,17.9532);%找到300km的一个点
%[X1,Y1]=m_ll2xy(114.0623-R/DIST,5);%找到300km的一个点
theta=0:pi/45:2*pi;%360度，91个点。0到2pi
r=sqrt((X0-X1)^2+(Y0-Y1)^2);%地图距离到图上距离转换
r=linspace(0,r,10);%等分成10份
%theta=0:pi/45:2*pi;%360度，91个点。0到2pi

%% 
%定义装位置和数值的空矩阵
circle_x=zeros(length(r),length(theta));
circle_y=zeros(length(r),length(theta));
circle_EDH=zeros(length(r),length(theta));
%% 
%N*2D方式，一共生成了91条线，每条线上有10个数据，
for i=1:length(r)
    x=X0+r(i)*cos(theta); %（X0,Y0）圆心，得到想要获取的EDH的位置（x,y）
    y=Y0+r(i)*sin(theta);
    for j=1:length(theta)
        [min_lon,lon_index]=min(abs(lon-x(j)));%得到最匹配的x轴的值
        [min_lat,lat_index]=min(abs(lat-y(j)));%得到最匹配的y轴的值
        circle_EDH(i,j)=EDH_south_sea(lat_index,lon_index);%存储蒸发波导高度
    end
    circle_x(i,:)=x;
    circle_y(i,:)=y;
    %plot(x,y,'.','markersize',6,'color','r')
end
%% 
%画图
colormap('jet');
shading flat;%平滑方式
gca=pcolor(circle_x,circle_y,circle_EDH)
set(gca, 'LineStyle','none');%去除网格
axis equal tight %按比例展示
% colorbar %颜色条

c=colorbar('eastoutside','fontsize',12)
%caxis([0,14])%颜色条范围限制
title( '2008年南海平均蒸发波导高度','fontsize',15);%[a b c d]确定绘图区域的位置, [a b]为绘图区域左下点的坐标。
m_coast('patch',[0.1 0.7 0.4]);%绿色填充
m_grid('box','fancy','tickdir','in');         %c，d分别为绘图区域的宽和高。

这样可以看到，靠近岸的地方温度低，靠近海的地方温度高。

7、获取路径损失数据图

%获取圆形区域路径损失数据
clear all;
% load testUsrDef.mat
tic
freq = 8000;    % in MHz频率
thetabw = 2;    % in deg 3分贝波束宽度
thetae = 0; % in deg 仰角度数
polrz = 1;  % 1 for horizontal polarization 水平极化, 2 for vertical polarization 极化方式
tx_height = 3;  % tx_antenna height, in m，天线高度
range = 300; % calcu range, in km，水平传播距离
zmax_user = 100;    % maximum height (max desired calculation height), in m，电波限高
rx_height=3;%接收天线高度
delx = 100; % range step, in m，距离步进值迭代设置
delz = 0.1;    % altitude step, in m，高度步进值
%%
%%地形设置
edge_range = [ 20 30 50 ];   % array of edge range, in km，地形限制
edge_height = [ 5 20 7 ];    % array of edge heights at edge ranges, in m，地形高度；？？？无效？
terrain_type = 1;   % tarrain type, 1 for no terrain case, 2 for terrain case 
interp_type = 2;    % edge interp type, 2 for linear, 3 for cubic spline
backward = 1;   % 1 for 1-way SSPE, 2 for 2-way SSPE

%% 海洋表面参数
ground_type = 1;    % ground type, 1 for PEC, 2 for mixed
epsilon = 69.13;    % dielectric constant 介电常数
sigma = 7.146;  % conductivity 电导率
%% 
%预留数据空间
load CIRCLE_EDH  % 10*91 91次循环，91条线，10个蒸发波导高度
circle_path_loss=zeros(10*range+1,length(circle_EDH));
circle_pl_x=zeros(10*range+1,length(circle_EDH));
circle_pl_y=zeros(10*range+1,length(circle_EDH));

%% 波导条件设置
% duct_type = 5;  % duct type, 1 for standard atmosphere, 3 for surface-based duct, 5 for evaporation duct
% duct_M =  [360 0];%必须带0；
% duct_height = [15 0]; %必须带0
% duct_range = 1; %？？？水平不均匀性的设置

% duct_type = 5;  % duct type, 1 for standard atmosphere, 3 for surface-based duct, 5 for evaporation duct
% duct_M =  [360 0 ;300 0;344 0];
% duct_height = [15 0;18 0; 20 0]; 
% duct_range =[0 100 200]; %水平不均匀性的设置,0km处蒸发波导高度15m,100km处蒸发波导高度5m,200km处蒸发波导高度20m.
for i=1:length(circle_EDH)
    duct_type = 5;  % duct type, 1 for standard atmosphere, 3 for surface-based duct, 5 for evaporation duct
    
    duct_M=zeros(10,2);
    
    duct_height = zeros(10,2);
    duct_height(:,1)=circle_EDH(:,i);%获取第i列（线）波导高度
    duct_range =linspace(0,300,10); %水平不均匀性的设置,0km处蒸发波导高度15m,100km处蒸发波导高度5m,200km处蒸发波导高度20m.
    %计算路径损失
    [umat, path_loss, prop_fact, free_space_loss, range_vec, z_user, z, stopflag] =  SSPE_function(freq, ...
        thetabw, thetae, polrz, tx_height, range, zmax_user, edge_range, edge_height, ...
        duct_type, duct_M, duct_height, duct_range, terrain_type, interp_type, backward, ...
        ground_type, epsilon, sigma, delx, delz);
    circle_path_loss(:,i)=path_loss(rx_height*10,:);
    
end
%计算坐标
for i=1:length(circle_EDH)
   circle_pl_x(:,i)=linspace(circle_x(1,i),circle_x(10,i),10*range+1);
   circle_pl_y(:,i)=linspace(circle_y(1,i),circle_y(10,i),10*range+1);
end
toc
save CIRCLE_PATH_LOSS circle_path_loss circle_pl_x circle_pl_y

8、圆形路径损失图

%圆形路径损失绘制
clear all
LATLIMS=[14 22];
LONLIMS=[108 118];%南海边界范围
m_proj('lambert','lon',LONLIMS,'lat',LATLIMS);%投影模式
hold on %一定要有，否则地图会被覆盖
%% 
load CIRCLE_PATH_LOSS
%% 
%画图
colormap('jet');
shading flat;%平滑方式
gca=pcolor(circle_pl_x,circle_pl_y,circle_path_loss);
set(gca, 'LineStyle','none');%去除网格
axis equal tight ;%按比例展示
%colorbar %颜色条

c=colorbar('eastoutside','fontsize',12);
caxis([100,220]);%颜色条范围限制
m_coast('patch',[0.1 0.7 0.4]);%绿色填充
m_grid('box','fancy','tickdir','in');         %c，d分别为绘图区域的宽和高。