Unverified Commit c7e9fec1 authored by Kelvin Rodriguez's avatar Kelvin Rodriguez Committed by GitHub
Browse files

added dem generation code (#521) 

* added topo generation code

* comments

* added option to return changes as polygon
parent 186c6595

autocnet/cg/cg.py

+102 −0
Original line number Diff line number Diff line
@@ -526,3 +526,105 @@ def distribute_points_in_geom(geom, method="classic", 
    else:

        print('WTF Willy')

    return valid





def alpha_shape(points, alpha):

    """

    Compute a convex hull from a set of points. 



    credit: https://gist.github.com/dwyerk/10561690



    Parameters

    ----------

    

    points : np.array

             points in the format [[x0,y0], [x1, y,1], ... [xn, yn]]

    

    alpha : float 

            Higher alphas creater a tighter the boundery around input points, alphas approaching 0 create a convex hull. 

            Best to keep in the range (0, 1]





    Returns 

    -------



    convex_hull : shapely.geometry.Polygon 

                  Shapely polygon of the convex hull 



    """

    tri = Delaunay(points)

    triangles = points[tri.vertices]

    a = ((triangles[:,0,0] - triangles[:,1,0]) ** 2 + (triangles[:,0,1] - triangles[:,1,1]) ** 2) ** 0.5

    b = ((triangles[:,1,0] - triangles[:,2,0]) ** 2 + (triangles[:,1,1] - triangles[:,2,1]) ** 2) ** 0.5

    c = ((triangles[:,2,0] - triangles[:,0,0]) ** 2 + (triangles[:,2,1] - triangles[:,0,1]) ** 2) ** 0.5

    s = ( a + b + c ) / 2.0

    areas = (s*(s-a)*(s-b)*(s-c)) ** 0.5

    circums = a * b * c / (4.0 * areas)

    

    # avoid devide by zero

    thresh = 1.0/alpha if alpha is not 0 else circums.max()

    filtered = triangles[circums < thresh]

    

    edge1 = filtered[:,(0,1)]

    edge2 = filtered[:,(1,2)]

    edge3 = filtered[:,(2,0)]

    edge_points = np.unique(np.concatenate((edge1,edge2,edge3)), axis = 0).tolist()

    m = geometry.MultiLineString(edge_points)

    triangles = list(polygonize(m))

    return cascaded_union(triangles)





def rasterize_polygon(shape, vertices, dtype=bool):

    """

    Simple tool to convert poly into a boolean numpy array.

    

    source: https://stackoverflow.com/questions/37117878/generating-a-filled-polygon-inside-a-numpy-array

    

    Parameters

    ----------

    

    shape : tuple 

            size of the array in (y,x) format

    

    vertices : np.array, list

               array of vertices in [[x0, y0], [x1, y1]...] format

    

    dtype : type

            datatype of output mask

    

    Returns

    -------

    

    mask : np.array

           mask with filled polygon set to true

    

    """

    def check(p1, p2, base_array):

        idxs = np.indices(base_array.shape) # Create 3D array of indices



        p1 = p1.astype(float)

        p2 = p2.astype(float)



        # Calculate max column idx for each row idx based on interpolated line between two points

        if p1[0] == p2[0]:

            max_col_idx = (idxs[0] - p1[0]) * idxs.shape[1]

            sign = np.sign(p2[1] - p1[1])

        else:

            max_col_idx = (idxs[0] - p1[0]) / (p2[0] - p1[0]) * (p2[1] - p1[1]) + p1[1]

            sign = np.sign(p2[0] - p1[0])

            

        return idxs[1] * sign <= max_col_idx * sign



    base_array = np.zeros(shape, dtype=dtype)  # Initialize your array of zeros



    fill = np.ones(base_array.shape) * True  # Initialize boolean array defining shape fill



    # Create check array for each edge segment, combine into fill array

    for k in range(vertices.shape[0]):

        fill = np.all([fill, check(vertices[k-1], vertices[k], base_array)], axis=0)

    

    print(fill.any())

    # Set all values inside polygon to one

    base_array[fill] = 1

    return base_array

autocnet/cg/change_detection.py

+253 −6
Original line number Diff line number Diff line 
import numpy as np



import matplotlib

from matplotlib.path import Path

from matplotlib import pyplot as plt



import cv2



from sklearn.cluster import  OPTICS

from plio.io.io_gdal import GeoDataset

from scipy.spatial import cKDTree

from sklearn.cluster import DBSCAN

from sklearn.neighbors import NearestNeighbors



from skimage.feature import blob_log, blob_doh

from math import sqrt, atan2, pi



import matplotlib

from matplotlib import pyplot as plt

from scipy.spatial import cKDTree



from plio.io.io_gdal import GeoDataset



from shapely import wkt

from shapely.geometry import Point, MultiPoint, Polygon



import pandas as pd

import geopandas as gpd



from math import sqrt, atan2, pi



import pysis



from autocnet.utils.utils import bytescale

from autocnet.matcher.cpu_extractor import extract_features



from autocnet import cg



def image_diff(arr1, arr2):

     arr1 = arr1.astype("float32")
@@ -461,6 +470,243 @@ def blob_detector(image1, image2, sub_solar_azimuth, image_func=image_diff_sq, 
     return changes, bdiff





def compute_depression(input_dem, scale_factor=1, curvature_percentile=75, return_polygon=True, alpha=0.5):

    """

    Compute depressions and return a new image with largest depressions filled in. 

    

    Parameters

    ----------

    

    input_dem : np.array, rd.rdarray

                2d array of elevation DNs, a DEM

    

    scale_factor : float

                   Value to scale the erotion of planform curvatures by

                   

    curvature_percentile : float 

                           what percentile of the curvature to keep, lower values

                           results in bigger blobs 

                   

    

    Returns

    -------

    dem : rd.rdarray

          Dem with filled depressions

    

    mask : np.array

           Change mask, true on pixels that have been changed 

    

    

    """

    if isinstance(input_dem, np.ndarray):

        dem = rd.rdarray(input_dem.copy(), no_data=0)

    elif isinstance(input_dem, rd.rdarray):

        # take ownership of the reference

        dem = input_dem.copy()



    # create filled DEM

    demfilled = rd.FillDepressions(dem, epsilon=True, in_place=False, topology="D8")

    

    # Mask out filled areas

    mask = np.abs(dem-demfilled)

    thresh = np.percentile(mask, 95)

    mask[mask <= thresh] = False

    mask[mask > thresh] = True

    

    curvatures = rd.TerrainAttribute(dem, attrib='planform_curvature')

    curvatures = (curvatures - np.min(curvatures))/np.ptp(curvatures) 

    curvatures[curvatures < np.percentile(curvatures, curvature_percentile)] = 0

    curvatures[mask.astype(bool)] = 0

    

    demfilled -= curvatures * scale_factor

    

    mask = (curvatures+mask).astype(bool)

    

    # Get 3rd nn distance 

    coords = np.argwhere(mask)

    nbrs = NearestNeighbors(n_neighbors=3, algorithm='kd_tree').fit(coords)

    dists, _ = nbrs.kneighbors(coords)

    eps = np.percentile(dists, 95)

    

    # Cluster

    db = DBSCAN(eps=eps, min_samples=3).fit(coords)

    labels = db.labels_

    unique, counts = np.unique(labels, return_counts=True)

    

    # First count are outliers, ignore

    counts = counts[1:]

    unique = unique[1:]

    

    index = np.argwhere(counts == counts.max())

    group = unique[index][0][0]

    cluster = coords[labels == group]

    

    # mask out depression

    dmask = np.full(dem.shape, False)

    dmask[[*cluster.T]] = True

    

    dem[dmask] = 0

    demfilled[~dmask] = 0

    dem = dem+demfilled



    if return_polygon: 

        concave_hull = cg.alpha_shape(np.argwhere(dmask), alpha=alpha)

        return dem, concave_hull 



    return dem, dmask





def generate_dem(alpha=1.0, size=800, scales=[160,80,32,16,8,4,2,1], scale_factor=5):

    """

    Produces a random DEM

    

    Parameters

    ----------

    

    alpha : float 

            Controls height variation. Lower number makes a shallower and noisier DEM, 

            higher values create smoother DEM with large peaks and valleys. 

            Reccommended range = (0, 1.5]

    

    size : int

           size of DEM, output DEM is in the shape of (size, size)

    

    scale_factor : float 

                   scalar to multiply the slope degridation by, higher values = more erotion. 

                   Reccomended to increase proportionately with alpha 

                   (higher alphas mean you might want higher scale_factor)

    

    Returns 

    -------

    

    dem : np.array 

          DEM array in the shape (size, size)

    

    """

    

    topo=np.zeros((2,2))+random.rand(2,2)*(200/(2.**alpha))



    for k in range(len(scales)):

        nn = size/scales[k]

        topo = scipy.misc.imresize(topo, (int(nn), int(nn)), "cubic", mode="F")

        topo = topo + random.rand(int(nn), int(nn))*(200/(nn**alpha))

    

    topo = rd.rdarray(topo, no_data=0)

    

    curvatures = rd.TerrainAttribute(topo, attrib='slope_riserun')

    curvatures = (curvatures - np.min(curvatures))/np.ptp(curvatures) * scale_factor

    return topo - curvatures





def hillshade(img, azi=255, alt=60, min_slope=20, max_slope=100, min_bright=0, grayscale=False):

    """

    hillshade a DEM, based on IDL code by Colin Dundas translated by Adam Paquette 

    

    Parameters

    ----------

    

    img : np.array

          DEM to hillshade

    

    azi : float 

          Sun azimuth in degrees 

    

    alt: float 

         base alt

    

    min_slope : float 

                minimum slope value 

    

    max_slope : float 

                maximum slope value 

    

    min_bright : float 

                 minimum brightness 

    

    grayscale : bool 

                whether or not to produce grayscale image 

    

    

    Returns

    -------

    

    dem : np.array 

          hillshaded DEM 

    

    """

    dem = np.array(np.flip(bytescale(img), axis = 0), dtype=int)

    emax = np.max(dem)

    emin = np.min(dem)



    indices = np.linspace(0, 255, 256) / 25.5



    red_array = [0,25,50,101,153,204,255,255,255,255,255,255]

    red_index = np.arange(len(red_array))

    red_vec = np.interp(indices, red_index, red_array)



    green_array = [42,101,153,204,237,255,255,238,204,153,102,42]

    green_index = np.arange(len(green_array))

    green_vec = np.interp(indices, green_index, green_array)



    blue_array = [255,255,255,255,255,255,204,153,101,50,25,0]

    blue_index = np.arange(len(blue_array))

    blue_vec = np.interp(indices, blue_index, blue_array)



    zz = (255.0/(emax-emin))*(dem-emin)

    zz = zz.astype(int)



    nx = (np.roll(dem, 1, axis = 1) - dem)

    ny = (np.roll(dem, 1, axis = 0) - dem)

    sz = np.shape(nx)

    nz = np.ones(sz)

    nl = np.sqrt(np.power(nx, 2.0) + np.power(ny, 2.0) + np.power(nz, 2.0))

    nx = nx/nl

    ny = ny/nl

    nz = nz/nl



    azi_rad = math.radians(azi)

    alt_rad = math.radians(alt)

    lx = math.sin(azi_rad)*math.cos(alt_rad)

    ly = math.cos(azi_rad)*math.cos(alt_rad)

    lz = math.sin(alt_rad)



    dprod = nx*lx + ny*ly + nz*lz



    if min_slope is not None:

        min_dprod = math.cos(math.radians(max_slope + 90.0 - alt))

    else:

        min_dprod = np.min(dprod)



    if max_slope is not None:

        max_dprod = math.cos(math.radians(90.0 - alt - max_slope))

    else:

        max_dprod = np.max(dprod)



    bright = ((dprod - min_dprod) + min_bright)/((max_dprod - min_dprod) + min_bright)



    if grayscale:

        qq=(255*bright)

    else:

        qq = red_vec[zz]*bright



    if grayscale:

        rr = (255*bright)

    else:

        rr = green_vec[zz]*bright



    if grayscale:

        ss=(255*bright)

    else:

        ss = blue_vec[zz]*bright



    arrforout = np.dstack((qq, rr ,ss))

    arrforout = np.flip(arrforout.astype(int), axis = 0)

    arrfotout = bytescale(arrforout)

    arrforout.shape

    return arrforout







def generate_boulder(dem, radius, height=None, x=None, y=None):

    '''

    Generates a half dome with a given radius, at a given height,
@@ -609,3 +855,4 @@ def generate_boulder_field(dem, num_boulders, x_shift_min = 5, x_shift_max = 10, 
            after_polys.append(after_geom)



    return before_dem, before_polys, after_dem, after_polys

autocnet/utils/utils.py

+10 −1
Original line number Diff line number Diff line
@@ -10,6 +10,15 @@ import networkx as nx 


from osgeo import ogr



from sklearn.cluster import DBSCAN

from sklearn.neighbors import NearestNeighbors



from scipy.spatial import Delaunay



from shapely import geometry

from shapely.geometry import MultiPoint

from shapely.ops import cascaded_union, polygonize



def tile(array_size, tilesize=1000, overlap=500):

    stepsize = tilesize - overlap

    if stepsize < 0: