基於DEM的坡度坡向分析

来源:https://www.cnblogs.com/XiaoRser/archive/2022/05/14/16269417.html
-Advertisement-
Play Games

一、分頁查詢(引用了element框架) 思路: Sql語句:select * from <表名> limit <從第幾條開始查詢>,<查詢多少條數據> 例子:select * from user limit 2,10; 前端:傳遞兩個數據給後端:begin,size 但是我們分頁查詢的頁面肯定不會 ...


坡度坡向分析方法

坡度(slope)是地面特定區域高度變化比率的量度。坡度的表示方法有百分比法、度數法、密位法和分數法四種,其中以百分比法和度數法較為常用。本文計算的為坡度百分比數據。如當角度為45度(弧度為π/4)時,高程增量等於水平增量,高程增量百分比為100%。

 

坡向(aspect)是指地形坡面的朝向。坡向用於識別出從每個像元到其相鄰像元方向上值的變化率最大的下坡方向。坡向可以被視為坡度方向。坡向是一個角度,將按照順時針方向進行測量,角度範圍介於 0(正東)到 360(仍是正東)之間,即完整的圓。不具有下坡方向的平坦區域將賦值為-1(arcgis處理時為-1,其他可能為0)。

坡度、坡向計算一般採用擬合曲面法。擬合曲面一般採用二次曲面,即3×3的視窗,如下圖所示。每個視窗的中心為一個高程點。圖中的中心點e坡度和坡向計算過程如下。

參考鏈接:

[1]https://blog.csdn.net/zhouxuguang236/article/details/40017219

[2]https://blog.csdn.net/weixin_45561357/article/details/106677574

[3]https://www.cnblogs.com/gispathfinder/p/5790469.html

註意:DEM的空間坐標系一定要為投影坐標系

ArcGIS坡度坡向分析

打開DEM數據

坡度分析

 

坡度結果

坡向分析

 

坡向結果

python-gdal坡度坡向分析

from osgeo import gdal

demfile = r"D:\微信公眾號\坡度坡向\N40E117_Albers.tif"

# 獲取DEM信息
infoDEM = gdal.Info(demfile)

# 計算坡度
slopfile = r"D:\微信公眾號\坡度坡向\N40E117_Albers_gdal_Slope.tif"
slope = gdal.DEMProcessing(slopfile, demfile, "slope", format='GTiff', slopeFormat="percent", zeroForFlat=1, computeEdges=True)

# 計算坡向
aspectfile = r"D:\微信公眾號\坡度坡向\N40E117_Albers_gdal_Aspect.tif"
b = gdal.DEMProcessing(aspectfile, demfile, "aspect", format='GTiff', trigonometric=0, zeroForFlat=1, computeEdges=True)

 

坡度結果

 

坡向結果

python坡度坡向分析

import gdal
import numpy as np
from scipy import ndimage as nd
from copy import deepcopy

demfile = r"D:\微信公眾號\坡度坡向\N40E117_Albers.tif"
slopefile = r"D:\微信公眾號\坡度坡向\N40E117_Albers_python_Slope.tif"

#讀取DEM數據
ds = gdal.Open(demfile)
cols = ds.RasterXSize
rows = ds.RasterYSize
geo = ds.GetGeoTransform()
proj = ds.GetProjection()
dem_data = ds.ReadAsArray()
data = deepcopy(dem_data).astype(np.float32)
band = ds.GetRasterBand(1)
nodata = band.GetNoDataValue()
data[data == nodata] = np.nan
# data[data<-999]=np.nan
mask = np.isnan(data)
# 將無效值或背景值臨近像元填充
if np.sum(mask) > 0:
   ind = nd.distance_transform_edt(mask, return_distances=False, return_indices=True)
   data = data[tuple(ind)]

# 計算坡度
xsize = np.abs(geo[1])
ysize = np.abs(geo[5])
x = ((data[:-2, 2:] - data[:-2, :-2]) + 2 * (data[1:-1, 2:] - data[1:-1, :-2]) + (data[2:, 2:] - data[2:, :-2])) / (8 * xsize)
y = ((data[2:, :-2] - data[:-2, :-2]) + 2 * (data[2:, 1:-1] - data[:-2, 1:-1]) + (data[2:, 2:] - data[:-2, 2:])) / (8 * ysize)
s_data = np.full((rows, cols), -999, dtype=np.float32)
s_data[1:-1, 1:-1] = (np.arctan(np.sqrt((np.power(x, 2) + np.power(y, 2)))))
s_data[1:-1, 1:-1] = np.abs(np.tan(s_data[1:-1, 1:-1])) * 100
s_mask = s_data==-999
# 邊緣填充
if np.sum(s_mask) > 0:
   ind = nd.distance_transform_edt(s_mask, return_distances=False, return_indices=True)
   s_data = s_data[tuple(ind)]
# 掩膜
s_data[dem_data==nodata] = -999
# 寫出結果
driver = gdal.GetDriverByName("gtiff")
outds = driver.Create(slopefile, cols, rows, 1, gdal.GDT_Float32)
outds.SetGeoTransform(geo)
outds.SetProjection(proj)
outband = outds.GetRasterBand(1)
outband.WriteArray(s_data)
outband.SetNoDataValue(-999)

 

坡度結果

import gdal
import numpy as np
from scipy import ndimage as nd
from copy import deepcopy

demfile = r"D:\微信公眾號\坡度坡向\N40E117_Albers.tif"
aspectfile = r"D:\微信公眾號\坡度坡向\N40E117_Albers_python_Aspect.tif"

#讀取DEM數據
ds = gdal.Open(demfile)
cols = ds.RasterXSize
rows = ds.RasterYSize
geo = ds.GetGeoTransform()
proj = ds.GetProjection()
dem_data = ds.ReadAsArray()
data = deepcopy(dem_data).astype(np.float32)
band = ds.GetRasterBand(1)
nodata = band.GetNoDataValue()
data[data == nodata] = np.nan
# data[data<-999]=np.nan
mask = np.isnan(data)
# 將無效值或背景值臨近像元填充
if np.sum(mask) > 0:
   ind = nd.distance_transform_edt(mask, return_distances=False, return_indices=True)
   data = data[tuple(ind)]

# 計算坡向
xsize = np.abs(geo[1])
ysize = np.abs(geo[5])
x = ((data[:-2, 2:] - data[:-2, :-2]) + 2 * (data[1:-1, 2:] - data[1:-1, :-2]) + (data[2:, 2:] - data[2:, :-2])) / (8 * xsize)
y = ((data[2:, :-2] - data[:-2, :-2]) + 2 * (data[2:, 1:-1] - data[:-2, 1:-1]) + (data[2:, 2:] - data[:-2, 2:])) / (8 * ysize)
a_data = np.full((rows, cols), -999, dtype=np.float32)
a_data[1:-1, 1:-1] = np.arctan2(y, -1 * x) * 57.29578
a_data_ = deepcopy(a_data[1:-1, 1:-1])
a_data[1:-1, 1:-1][a_data_ < 0] = 90 - a_data[1:-1, 1:-1][a_data_ < 0]
a_data[1:-1, 1:-1][a_data_ >90] = 450 - a_data[1:-1, 1:-1][a_data_ > 90]
a_data[1:-1, 1:-1][(a_data_ >= 0) & (a_data_ <= 90)] = 90 - a_data[1:-1, 1:-1][(a_data_ >= 0) & (a_data_ <= 90)]
a_data[1:-1, 1:-1][(x==0.)& (y==0.)] = -1
a_data[1:-1, 1:-1][(x==0.)& (y>0.)] = 0
a_data[1:-1, 1:-1][(x==0.)& (y<0.)] = 180
a_data[1:-1, 1:-1][(x>0.)& (y==0.)] = 90
a_data[1:-1, 1:-1][(x<0.)& (y==0.)] = 270.

# 邊緣填充
a_mask = a_data==-999
if np.sum(a_mask) > 0:
   ind = nd.distance_transform_edt(a_mask, return_distances=False, return_indices=True)
   a_data = a_data[tuple(ind)]

# 掩膜
a_data[dem_data==nodata] = -999
# 寫出結果
driver = gdal.GetDriverByName("gtiff")
outds = driver.Create(aspectfile, cols, rows, 1, gdal.GDT_Float32)
outds.SetGeoTransform(geo)
outds.SetProjection(proj)
outband = outds.GetRasterBand(1)
outband.WriteArray(a_data)
outband.SetNoDataValue(-999)

 

坡向結果

測試數據:

鏈接:https://pan.baidu.com/s/1PODbTJn1JOqOA4qeaJq4Gg 

提取碼:l3fw 

歡迎關註個人wx_gzh: 小Rser

 

您的分享是我們最大的動力!

-Advertisement-
Play Games
更多相關文章
  • 今天發現之前學的愛前端的課中JS部分函數等不全,果斷換了一個課——渡一的《Web前端開發JavaScript高薪課堂》接著學習,不過廢話有點多 1、條件語句 語法: 1、單if,條件成立,執行語句體 if (條件){ 語句體; } 2、if else,條件成立,執行if後的語句體,否則執行else的 ...
  • JS 頁面演示背景 <!DOCTYPE html> <html lang="en"> <head> <meta charset="UTF-8"> <title>ivanlee</title> <link rel="shortcut icon" href="ab_favicon.ico"> <styl ...
  • 一、概述 Scala是一門多範式的編程語言,一種類似java的編程語言 ,設計初衷是實現可伸縮的語言 、並集成面向對象編程和函數式編程的各種特性。Spark就是使用Scala編寫的。因此為了更好的學習大數據開發, 需要掌握Scala這門語言,當然Spark的興起,也帶動Scala語言的發展!官方文檔 ...
  • 實踐是唯一的真理。 #變數 ##變數的定義 變數就是可以變化的量。 JAVA變數是程式中最基礎的程式單元,其要素包括變數名,變數類型及作用域。 寫程式要註意程式的可讀性 如圖所示,可以一行寫多個對象,但是不建議。 代碼也能使用Ctrl+F搜索,方便排錯。 ##註意事項: 每個變數都有類型,類型可以是 ...
  • GO的環境配置? GOPATH GOROOT 都是幹嘛用的? 配置環境跟java對比有點奇怪 https://blog.csdn.net/weixin_40563757/article/details/115476327 語言特性 協程? 建立一個協程很簡單 加一個go關鍵字就可以 package ...
  • mybatis層編寫完畢後的項目目錄 1.右鍵SpringMVC2項目-》new-》Modual-》選擇maven項目(我的項目名為Study09_ssm),輸入模塊名,點擊Finish 2.第二部的操作就是將idea的基本運行環境搞定,包括:添加web支持,配置tomcat,配置project s ...
  • #類型轉換 由於java是強類型語言,所以在進行某些運算的時候,需要用到類型轉換。 低-->高指的是位元組大小,從小到大。 小數的優先順序大於整數 數值進行類型轉換時不要讓數據溢出 由低到高可以直接轉換,無需額外代碼。 註意點: 1 不能對布爾值進行轉換 2 不能把對象類型轉換為不相干的類型 3 在把高 ...
  • 開發環境: SpringBoot: 2.6.5 SpringCloud: 2021.0.0 SpringCloudAlibaba: 2021.0.1.0 Nacos: 2.1.0 代碼: @Slf4j @Component public class MyInMemoryRouteDefinition ...
一周排行
    -Advertisement-
    Play Games
  • 分組和樹形結構是不一樣的。 樹形結構是以遞歸形式存在。分組是以鍵值對存在的形式,類似於GroupBy這樣的形式。 舉個例子 ID NAME SEX Class 1 張三 男 1 2 李四 女 2 3 王二 男 1 當以Sex為分組依據時則是 Key Value 男 1 張三 男 1 3 王二 男 1 ...
  • NetCore中將SQLServer資料庫備份為Sql腳本 描述: 最近寫項目收到了一個需求, 就是將SQL Server資料庫備份為Sql腳本, 如果是My Sql之類的還好說, 但是在網上搜了一大堆, 全是教你怎麼操作SSMS的, 就很d疼! 解決方案: 通過各種查找資料, 還有一些老哥的幫助, ...
  • 我的Notion Clowd.Squirrel Squirrel.Windows 是一組工具和適用於.Net的庫,用於管理 Desktop Windows 應用程式的安裝和更新。 Squirrel.Windows 對 Windows 應用程式的實現語言沒有任何要求,甚至無需服務端即可完成增量更新。 ...
  • 轉載請註明來源 https://www.cnblogs.com/brucejiao/p/16188865.html 謝謝! 轉載請註明來源 https://www.cnblogs.com/brucejiao/p/16188865.html 謝謝! 轉載請註明來源 https://www.cnblog ...
  • 1. Netty源碼研究筆記(3)——Channel系列 依舊是通過先縱向再橫向的研究方法,在開篇中,我們發現不管是Sever還是Client,最終的啟動是通過調用channel的對應方法來完成的,而這個動作實際在channel綁定的eventLoop中執行。 接下來,我們繼續EchoSever、E ...
  • 大家好,今天給大家介紹一款輕量、快速、穩定可編排的組件式規則引擎框架LiteFlow。 一、LiteFlow的介紹 LiteFlow官方網站和代碼倉庫地址 官方網站:https://yomahub.com/liteflow Gitee托管倉庫:https://gitee.com/dromara/li ...
  • 我使用Spring AOP實現了用戶操作日誌功能 今天答辯完了,復盤了一下系統,發現還是有一些東西值得拿出來和大家分享一下。 需求分析 系統需要對用戶的操作進行記錄,方便未來溯源 首先想到的就是在每個方法中,去實現記錄的邏輯,但是這樣做肯定是不現實的,首先工作量大,其次違背了軟體工程設計原則(開閉原 ...
  • 《零基礎學Java》 繪製幾何圖形 Java可以分別使用 Graphics 和 Graphics2D 繪製圖形,Graphics類 使用不同的方法繪製不同的圖形(drawLine()方法可f以繪製線、drawRect()方法用於繪製矩形、drawOval()方法用於繪製橢圓形)。 Graphics類 ...
  • 本期教程人臉識別第三方平臺為虹軟科技,本文章講解的是人臉識別RGB活體追蹤技術,免費的功能很多可以自行搭配,希望在你看完本章課程有所收穫。 ...
  • 很多人都喜歡使用黑色的主題樣式,包括我自己,使用了差不多三年的黑色主題,但是個人覺得在進行視窗轉換的時候很廢眼睛。 比如IDEA是全黑的,然後需要看PDF或者WORD又變成白色的了,這樣來回切換導致眼睛很累,畢竟現在網頁以及大部分軟體的界面都是白色的。那麼還是老老實實的使用原來比較順眼的模式吧。 1 ...