python (scipy) TypeError

[chrischris201444 at gmail.com]

hello


i try to follow some tutorial but i have that error :

Traceback (most recent call last):
  File "C:\Python27\test\test\earth.py", line 42, in <module>
    slope_array = np.ones_like(data_array) * nodataval
TypeError: unsupported operand type(s) for *: 'int' and 'NoneType'


first i define that lines of code :

data_array = raster2array(filename)
nodataval = getNoDataValue(filename)

and the line with error :

slope_array = np.ones_like(data_array) * nodataval

how can i fix this error ?i change type field ?

on the four prints:

print(resolution)
print(nodataval)
print(type(data_array))
print(data_array.shape)

i take that exports :

{'east-west': 0.0002777777777777778, 'north-south': 0.0002777777777777778}
None
<type 'numpy.ndarray'>
(3601, 3601)

the full 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 < 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()