python - 使用 gdal 将数组转换为 tiff 光栅图像

Question

更新

我尝试关注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) 

任何想法

Answer 1

该错误准确地说明了问题所在。您正在尝试将 int 类型乘以 NoneType(无)。最可能的情况是 nodataval 为 None，这是因为文件名的第一个栅格波段的 NoDataValue 未定义。您的 print(nodataval) 命令应该证明这一点。 请记住，打印 None 不会显示为字符串“None”。它显示为空白或无字符。

您的编辑显示 nodataval 为 None。