爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 4677|回复: 1

[源程序] 雷达山体阻挡分析

[复制链接]

新浪微博达人勋

发表于 2017-4-18 22:08:24 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
各位大神,可以帮忙看下我的程序哪出问题了吗?
我是要做雷达山体阻挡分析的,可每次运行程序都出不来结果,我的循环结构还可以改进吗?


%DEM read
DEM=imread('E:\MATLAB\DEM\Guangdong 90\DEM_Guangdong90m.tif');
lonDEM=109.999583817611++0.000833333332+(0:1:12002)*0.00083333333;
latDEM=25.0004171267299-0.00083333333/2-(0:1:6001)*0.00083333333;

%Radar site informationfor guangzhou
RDST.Name='guangzhou';
RDST.lon=113.3550;
RDST.lat=23.00389;
RDST.El=180.6;
h0=RDST.El/1000;%雷达高程

%雷达山体阻挡分析
Stzd=zeros(9,360,230);
Re=8500.0;
%elv=FxAgl;
elv=[0.5,1.5,2.4,3.4,4.3,6.0, 9.9,14.6,19.5];
Els=elv*pi/180;
r=(1:1:230)-0.5;
az=0:1:359;

for i=1:9
     iEl=Els(i); % Do the first elvation scan analysis

  for j=1:360 %interate azimuth angle
     iAz=az(j);

     for k=1:230;  %interate the radar gate along the radar beam

         kR=r(k);
         h=sqrt(kR*kR+Re*Re+2*kR*Re*sin(iEl))-Re;   
         kH=h+h0;   % Now get the elevation for kth radar gate
         Rs=Re*asin(kR*cos(iEl)/(Re+kH));   %obtain the surface distance from the radar site and the point of kth radar gate

         [iLon,iLat]=fnt_RangeAzToLonLat(RDST.lon,RDST.lat, iAz,Rs);   %obtain the DEM at location (iLon,iLat)
         iDEM=fnt_GetDEM(iLon,iLat);
         %detLat=abs(Lat-ilat);
         %detLon=abs(Lon-ilon);
         %m=find(detLat==min(detLat));
         %n=find(detLon==min(detLon));
         %iDEM=DEM(m,n);

         if kH< iDEM    %block
          Stzd(i,j,k:230)=1;
         else        %  kH>iDEM  not block
          Stzd(i,j,k)=0;
         end
      end
   end
end
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表