array to tiff raster image with gdal
<p>Update </p>

<p>i try to follow <a href="https://earthlab.github.io/python/DEM-slope-aspect-python/" rel="nofollow">this tutorial</a> </p>

<p>but i dont know how can export to new tiff images slope/aspect with GDAL ?</p>

<p>the full code :</p>

<pre><code>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 &lt; 0:
                aspect_array[i][j] = 90 - aspect
            elif aspect &gt; 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()
</code></pre>

<p>first trial to export i follow this <a href="https://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html" rel="nofollow">tutorial</a>,but the slope in not correct
have lower degrees from the original(using gis program)</p>

<pre><code>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)
</code></pre>

<p>second trial i follow this <a href="http://stackoverflow.com/questions/37648439/simplest-way-to-save-array-into-raster-file-in-python">quest</a> but i dont take some export</p>

<pre><code>def array_to_raster(slope_array):
    """Array &gt; 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) 
</code></pre>

<p>any idea</p>

<p>The error states exactly what the problem is. You are trying to multiply an <code>int</code> type by a <code>NoneType</code> (None). The most likely case is that <code>nodataval</code> is None, which would occur because the NoDataValue for filename's first raster band is not defined. Your print(nodataval) command should demonstrate that. <strike>Remember, printing None doesn't appear as the string 'None'. It appears as a blank or no character.</strike></p>

<p>Your edits show that nodataval is None.</p>

python python-2.7 numpy scipy