利用MATLAB制作DEM山体阴影

这篇具有很好参考价值的文章主要介绍了利用MATLAB制作DEM山体阴影。希望对大家有所帮助。如果存在错误或未考虑完全的地方,请大家不吝赐教,您也可以点击"举报违法"按钮提交疑问。

在地理绘图中,我们使用的DEM数据添加山体阴影使得绘制的图件显得更加的美观。

GIS中使用ArcGIS软件就可以达到这一目的,或者使用GMT,同样可以得到山体阴影的效果。

本文提供了一个MATLAB的函数,可以得到山体阴影。

clear all;clc;close all
load demo.mat
%% draw hillshade
x=x(:,1);
y=y(1,:);
hs=hillshade_esri(z,x,y);
subplot(1,2,1)
imagesc(x,y,z')
axis image
set(gca,'ydir','normal')
title('DEM')
colormap(gray)
caxis([min(z(:)) max(z(:))])
subplot(1,2,2)
imagesc(x,y,hs')
axis image
set(gca,'ydir','normal')
title('Hillshade')
colormap(gray)
hold on
caxis([min(hs(:)) max(hs(:))])
drawnow

利用MATLAB制作DEM山体阴影,matlab,javascript,开发语言

其中调用的函数 hillshade_esri.m如下:

function h = hillshade_esri(dem,X,Y,varargin)
% PUPROSE: Calculate hillshade for a digital elevation model (DEM) based on
%          the algorithm posted on http://edndoc.esri.com/arcobjects/9.2/net/shared/geoprocessing/spatial_analyst_tools/how_hillshade_works.htm
% -------------------------------------------------------------------
% USAGE: h = hillshade_esri(dem,X,Y,varagin)
% where: dem is the DEM to calculate hillshade for
%        X and Y are the DEM coordinate vectors
%        varargin are parameters options
%
% OPTIONS: 
%        'azimuth'  is the direction of lighting in deg (default 315)
%        'altitude' is the altitude of the lighting source in
%                   in degrees above horizontal (default 45)
%        'zfactor'  is the DEM altitude scaling z-factor (default 1)
%        'plotit'   creates a simple plot of the hillshade
%
% EXAMPLE:
%       h=hillshade_esri(peaks(50),1:50,1:50,'azimuth',45,'altitude',100,'plotit')
%       - calculates the hillshade for an example 50x50 peak surface.
%       - changes the default settings for azimuth and altitude.
%       - creates an output hillshade plot

% See also: GRADIENT, CART2POL
%
% Note: Uses ESRIs hillshade algorithm, the output will be the same as the
% output with ArcGIS Hillshade Function.
%
% Felix Hebeler, Dept. of Geography, University Zurich, February 2007.
% modified by Andrew Stevens (astevens@usgs.gov), 5/04/2007
% modified by Wenbin Jiang (jwbalbert@gmail.com), 7/06/2011

%% configure inputs
%default parameters
azimuth=315;
altitude=45;
zf=1;
plotit=0;

%parse inputs
if isempty(varargin)~=1     % check if any arguments are given
    [m1,n1]=size(varargin);
    opts={'azimuth';'altitude';'zfactor';'plotit'};
    for i=1:n1;             % check which parameters are given
        indi=strcmpi(varargin{i},opts);
        ind=find(indi==1);
        if isempty(ind)~=1
            switch ind
                case 1
                    azimuth=varargin{i+1};
                case 2
                    altitude=varargin{i+1};
                case 3
                    zf=varargin{i+1};
                case 4
                    plotit=1;
            end
        end
    end
end

%% Initialize paramaters
dx=abs(X(2)-X(1));  % get cell spacing in x and y direction
dy=abs(Y(2)-Y(1));  % from coordinate vectors

% lighting azimuth
azimuth = 360.0-azimuth+90; %convert to mathematic unit 
azimuth(azimuth>=360)=azimuth-360;
azimuth = azimuth * (pi/180); %  convert to radians

%lighting altitude
altitude = (90-altitude) * (pi/180); % convert to zenith angle in radians

%% calc slope and aspect (radians)
im=length(X);
jm=length(Y);
fx=zeros(im,jm);
fy=zeros(im,jm);
asp=zeros(im,jm);
for i=2:im-1
    for j=2:jm-1
        fx(i,j)=((dem(i+1,j+1)+2*dem(i+1,j)+dem(i+1,j-1))-(dem(i-1,j+1)+2*dem(i-1,j)+dem(i-1,j-1)))/8/dx;
        fy(i,j)=((dem(i-1,j-1)+2*dem(i,j-1)+dem(i+1,j-1))-(dem(i-1,j+1)+2*dem(i,j+1)+dem(i+1,j+1)))/8/dy;
        if fx(i,j)~=0
            asp(i,j) = atan2(fy(i,j),-fx(i,j));
            if asp(i,j)<0
                asp(i,j)=asp(i,j)+2*pi;
            end
        else
            if fy(i,j)>0
                asp(i,j)=pi/2;
            elseif fy(i,j)<0
                asp(i,j)=3*pi/2;
            end
        end     
    end
end
grad = hypot(fx,fy);
grad=atan(zf*grad); %steepest slope

%% hillshade calculation
h = 255.0*( (cos(altitude).*cos(grad) ) + ( sin(altitude).*sin(grad).*cos(azimuth-asp)) ); % ESRIs algorithm
h(h<0)=0; % set hillshade values to min of 0.
h=setborder(h,1,NaN); % set border cells to NaN

%% plot results if requested
if plotit==1
    figure
    imagesc(X,Y,h')
    axis image
    set(gca,'ydir','normal')
    colormap(gray)
end

%% -- Subfunction--------------------------------------------------------------------------
function grid = setborder(grid,bs,bvalue)
grid(1:bs,:)=bvalue; %toprows
grid(size(grid,1)-bs+1:size(grid,1),:)=bvalue; %bottom rows
grid(:,1:bs)=bvalue; %left cols
grid(:,size(grid,2)-bs+1:size(grid,2))=bvalue;

其中有三个参数可以修改:azimuth=315;altitude=45;zf=1;

1.修改 azimuth,the direction of lighting in deg,下图的变化范围为0:360:

利用MATLAB制作DEM山体阴影,matlab,javascript,开发语言

2.修改 altitude,the altitude of the lighting source in degrees above horizontal,下图变化范围为0:180:

利用MATLAB制作DEM山体阴影,matlab,javascript,开发语言

 

3.修改 zf,the DEM altitude scaling z-factor ,下图变化范围为1:50:

利用MATLAB制作DEM山体阴影,matlab,javascript,开发语言文章来源地址https://www.toymoban.com/news/detail-612628.html

到了这里,关于利用MATLAB制作DEM山体阴影的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处: 如若内容造成侵权/违法违规/事实不符,请点击违法举报进行投诉反馈,一经查实,立即删除!

领支付宝红包 赞助服务器费用

相关文章

  • matlab去除图片中的阴影(用高斯模糊估计阴影模式)

    1、思路 2、阴影模式的估计方法 3、代码实现 3.1、通过现成的函数来构造高斯滤波器 3.2、通过原理构建高斯低通滤波器 由于光照等原因造成的图片存在阴影且影响视觉效果的情况,带阴影的图像可以表示为: 其中,是最终的图像,和分别是没有受到影响的原图像和阴影图像

    2024年02月08日
    浏览(39)
  • DEM地形图的制作

    1.数据下载 下载数据主要基于“地理空间数据云”,目前DEM的空间分辨率有两种:30M、90M空间分辨率。本次制作选用 GDEMV2 30M 数据集,空间范围根据行政区选择”云南省“,也可以通过经纬度等其他方式进行选择。然后,对其数据检索并下载。 2.将下载的压缩包文件进行解压

    2024年02月07日
    浏览(49)
  • 39.利用matlab寻找素数(matlab程序)

    1. 简述        MATLAB嵌套循环允许使用一个循环在另一循环内,下面用一个嵌套循环来把所有从1到100的素数显示出来。 2. 代码 %%  学习目标:寻找素数 clear sum=5;         %求0~100素数之和 ss=0;          %用来标定是否是素数,0表示不是 prime=[2 3];     %用来存放素数,2,

    2024年02月14日
    浏览(32)
  • 41.利用matlab 平衡方程用于图像(matlab程序)

    1. 简述        白平衡 白平衡的英文为White Balance,其基本概念是“不管在任何光源下,都能将白色物体还原为白色”,对在特定光源下拍摄时出现的偏色现象,通过加强对应的补色来进行补偿。 所谓的白平衡是通过对白色被摄物的颜色还原(产生纯白的色彩效果),进而达

    2024年02月14日
    浏览(46)
  • 使用matlab制作音乐

    简谱中最重要的信息就是曲调、节拍,位于简谱左上角, 如图中的1=G,是以G调为基准频率,即1对应G调,其他常见的还有1=C等, 4/4为一节4个1/4拍,一节则为一个短竖线隔开的,相似的还有3/4排,1/2拍等。  对应数字是在以1为基频下的不同音频 0拍表示占位,不出声 数字左上

    2024年02月10日
    浏览(32)
  • 4.利用matlab符号矩阵的四则运算(matlab程序)

    1. 简述     符号对象的建立 sym函数 sym函数用于建立单个符号对象,其常用调用格式为: 符号对象名=sym(A) 1 将由A来建立符号对象,其中,A可以是一个数值常量、数值矩阵或数值表达式(不加单引号),此时符号对象为一个符号常量;A也可以是一个变量名(加单引号),这时符号对

    2024年02月13日
    浏览(46)
  • Matlab根据实验照片制作视频

    关注 M r . m a t e r i a l   , color{Violet} rm Mr.material , Mr.material  

    2023年04月08日
    浏览(31)
  • matlab appdesigner制作UI

    下载工程 图形用户界面(GUI)也称为应用程序,提供对软件应用程序的点击控制,从而无需学习其他语言或键入命令来运行应用程序, 可以打包成matlab内部app或桌面应用. 传统的guide也可以做ui, 但用起来不如appdesigner好用, 已经即将要被淘汰了. 在matlab命令行输入 appdesigner 即可

    2024年02月05日
    浏览(35)
  • Matlab动图保存——GIF制作与视频制作

    在Matlab绘制动图时,若想保存成GIF或视频,可参考以下代码。 (1)GIF格式 (2)视频格式

    2024年02月07日
    浏览(50)
  • MATLAB | 我用MATLAB制作了一款伪3D第一视角迷宫小游戏

    使用键盘上方向键 ↑ 向前移动 使用键盘左右方向键调整 ← → 朝向 原理很简单,如效果图所示,主要就是以角色视角方向发射大量的直线模拟视线,并计算直线与墙壁交点,获取每一条视线方向下,角色到墙壁的距离,然后根据近大远小的原理绘制不同长度的竖向直线模拟

    2024年02月07日
    浏览(53)

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包,优惠每天领

二维码1

领取红包

二维码2

领红包