利用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

其中调用的函数 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 givenindi=strcmpi(varargin{i},opts);ind=find(indi==1);if isempty(ind)~=1switch indcase 1azimuth=varargin{i+1};case 2altitude=varargin{i+1};case 3zf=varargin{i+1};case 4plotit=1;endendend
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-1for j=2:jm-1fx(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)~=0asp(i,j) = atan2(fy(i,j),-fx(i,j));if asp(i,j)<0asp(i,j)=asp(i,j)+2*pi;endelseif fy(i,j)>0asp(i,j)=pi/2;elseif fy(i,j)<0asp(i,j)=3*pi/2;endend     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==1figureimagesc(X,Y,h')axis imageset(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:

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

 

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

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mzph.cn/news/11915.shtml

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈email:809451989@qq.com,一经查实,立即删除!

相关文章

C#中窗体之间的传值操作(传递值和获取返回值)

在C#的应用程序开发中&#xff0c;窗体之间的传值操作是不可避免会经常遇到的操作。   比如&#xff0c;在一个窗体中调用另外一个窗体、传递值并且得到返回值&#xff0c;通常情况下有哪些方法呢&#xff1f;   一般情况下&#xff0c;通过工程项目的公有变量、窗体的公有…

React 中 ref 如何使用?

在React 中要使用 ref &#xff0c;首先要创建一个新的对象 // 创建 ref 对象于jsx 绑定const inputRer useRef(null);在使用 ref 时不要在组件渲染时使用 ref 进行 dom 操作&#xff0c;因为此时ref 还没有值&#xff0c;会报错 &#xff08;注意&#xff1a;操作文本框时尽…

AI面试官:LINQ和Lambda表达式(二)

AI面试官&#xff1a;LINQ和Lambda表达式&#xff08;二&#xff09; 当面试官面对C#中关于LINQ和Lambda表达式的面试题时&#xff0c;通常会涉及这两个主题的基本概念、用法、实际应用以及与其他相关技术的对比等。以下是一些可能的面试题目&#xff0c;附带简要解答和相关案…

解放Linux内存:释放缓存(linux释放缓存)

随着软件越来越复杂&#xff0c;内存变得越来越宝贵。尤其是在Linux系统上&#xff0c;内存管理策略十分重要。它不仅可以帮助系统保持高效运行&#xff0c;而且也能够让程序有更多的空间来运行&#xff0c;避免系统出现假死和其他性能问题。 在Linux系统中&#xff0c;释放缓…

《面试1v1》如何能从Kafka得到准确的信息

&#x1f345; 作者简介&#xff1a;王哥&#xff0c;CSDN2022博客总榜Top100&#x1f3c6;、博客专家&#x1f4aa; &#x1f345; 技术交流&#xff1a;定期更新Java硬核干货&#xff0c;不定期送书活动 &#x1f345; 王哥多年工作总结&#xff1a;Java学习路线总结&#xf…

安防视频管理平台GB设备接入EasyCVR, 如何获取RTMP与RTSP视频流

安防视频监控平台EasyCVR可拓展性强、视频能力灵活、部署轻快&#xff0c;可支持的主流标准协议有国标GB28181、RTSP/Onvif、RTMP等&#xff0c;以及支持厂家私有协议与SDK接入&#xff0c;包括海康Ehome、海大宇等设备的SDK等。平台既具备传统安防视频监控的能力&#xff0c;比…

基于粒子群优化算法的分布式电源选址与定容【多目标优化】【IEEE33节点】(Matlab代码实现)

目录 &#x1f4a5;1 概述 1.1 目标函数 2.2 约束条件 &#x1f4da;2 运行结果 &#x1f389;3 参考文献 &#x1f308;4 Matlab代码实现 &#x1f4a5;1 概述 分布式电源接入配电网&#xff0c;实现就地消纳&#xff0c;可以提高新能源的利用率、提高电能质量和降低系统网损…

网络渗透入门指南

目录 简介&#xff1a; 1. 什么是网络渗透&#xff1f; 2. 渗透测试类型 1.黑盒测试 2.白盒测试 3.灰盒测试 3. 渗透测试步骤 3.1 信息收集 3.2 漏洞扫描与评估 3.3 漏洞利用 3.4 特权升级与持久性访问 3.5 横向移动 3.6 数据获取与报告 实例&#xff1a; 总结 简…

ResponseEntity

ResponseEntity : 是 spring framework中的一个类&#xff0c;用于封装 http响应的实体部分&#xff0c;包含&#xff0c;主题内容&#xff0c;http响应头&#xff0c;http状态码&#xff0c;等信息&#xff0c;&#xff0c; 用来返回相应给客户端 HttpStatus 是个枚举类&…

出海周报|Temu在美状告shein、ChatGPT安卓版上线、小红书回应闪退

工程机械产业“出海”成绩喜人&#xff0c;山东相关企业全国最多Temu在美状告shein&#xff0c;跨境电商战事升级TikTok将在美国推出电子商务计划&#xff0c;售卖中国商品高德即将上线国际图服务&#xff0c;初期即可覆盖全球超200个国家和地区ChatGPT安卓版正式上线&#xff…

echarts遇到的问题

文章目录 折线图-区域面积图 areaStyley轴只有整数y轴不从0开始y轴数值不确定&#xff0c;有大有小&#xff0c;需要动态处理折线-显示label标线legend的格式化和默认选中状态x轴的lable超长处理x轴的相关设置 echarts各个场景遇到的问题 折线图-区域面积图 areaStyle areaStyl…

node.js的优点

提示&#xff1a;node.js的优点 文章目录 一、什么是node.js二、node.js的特性 一、什么是node.js 提示&#xff1a;什么是node.js? Node.js发布于2009年5月&#xff0c;由Ryan Dahl开发&#xff0c;是一个基于ChromeV8引擎的JavaScript运行环境&#xff0c;使用了一个事件驱…

Redis 场景

1.分布式锁实现 上锁&#xff1a;setnx key value 解锁&#xff1a;del key 超时机制&#xff1a;set key value nx ex time 基于此思路实现伪代码&#xff1a; public class RedisLock{Autowiredprivate static RedisTemplate rt; private static final LOCK_WORD &quo…

【c语言进阶】字符函数和字符串函数知识总结

字符函数和字符串函数 前期背景求字符串长度函数strlen函数strlen函数三种模拟实现 长度不受限制的字符串函数strcpy函数strcpy函数模拟实现strcat函数strcat函数模拟实现strcmp函数strcmp函数模拟实现 长度受限制的字符串函数strncpy函数strncpy函数模拟实现strncat函数strnca…

粘包处理的方式

为什么出现粘包&#xff1a; 发送端在发送的时候由于 Nagel 算法的存在会将字节数较小的数据整合到一起发送&#xff0c;导致粘包&#xff1b;接收端不知道发送端数据的长度&#xff0c;导致接收时无法区分数据&#xff1b; 粘包处理的方式&#xff1a; 通过在数据前面加上报…

最新版本docker 设置国内镜像源 加速办法

解决问题:加速 docker 设置国内镜像源 目录: 国内加速地址 修改方法 国内加速地址 1.Docker中国区官方镜像 https://registry.docker-cn.com 2.网易 http://hub-mirror.c.163.com 3.ustc https://docker.mirrors.ustc.edu.cn 4.中国科技大学 https://docker.mirrors…

windows11打不开任务管理器,

目录 第一章、win11系统任务管理器打不开&#xff1f;第二章、解决方式修改注册表 友情提醒&#xff1a; 先看文章目录&#xff0c;大致了解文章知识点结构&#xff0c;点击文章目录可直接跳转到文章指定位置。 第一章、win11系统任务管理器打不开&#xff1f; Win11任务管理…

Docker 安全及日志管理与https部署

容器的安全性问题的根源在于容器和宿主机共享内核。如果容器里的应用导致Linux内核崩溃&#xff0c;那么整个系统可能都会崩溃。与虚拟机是不同的&#xff0c;虚拟机并没有与主机共享内核&#xff0c;虚拟机崩溃一般不会导致宿主机崩溃。 Docker 容器与虚拟机的区别 虚拟机通…

FPGA+EMMC 8通道存储小板

FPGA 采用XILINX公司A7100作为主芯片 AD采用AD7606及一款陀螺仪传感器&#xff0c;可以实时存储到EMMC&#xff0c;系统分为采集模式及回放模式 通过232接口对工作模式进行配置&#xff0c;采样率可以动态配置 回放采用W5100S通过TCP协议进行回放数据

二、使用运行自己的docker python容器环境

第一篇参考&#xff1a; https://blog.csdn.net/weixin_42357472/article/details/131953866 运行容器同时执行命令或脚本 1&#xff09;这是打开一个对外的jupyter notebook容器环境 docker run -d --name my_container -p 8090:8888 mynewpythonimage jupyter notebook --…