当前位置: 首页 > news >正文

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

相关文章:

利用MATLAB制作DEM山体阴影

在地理绘图中&#xff0c;我们使用的DEM数据添加山体阴影使得绘制的图件显得更加的美观。 GIS中使用ArcGIS软件就可以达到这一目的&#xff0c;或者使用GMT&#xff0c;同样可以得到山体阴影的效果。 本文提供了一个MATLAB的函数&#xff0c;可以得到山体阴影。 clear all;c…...

ubuntu 使用 rsync 的 SSH 方式同步备份远程WEB服务器

ubuntu 20.04 自带 rsync &#xff0c;对于 WEB 服务器这种更新频率不高的情况&#xff0c;直接使用定时同步复制远程服务器的方法&#xff0c;比较直接和简单&#xff01; $ rsync --version rsync version 3.1.3 protocol version 31 参考&#xff1a; Ubuntu20.04中的rsyn…...

机器学习 | Python实现NARX模型预测控制

机器学习 | Python实现NARX模型预测控制 目录 机器学习 | Python实现NARX模型预测控制效果一览基本介绍研究内容程序设计参考资料效果一览 基本介绍 机器学习 | Python实现NARX模型预测控制 研究内容 贝叶斯黑盒模型预测控制,基于具有外源输入的非线性自回归模型的预期自由能最…...

M5ATOMS3基础03给ROS1发一个问候(rosserial)

引出问题 关于之前2020年的博客&#xff1a; 01. ESP8266和ROS调试一些问题汇总 02. ESP8266和ESP32配置&#xff08;需使用ROS1和ROS2&#xff09; 效果展示 使用M5ATOMS3与ROS1&#xff08;kinetic&#xff0c;melodic&#xff0c;noetic&#xff09;版本通信比较通用的是…...

基于Vue3实现鼠标按下某个元素进行移动,实时改变左侧或右侧元素的宽度,以及点击收起或展开的功能

其原理主要是利用JavaScript中的鼠标事件来控制CSS样式。大致就是监听某个DOM元素的鼠标按下事件&#xff0c;以及按下之后的移动事件和松开事件。在鼠标按下且移动过程中&#xff0c;可实时获得鼠标的X轴坐标的值&#xff0c;通过简单计算&#xff0c;可计算出目标元素的宽度&…...

使用MyBatis(2)

目录 一、定义接口、实体类、创建XML文件实现接口&#xff09; 二、MyBatis的增删改查 &#x1f345;1、MyBatis传递参数查询 &#x1f388;写法一 &#x1f388;写法二 &#x1f388;两种方式的区别 &#x1f345;2、删除操作 &#x1f345;3、根据id修改用户名 &#x…...

【FPGA/D6】

2023年7月25日 VGA控制器 视频23notecodetb 条件编译error时序图保存与读取&#xff1f;&#xff1f;RGBTFT显示屏 视频24PPI未分配的引脚或电平的解决方法 VGA控制器 视频23 note MCU单片机 VGA显示实时采集图像 行消隐/行同步/场同步/场消隐 CRT&#xff1a;阴极射线管 640…...

【WebGIS实例】(10)Cesium开场效果(场景、相机旋转,自定义图片底图)

效果 漫游效果视频&#xff1a; 【WebGIS实例】&#xff08;10&#xff09;Cesium开场效果&#xff08;场景、相机 点击鼠标后将停止旋转并正常加载影像底图&#xff1a; 代码 可以直接看代码&#xff0c;注释写得应该比较清楚了&#xff1a; /** Date: 2023-07-28 16:21…...

【Spring】IOC的原理

一、 IOC 的概念 Spring 的 IOC &#xff0c;即控制反转&#xff0c;所谓控制反转 —— 本来管理业务对象&#xff08;bean&#xff09;的操作是由我们程序员去做的&#xff0c;但是有了 Spring 核心容器后&#xff0c;这些 Bean 对象的创建和管理交给我们Spring容器去做了&am…...

AI加速游戏开发 亚马逊云科技适配3大场景,打造下一代游戏体验

随着疫情的消散&#xff0c;中国游戏产业正在快速前进。在伴随着游戏产业升级的同时&#xff0c;整个行业都在面临着新的挑战与新的诉求。亚马逊云科技游戏研发解决方案和服务&#xff0c;覆盖端到端3大场景&#xff0c;为游戏公司与游戏开发人员赋能。 场景1&#xff1a;AI辅助…...

C++ | 继承(基类,父类,超类),(派生类,子类)

文章参考&#xff1a;https://blog.csdn.net/war1111886/article/details/8609957 一 .继承中的访问权限关系 &#xff11;&#xff0e;基类&#xff0c;父类&#xff0c;超类是指被继承的类&#xff0c;派生类&#xff0c;子类是指继承于基类的类&#xff0e; &#xff12;…...

Commands Of Hadoop

序言 持续整理下常用的命令cuiyaonan2000163.com Command 文件拷贝 当从多个源拷贝时&#xff0c;如果两个源冲突&#xff0c;distcp会停止拷贝并提示出错信息&#xff0c;. 如果在目的位置发生冲突&#xff0c;会根据选项设置解决。 默认情况会跳过已经存在的目标文件&am…...

SQL-每日一题【620.有趣的电影】

题目 某城市开了一家新的电影院&#xff0c;吸引了很多人过来看电影。该电影院特别注意用户体验&#xff0c;专门有个 LED显示板做电影推荐&#xff0c;上面公布着影评和相关电影描述。 作为该电影院的信息部主管&#xff0c;您需要编写一个 SQL查询&#xff0c;找出所有影片…...

linux 精华总结

...

Eureka 学习笔记2:客户端 DiscoveryClient

版本 awsVersion ‘1.11.277’ DiscoveryClient # cacheRefreshTask // 配置shouldFetchRegistry if (clientConfig.shouldFetchRegistry()) {// 配置client.refresh.intervalint registryFetchIntervalSeconds clientConfig.getRegistryFetchIntervalSeconds();// 配置expB…...

okhttp原理分析

工程目录图 请点击下面工程名称&#xff0c;跳转到代码的仓库页面&#xff0c;将工程 下载下来 Demo Code 里有详细的注释 01okhttp module里 包含的设计模式&#xff1a;建造者设计模式、责任链设计模式 CustomInject 演示自定义注解 代码&#xff1a;okhttp原理分析、Andro…...

freeswitch的mod_xml_curl模块

概述 freeswitch是一款简单好用的VOIP开源软交换平台。 随着fs服务的增多&#xff0c;每一台fs都需要在后台单独配置&#xff0c;耗时耗力&#xff0c;心力憔悴。 如果有一个集中管理配置的配置中心&#xff0c;统一管理所有fs的配置&#xff0c;并可以实现动态的修改配置就…...

高速数据采集专家-FMC140【产品手册】

FMC140是一款具有缓冲模拟输入的低功耗、12位、双通道&#xff08;5.2GSPS/通道&#xff09;、单通道10.4GSPS、射频采样ADC模块&#xff0c;该板卡为FMC标准&#xff0c;符合VITA57.1规范&#xff0c;该模块可以作为一个理想的IO单元耦合至FPGA前端&#xff0c;8通道的JESD204…...

【SSM】知识集锦

项目一&#xff1a;狂神JAVA 功能1&#xff1a;实现全部书籍查询 1.思路&#xff1a;首页index.jsp ——>Controller——>hello.jsp 2.步骤&#xff1a; step1:index.jsp <% page language"java" contentType"text/html; charsetUTF-8" page…...

Flowable-中间事件-信号中间抛出事件

定义 当流程执行到达信号抛出事件时&#xff0c;流程引擎会直接抛出信号&#xff0c;其他引用了与其相同的信号捕获 事件会被触发&#xff0c;信号发出后事件结束&#xff0c;流程沿后继路线继续执行。其抛出的信号可以被信号开始事 件&#xff08;Signal Start Event&#xf…...

从WWDC看苹果产品发展的规律

WWDC 是苹果公司一年一度面向全球开发者的盛会&#xff0c;其主题演讲展现了苹果在产品设计、技术路线、用户体验和生态系统构建上的核心理念与演进脉络。我们借助 ChatGPT Deep Research 工具&#xff0c;对过去十年 WWDC 主题演讲内容进行了系统化分析&#xff0c;形成了这份…...

Java 8 Stream API 入门到实践详解

一、告别 for 循环&#xff01; 传统痛点&#xff1a; Java 8 之前&#xff0c;集合操作离不开冗长的 for 循环和匿名类。例如&#xff0c;过滤列表中的偶数&#xff1a; List<Integer> list Arrays.asList(1, 2, 3, 4, 5); List<Integer> evens new ArrayList…...

在rocky linux 9.5上在线安装 docker

前面是指南&#xff0c;后面是日志 sudo dnf config-manager --add-repo https://download.docker.com/linux/centos/docker-ce.repo sudo dnf install docker-ce docker-ce-cli containerd.io -y docker version sudo systemctl start docker sudo systemctl status docker …...

【磁盘】每天掌握一个Linux命令 - iostat

目录 【磁盘】每天掌握一个Linux命令 - iostat工具概述安装方式核心功能基础用法进阶操作实战案例面试题场景生产场景 注意事项 【磁盘】每天掌握一个Linux命令 - iostat 工具概述 iostat&#xff08;I/O Statistics&#xff09;是Linux系统下用于监视系统输入输出设备和CPU使…...

鸿蒙中用HarmonyOS SDK应用服务 HarmonyOS5开发一个医院挂号小程序

一、开发准备 ​​环境搭建​​&#xff1a; 安装DevEco Studio 3.0或更高版本配置HarmonyOS SDK申请开发者账号 ​​项目创建​​&#xff1a; File > New > Create Project > Application (选择"Empty Ability") 二、核心功能实现 1. 医院科室展示 /…...

在四层代理中还原真实客户端ngx_stream_realip_module

一、模块原理与价值 PROXY Protocol 回溯 第三方负载均衡&#xff08;如 HAProxy、AWS NLB、阿里 SLB&#xff09;发起上游连接时&#xff0c;将真实客户端 IP/Port 写入 PROXY Protocol v1/v2 头。Stream 层接收到头部后&#xff0c;ngx_stream_realip_module 从中提取原始信息…...

现代密码学 | 椭圆曲线密码学—附py代码

Elliptic Curve Cryptography 椭圆曲线密码学&#xff08;ECC&#xff09;是一种基于有限域上椭圆曲线数学特性的公钥加密技术。其核心原理涉及椭圆曲线的代数性质、离散对数问题以及有限域上的运算。 椭圆曲线密码学是多种数字签名算法的基础&#xff0c;例如椭圆曲线数字签…...

智能分布式爬虫的数据处理流水线优化:基于深度强化学习的数据质量控制

在数字化浪潮席卷全球的今天&#xff0c;数据已成为企业和研究机构的核心资产。智能分布式爬虫作为高效的数据采集工具&#xff0c;在大规模数据获取中发挥着关键作用。然而&#xff0c;传统的数据处理流水线在面对复杂多变的网络环境和海量异构数据时&#xff0c;常出现数据质…...

Java多线程实现之Thread类深度解析

Java多线程实现之Thread类深度解析 一、多线程基础概念1.1 什么是线程1.2 多线程的优势1.3 Java多线程模型 二、Thread类的基本结构与构造函数2.1 Thread类的继承关系2.2 构造函数 三、创建和启动线程3.1 继承Thread类创建线程3.2 实现Runnable接口创建线程 四、Thread类的核心…...

论文笔记——相干体技术在裂缝预测中的应用研究

目录 相关地震知识补充地震数据的认识地震几何属性 相干体算法定义基本原理第一代相干体技术&#xff1a;基于互相关的相干体技术&#xff08;Correlation&#xff09;第二代相干体技术&#xff1a;基于相似的相干体技术&#xff08;Semblance&#xff09;基于多道相似的相干体…...