Python计算坡度坡向并输出二维、三维图
作者：互联网
在上次代码的基础上做了一点儿修改，将定义的函数单独放在一个模块里面，主函数去单独调用该模块。
DEMslopeAspect模块
from osgeo import gdal,ogr,osr
import numpy as np
import math
import datetime
# Python matplotlib模块代码示例 https://vimsky.com/examples/detail/python-module-matplotlib.html
# Axes3D是mpl_toolkits.mplot3d中的一个绘图函数，mpl_toolkits.mplot3d;是Matplotlib里面专门用来画三维图的工具包。
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cbook, cm
from matplotlib.colors import LightSource
import matplotlib.pyplot as plt
# 正则表达式(regular expression)描述了一种字符串匹配的模式，可以用来检查一个串是否含有某种子串、将匹配的子串做替换或者从某个串中取出符合某个条件的子串等。
import re
# 读取TIFF遥感影像
def read_img(filename):
# dataset = gdal.Open(r'D:\ProfessionalProfile\DEMdata\slopeAspectPython0322\test2_3-0326\test.tif') # 打开文件
# dataset = gdal.Open(r'D:\ProfessionalProfile\DEMdata\slopeAspectPython0322\proUTM50python.tif')
dataset = gdal.Open(filename)
im_width = dataset.RasterXSize # 栅格矩阵的列数
im_height = dataset.RasterYSize # 栅格矩阵的行数
im_bands = dataset.RasterCount # 波段数
im_geotrans = dataset.GetGeoTransform() # 仿射矩阵，左上角像素的大地坐标和像素分辨率
im_proj = dataset.GetProjection() # 地图投影信息，字符串表示
im_data = dataset.ReadAsArray(0, 0, im_width, im_height)
datatype = im_data.dtype
del dataset # 关闭对象dataset，释放内存
return im_data, im_proj, im_geotrans, im_height,im_width, im_bands, datatype
# 为便于后续坡度计算，需要在原图像的周围添加一圈数值
def AddRound(npgrid):
nx, ny = npgrid.shape[0], npgrid.shape[1] # ny:行数，nx:列数;此处注意顺序
# np.zeros()返回来一个给定形状和类型的用0填充的数组；
zbc=np.zeros((nx+2,ny+2))
# 填充原数据数组
zbc[1:-1,1:-1]=npgrid
#四边填充数据
zbc[0,1:-1]=npgrid[0,:] #上边；0行，所有列；
zbc[-1,1:-1]=npgrid[-1,:] #下边；最后一行，所有列；
zbc[1:-1,0]=npgrid[:,0] #左边；所有行，0列。
zbc[1:-1,-1]=npgrid[:,-1] #右边；所有行，最后一列
#填充剩下四个角点值
zbc[0,0]=npgrid[0,0]
zbc[0,-1]=npgrid[0,-1]
zbc[-1,0]=npgrid[-1,0]
zbc[-1,-1]=npgrid[-1,0]
return zbc
#####计算xy方向的梯度
def Cacdxdy(npgrid,sizex,sizey):
nx, ny = npgrid.shape
s_dx = np.zeros((nx,ny))
s_dy = np.zeros((nx,ny))
a_dx = np.zeros((nx, ny))
a_dy = np.zeros((nx, ny))
# 忘记加range报错：object is not iterable
# 坡度、坡向变化率的计算：https://help.arcgis.com/zh-cn/arcgisdesktop/10.0/help/index.html#/na/009z000000vz000000/
for i in range(1,nx-1):
for j in range(1,ny-1):
s_dx[i,j] = ((npgrid[i-1,j+1]+2*npgrid[i,j+1]+npgrid[i+1,j+1])-(npgrid[i-1,j-1]+2*npgrid[i,j-1]+npgrid[i+1,j-1])) / (8 * sizex)
s_dy[i, j] = ((npgrid[i+1, j-1] + 2 * npgrid[i+1, j] + npgrid[i+1,j+1])-(npgrid[i-1,j-1]+2 * npgrid[i-1,j] + npgrid[i-1,j+1])) / (8 * sizey)
a_dx=s_dx*sizex
a_dy=s_dy*sizey
# 保留原数据区域的梯度值
s_dx = s_dx[1:-1,1:-1]
s_dy = s_dy[1:-1,1:-1]
a_dx = a_dx[1:-1, 1:-1]
a_dy = a_dy[1:-1, 1:-1]
# np.savetxt(r"D:\ProfessionalProfile\DEMdata\slopeAspectPython0322\1dxdy.csv",dx,delimiter=",")
return s_dx,s_dy,a_dx,a_dy
####计算坡度/坡向
def CacSlopAsp(s_dx,s_dy,a_dx,a_dy):
# 坡度
slope=(np.arctan(np.sqrt(s_dx*s_dx+s_dy*s_dy)))*180/math.pi #转换成°
#坡向
# #出错：TypeError: only size-1 arrays can be converted to Python scalars
# a2 = math.atan2(a_dy,-a_dx)*180/math.pi
a=np.zeros((a_dy.shape[0],a_dy.shape[1]))
for i in range(0,a_dx.shape[0]):
for j in range(0,a_dx.shape[1]):
a[i,j] = math.atan2(a_dy[i,j], -a_dx[i,j]) * 180 / math.pi
# 输出
aspect = a
# 坡向值将根据以下规则转换为罗盘方向值（0 到 360 度）：
# https://help.arcgis.com/zh-cn/arcgisdesktop/10.0/help/index.html#/na/009z000000vp000000/
x, y = a.shape[0],a.shape[1]
for m in range(0,x):
for n in range(0,y):
if a[m,n] < 0:
aspect[m,n] = 90-a[m,n]
elif a[m,n] > 90:
aspect[m,n] = 360.0 - a[m,n] + 90.0
else:
aspect[m,n] = 90.0 - a[m,n]
return slope,aspect
# 遥感影像的存储，写GeoTiff文件
def write_img(filename, tar_proj, im_geotrans, im_data, datatype):
# 判断栅格数据的数据类型
if 'int8' in im_data.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in im_data.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
# 判读数组维数
if len(im_data.shape) == 3:
# 注意数据的存储波段顺序：im_bands, im_height, im_width
im_bands, im_height, im_width = im_data.shape
else:
im_bands, (im_height, im_width) = 1, im_data.shape
# 创建文件时 driver = gdal.GetDriverByName("GTiff")，数据类型必须要指定，因为要计算需要多大内存空间。
driver = gdal.GetDriverByName("GTiff")
dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)
dataset.SetGeoTransform(im_geotrans) # 写入仿射变换参数
dataset.SetProjection(tar_proj) # 写入投影
if im_bands == 1:
dataset.GetRasterBand(1).WriteArray(im_data) # 写入数组数据
else:
for i in range(im_bands):
dataset.GetRasterBand(i + 1).WriteArray(im_data[i])
del dataset
# 定义投影函数（此次运行没有用到）
def SetPro(filename,tar_proj,outputfilename):
ds = gdal.Open(filename)
im_geotrans = ds.GetGeoTransform() # 仿射矩阵信息
im_proj = ds.GetProjection() # 地图投影信息
im_width = ds.RasterXSize # 栅格矩阵的列数
im_height = ds.RasterYSize # 栅格矩阵的行数
im_bands = ds.RasterCount
ds_array = ds.ReadAsArray(0, 0, im_width, im_height) # 获取原数据信息，包括数据类型int16，维度，数组等信息
# 设置数据类型(原图像有负值)
datatype = gdal.GDT_Float32
# 目标投影
img_proj = tar_proj
# 输出影像路径及名称
name = outputfilename
driver = gdal.GetDriverByName("GTiff") # 创建文件驱动
dataset = driver.Create(name, im_width, im_height, im_bands, datatype)
dataset.SetGeoTransform(im_geotrans) # 写入原图像的仿射变换参数
dataset.SetProjection(img_proj) # 写入目标投影
# 写入影像数据
dataset.GetRasterBand(1).WriteArray(ds_array)
del dataset
####绘制平面栅格图
def Drawgrid(judge,pre=[],A=[],strs=""):
if judge==0:
if strs == "":
plt.imshow(A, interpolation='nearest', cmap=plt.cm.hot, origin='lower') # cmap='bone' cmap=plt.cm.hot
# plt.imshow(A, interpolation='nearest', cmap=plt.cm.hot, origin='lower') # cmap='bone' cmap=plt.cm.hot
plt.colorbar(shrink=0.8)
plt.xticks(())
plt.yticks(())
plt.show()
else:
plt.imshow(A, interpolation='nearest', cmap=strs, origin='lower') # cmap='bone' cmap=plt.cm.hot
plt.colorbar(shrink=0.8)
# 影像范围（原始图像的im_geotrans六参数有）
X = np.arange(113.99986111111112, 6113.999861111111, 30)
Y = np.arange(35.00013888888889, 3035.000138888889, 30)
plt.xticks(())
plt.yticks(())
plt.show()
# judge==1绘制三维DEM
elif judge==1:
fig = plt.figure()
ax = Axes3D(fig)
X = np.arange(113.99986111111112,6113.999861111111, 30)
Y = np.arange(35.00013888888889, 3035.000138888889, 30)
# xt=range(114.79763889,853584.79763889 , 30)
# yt=range(38.21347222, 413348.21347222, 30)
X, Y = np.meshgrid(X, Y)
Z = pre
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=plt.get_cmap('rainbow')) # cmap=plt.get_cmap('rainbow')
ax.contourf(X, Y, Z, zdir='z', offset=-2, cmap=plt.cm.hot)
ax.set_zlim(0, 200000)
plt.show()
主函数
import D1_DEMslopeAspect as dem
from D1_DEMslopeAspect import Drawgrid
import datetime
# 程序入口
if __name__ == "__main__":
startime = datetime.datetime.now() # 程序开始时间
# 读取ASTER GDEM遥感影像
demgrid, proj, geotrans, row, column, band, type =dem.read_img(r'D:\ProfessionalProfile\DEMdata\slopeAspectPython0322\test2_3-0326\test2.tif')
# geotrans = (114.79763889, 0.00027777777778, 0.0, 38.21347222, 0.0, -0.00027777777778)
# row = 13777
# col = 28449
demgridata = demgrid
# 为计算梯度给影像添加周围一圈数据
demgrid = dem.AddRound(demgrid)
# 梯度计算
dx1,dy1,dx2,dy2 = dem.Cacdxdy(demgrid,30,30)
# 坡度、坡向计算
slope,aspect =dem.CacSlopAsp(dx1,dy1,dx2,dy2)
# 设置要投影的投影信息，此处是WGS84-UTM-50N
tar_proj = '''PROJCS["WGS 84 / UTM zone 50N",
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.01745329251994328,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4326"]],
UNIT["metre",1,
AUTHORITY["EPSG","9001"]],
PROJECTION["Transverse_Mercator"],
PARAMETER["latitude_of_origin",0],
PARAMETER["central_meridian",117],
PARAMETER["scale_factor",0.9996],
PARAMETER["false_easting",500000],
PARAMETER["false_northing",0],
AUTHORITY["EPSG","32650"],
AXIS["Easting",EAST],
AXIS["Northing",NORTH]]'''
# 输出TIFF格式遥感影像,并设置投影坐标
slopeT = dem.write_img(r'D:\ProfessionalProfile\DEMdata\slopeAspectPython0322\test2_3-0326\slopetest2.tif', tar_proj, geotrans, slope, type)
aspectT = dem.write_img(r'D:\ProfessionalProfile\DEMdata\slopeAspectPython0322\test2_3-0326\aspecttest2.tif', tar_proj, geotrans, aspect, type)
endtime = datetime.datetime.now() # 程序结束时间
runtime = endtime - startime # 程序运行时间
print('运行时间为： %d 秒' % (runtime.seconds))
# 绘制三维DEM
Drawgrid(judge=1, pre=demgridata)
# 绘制二维DEM
Drawgrid(judge=0, A=demgridata, strs="bone")
# 绘制坡度图
Drawgrid(judge=0, A=slope, strs="rainbow")
# 绘制坡向图
Drawgrid(judge=0, A=aspect)
效果图
由于整个山东省的面积太大，故而选择了一小片区域（100*200）测试。
原始DEM图
TIF坡度图
TIF坡向图
三维DEM图
二维DEM图
坡度图
坡向图
参考
[1]坡度：https://help.arcgis.com/zh-cn/arcgisdesktop/10.0/help/index.html#/na/009z000000vz000000/
[2]坡向：https://help.arcgis.com/zh-cn/arcgisdesktop/10.0/help/index.html#/na/009z000000vp000000/
[3]博主锃光瓦亮的枕小路：https://blog.csdn.net/weixin_45561357/article/details/106677574
[4]https://blog.csdn.net/weixin_40501429/article/details/114894497
[5]博主箜_Kong：https://blog.csdn.net/liminlu0314/article/details/8498985?ops_request_misc=%257B%2522request%255Fid%2522%253A%2522161657597316780266219174%2522%252C%2522scm%2522%253A%252220140713.130102334.pc%255Fblog.%2522%257D&request_id=161657597316780266219174&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~blog~first_rank_v1~rank_blog_v1-1-8498985.pc_v1_rank_blog_v1&utm_term=%E5%9D%A1%E5%BA%A6%E5%9D%A1%E5%90%91
标签：坡向,plt,Python,dataset,三维图,npgrid,dx,dy,im 来源： https://blog.csdn.net/weixin_40501429/article/details/115302922