qunt_mesh.py

gistfile1.txt

import matplotlib.cm as cm
import matplotlib.pyplot as plt
import matplotlib.tri as tri
import numpy as np
import numpy.ma as ma
from osgeo import gdal
from scipy import ndimage
from scipy.ndimage.morphology import binary_dilation
from scipy.signal import fftconvolve



class ProcessRaster(object):
    def __init__(self):
        self.NODATA = -32768
        self.raster_array_in
        self.sobel_mass_array








def readraster(filename):
    raster = gdal.Open(filename)
    return raster


def roughness(d, ci):
    return np.abs((d - d[ci])).max()

def sobel(arr):
    return ndimage.filters.sobel(arr)


def reject_outliers(data, m=2):
    return data[abs(data - np.mean(data)) < m * np.std(data)]

def gaussian_blur1d(in_array, size):
        #check validity
        try:
            if 0 in in_array.shape:
                raise Exception("Null array can't be processed!")
        except TypeError:
            raise Exception("Null array can't be processed!")
        # expand in_array to fit edge of kernel
        padded_array = np.pad(in_array, size, 'symmetric').astype(float)
        # build kernel
        x, y = np.mgrid[-size:size + 1, -size:size + 1]
        g = np.exp(-(x**2 / float(size) + y**2 / float(size)))
        g = (g / g.sum()).astype(float)
        # do the Gaussian blur
        out_array = fftconvolve(padded_array, g, mode='valid')
        return out_array.astype(in_array.dtype)

#0= no data 100 = all data
def applyMask(dem, roughness_arr, level, bad_data_mask):
    roughness_arr_masked = ma.array(roughness_arr, mask=bad_data_mask)
    mean_value = roughness_arr.mean()
    #level = level-50
    index_to_mask = level#(float(level)/100)*mean_value*10

    output_mask = ma.masked_less(roughness_arr_masked, index_to_mask).mask
    masked_dem = ma.array(dem, mask=output_mask)

    plotTriang(masked_dem)
    plot(masked_dem)




def plot(arrOut):

    fig = plt.figure()

    ax = fig.add_subplot(111)
    ax.imshow(arrOut, cmap=cm.jet, interpolation='nearest')

    numrows, numcols = arrOut.shape
    def format_coord(x, y):
        col = int(x+0.5)
        row = int(y+0.5)
        if col>=0 and col<numcols and row>=0 and row<numrows:
            z = arrOut[row,col]
            return 'x=%1.4f, y=%1.4f, z=%1.4f'%(x, y, z)
        else:
            return 'x=%1.4f, y=%1.4f'%(x, y)

    ax.format_coord = format_coord
    plt.show()




def plotTriang(arr):
    x, y = np.nonzero(arr.T)
    z = arr.T.ravel()[np.flatnonzero(arr.T)]
    triang = tri.Triangulation(x, y)
    plt.triplot(triang, lw=0.5, color='black')

    plt.title("High-resolution tricontouring")

    plt.show()

    # for t in tri:
    # # t[0], t[1], t[2] are the points indexes of the triangle
    #     t_i = [t[0], t[1], t[2], t[0]]
    #     pylab.plot(x[t_i],y[t_i])
    #
    # pylab.show()






def main():



    inputfile = "/home/dev/SRTM/N43W123.hgt"
    outputfile = "/home/dev/out.ply"

    raster = readraster(inputfile)
    in_arr = raster.GetRasterBand(1).ReadAsArray()
    #in_arr = in_arr[1:1000, 1:1000]
    # get the data as a numpt array ensuring it is floating point data
    in_arr = in_arr.astype(float)

    in_arr = ma.masked_values(in_arr, -32768)
    mask = in_arr.mask
    mask = binary_dilation(mask, iterations=10)
    #
    # #arrOut = arr.filled(arr.mean())
    #
    smoothed = gaussian_blur1d(in_arr, 4)
    # # calculate TRI using a 9x9 moving window
    roughness_arr = sobel(smoothed)#generic_filter(smoothed, roughness, 18, mode='nearest', extra_arguments=(40,))

    arrOut = applyMask(in_arr, roughness_arr, 60, mask)



main()