在上次代码的基础上做了一点儿修改,将定义的函数单独放在一个模块里面,主函数去单独调用该模块。
DEMslopeAspect模块from osgeo import gdal,ogr,osrimport numpy as npimport mathimport 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 Axes3Dfrom matplotlib import cbook, cmfrom matplotlib.colors import lightSourceimport 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\slopeAspectPython0322dxdy.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 demfrom D1_DEMslopeAspect import DrawgrIDimport 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
总结以上是内存溢出为你收集整理的Python计算坡度坡向并输出二维、三维图全部内容,希望文章能够帮你解决Python计算坡度坡向并输出二维、三维图所遇到的程序开发问题。
如果觉得内存溢出网站内容还不错,欢迎将内存溢出网站推荐给程序员好友。
欢迎分享,转载请注明来源:内存溢出
评论列表(0条)