更新资料
我尝试遵循this tutorial
但是我不知道如何使用GDAL将坡度/纵横比导出到新的tiff图像？
完整代码：
from __future__ import division
from osgeo import gdal
from matplotlib.colors import ListedColormap
from matplotlib import colors
import sys
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import math
filename = 'dem.tif'
def getResolution(rasterfn):
raster = gdal.Open(rasterfn)
geotransform = raster.GetGeoTransform()
res = {"east-west": abs(geotransform[1]),
"north-south": abs(geotransform[5])}
return res
def raster2array(rasterfn):
raster = gdal.Open(rasterfn)
band = raster.GetRasterBand(1)
return band.ReadAsArray()
def getNoDataValue(rasterfn):
raster = gdal.Open(rasterfn)
band = raster.GetRasterBand(1)
return band.GetNoDataValue()
data_array = raster2array(filename)
nodataval = getNoDataValue(filename)
resolution = getResolution(filename)
print(resolution)
print(nodataval)
print(type(data_array))
print(data_array.shape)
num_rows = data_array.shape[0]
num_cols = data_array.shape[1]
slope_array = np.ones_like(data_array) * nodataval
aspect_array = np.ones_like(data_array) * nodataval
for i in range(1, num_rows - 1):
for j in range(1, num_cols - 1):
a = data_array[i - 1][j - 1]
b = data_array[i - 1][j]
c = data_array[i - 1][j + 1]
d = data_array[i][j - 1]
e = data_array[i][j]
f = data_array[i][j + 1]
g = data_array[i + 1][j - 1]
h = data_array[i + 1][j]
q = data_array[i + 1][j + 1]
vals = [a, b, c, d, e, f, g, h, q]
if nodataval in vals:
all_present = False
else:
all_present = True
if all_present == True:
dz_dx = (c + (2 * f) + q - a - (2 * d) - g) / (8 * resolution['east-west'])
dz_dy = (g + (2 * h) + q - a - (2 * b) - c) / (8 * resolution['north-south'])
dz_dx_sq = math.pow(dz_dx, 2)
dz_dy_sq = math.pow(dz_dy, 2)
rise_run = math.sqrt(dz_dx_sq + dz_dy_sq)
slope_array[i][j] = math.atan(rise_run) * 57.29578
aspect = math.atan2(dz_dy, (-1 * dz_dx)) * 57.29578
if aspect < 0:
aspect_array[i][j] = 90 - aspect
elif aspect > 90:
aspect_array[i][j] = 360 - aspect + 90
else:
aspect_array[i][j] = 90 - aspect
hist, bins = np.histogram(slope_array, bins=100, range=(0, np.amax(slope_array)))
width = 0.7 * (bins[1] - bins[0])
center = (bins[:-1] + bins[1:]) / 2
plt.bar(center, hist, align='center', width=width)
plt.xlabel('Slope (degrees)')
plt.ylabel('Frequency')
plt.show()
color_map = ListedColormap(['white', 'darkgreen', 'green', 'limegreen', 'lime',
'greenyellow', 'yellow', 'gold',
'orange', 'orangered', 'red'])
# range begins at negative value so that missing values are white
color_bounds = list(range(-3, math.ceil(np.amax(slope_array)), 1))
color_norm = colors.BoundaryNorm(color_bounds, color_map.N)
#Create the plot and colorbar
img = plt.imshow(slope_array, cmap = color_map, norm = color_norm)
cbar = plt.colorbar(img, cmap = color_map, norm = color_norm,
boundaries = color_bounds, ticks = color_bounds)
#Show the visualization
plt.axis('off')
plt.title("Slope (degrees)")
plt.show()
plt.close()
我要按照此tutorial进行首次出口试验，但是坡度不正确
拥有比原学位更低的学位（使用gis程序）
import gdal, ogr, os, osr
import numpy as np
def array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array):
cols = array.shape[1]
rows = array.shape[0]
originX = rasterOrigin[0]
originY = rasterOrigin[1]
driver = gdal.GetDriverByName('GTiff')
outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Byte)
outRaster.SetGeoTransform((originX, pixelWidth, 0, originY, 0, pixelHeight))
outband = outRaster.GetRasterBand(1)
outband.WriteArray(array)
outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromEPSG(4326)
outRaster.SetProjection(outRasterSRS.ExportToWkt())
outband.FlushCache()
def main(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array):
reversed_arr = slope_array # reverse array so the tif looks like the array
array2raster(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,reversed_arr) # convert array to raster
if __name__ == "__main__":
rasterOrigin = (-123.25745,45.43013)
pixelWidth = 10
pixelHeight = 10
newRasterfn = 'test.tif'
array = slope_array
main(newRasterfn,rasterOrigin,pixelWidth,pixelHeight,array)
第二次试用，我遵循此quest，但我不进行任何导出
def array_to_raster(slope_array):
"""Array > Raster
Save a raster from a C order array.
:param array: ndarray
"""
dst_filename = 'xxx.tiff'
x_pixels = num_rows # number of pixels in x
y_pixels = num_cols # number of pixels in y
driver = gdal.GetDriverByName('GTiff')
dataset = driver.Create(
dst_filename,
x_pixels,
y_pixels,
1,
gdal.GDT_Float32, )
dataset.GetRasterBand(1).WriteArray(array)
dataset.FlushCache() # Write to disk.
return dataset, dataset.GetRasterBand(1)
任何想法
最佳答案
该错误确切说明了问题所在。您正在尝试将int类型乘以NoneType（无）。最可能的情况是nodataval为None，这是因为未定义文件名的第一个栅格波段的NoDataValue。您的print（nodataval）命令应证明这一点。请记住，打印无不会显示为字符串“无”。它显示为空白或无字符。
您的编辑显示nodataval为None。
关于python - 阵列以gdal转换光栅图像，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/39835794/