Loading autocnet/cg/cg.py +26 −26 Original line number Diff line number Diff line Loading @@ -8,15 +8,16 @@ import geopandas as gpd import ogr from skimage import transform as tf from scipy.spatial import ConvexHull from scipy.spatial import Voronoi from scipy.spatial import Voronoi, Delaunay, ConvexHull import shapely.geometry from shapely.geometry import Polygon, MultiPolygon, Point from shapely.affinity import scale from shapely import wkt from shapely import wkt, geometry from autocnet.utils import utils from autocnet.cg import cg from shapely.ops import cascaded_union, polygonize def two_point_extrapolate(x, xs, ys): Loading Loading @@ -570,7 +571,7 @@ def alpha_shape(points, alpha): 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)) triangles = list(cg.polygonize(m)) return cascaded_union(triangles) Loading Loading @@ -627,4 +628,3 @@ def rasterize_polygon(shape, vertices, dtype=bool): # Set all values inside polygon to one base_array[fill] = 1 return base_array autocnet/cg/change_detection.py +146 −78 Original line number Diff line number Diff line Loading @@ -11,7 +11,11 @@ from sklearn.cluster import DBSCAN from sklearn.neighbors import NearestNeighbors from skimage.feature import blob_log, blob_doh from math import sqrt, atan2, pi from hoggorm.mat_corr_coeff import RVcoeff import math import scipy from scipy.spatial import cKDTree from plio.io.io_gdal import GeoDataset Loading @@ -19,11 +23,11 @@ from plio.io.io_gdal import GeoDataset from shapely import wkt from shapely.geometry import Point, MultiPoint, Polygon import richdem as rd import pandas as pd import geopandas as gpd from math import sqrt, atan2, pi import pysis from autocnet.utils.utils import bytescale Loading Loading @@ -469,6 +473,70 @@ def blob_detector(image1, image2, sub_solar_azimuth, image_func=image_diff_sq, changes = gpd.GeoDataFrame(geometry=polys) return changes, bdiff def rv_detector(im1, im2, search_size, pattern_size=None, threshold=.999): """ RV coefficient based change detection. Parameters ---------- image1 : GeoDataset Image representing the "before" state of the ROI, can be a 2D numpy array or plio GeoDataset. image2 : GeoDataset Image representing the "after" state of the ROI, can be a 2D numpy array or plio GeoDataset. search_size : int The size of the search space surrounding each pixel. pattern_size : int The size of the window used to calculate the RV score. This window slides through the search space. threshold : float The cutoff value for an RV value to be considered a change : pd.DataFrame A pandas dataframe containing a points of changed areas : np.ndarray A numpy array containing the RV values of each pixel. Note that the array is padded by NaN values for 1/2 window size on each size """ def get_window(arr, ulx, uly, size): return arr[ulx:ulx+size, uly:uly+size] if isinstance(im1, GeoDataset): im1 = im1.read_array() if isinstance(im2, GeoDataset): im2 = im2.read_array() if pattern_size is None: pattern_size = search_size if search_size < pattern_size: print("Pattern size must be <= search size. Setting pattern_size=search_size") search_size = pattern_size rv = np.empty(im1.shape) rv[:] = np.NaN for row in range(im1.shape[0] - search_size): for col in range(im1.shape[1] - search_size): best = -float("inf") # Windows are determined by ulx, uly, but rv corresponds to window's center center_x = col + (search_size//2) center_y = row + (search_size//2) for row_offset in range(search_size - pattern_size + 1): for col_offset in range(search_size - pattern_size + 1): pattern_uly = row + row_offset pattern_ulx = col + col_offset best = max(best, abs(RVcoeff([get_window(im1, row, col, pattern_size), get_window(im2, pattern_uly, pattern_ulx, pattern_size)])[0,1])) rv[center_y, center_x] = best # Get x/y coordinates of points with correlation <= threshold filtered_rv = np.asarray(np.where(rv<=threshold)).T change_geometries = gpd.GeoDataFrame(geometry=[Point(x[1],x[0]) for x in filtered_rv]) return change_geometries, rv def compute_depression(input_dem, scale_factor=1, curvature_percentile=75, return_polygon=True, alpha=0.5): """ Loading Loading @@ -550,7 +618,7 @@ def compute_depression(input_dem, scale_factor=1, curvature_percentile=75, retur dem = dem+demfilled if return_polygon: concave_hull = cg.alpha_shape(np.argwhere(dmask), alpha=alpha) concave_hull = cg.cg.alpha_shape(np.argwhere(dmask), alpha=alpha) return dem, concave_hull return dem, dmask Loading Loading @@ -584,12 +652,12 @@ def generate_dem(alpha=1.0, size=800, scales=[160,80,32,16,8,4,2,1], scale_facto """ topo=np.zeros((2,2))+random.rand(2,2)*(200/(2.**alpha)) topo=np.zeros((2,2))+np.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 = topo + np.random.rand(int(nn), int(nn))*(200/(nn**alpha)) topo = rd.rdarray(topo, no_data=0) Loading Loading @@ -710,6 +778,7 @@ def generate_boulder(dem, radius, height=None, x=None, y=None): ''' Generates a half dome with a given radius, at a given height, at a given x, y in 2D topology array Parameters ---------- dem : 2d array Loading Loading @@ -843,4 +912,3 @@ 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 environment.yml +1 −0 Original line number Diff line number Diff line Loading @@ -11,6 +11,7 @@ dependencies: - dill - geoalchemy2 - geopandas - hoggorm - imageio - ipykernel - jupyter Loading Loading
autocnet/cg/cg.py +26 −26 Original line number Diff line number Diff line Loading @@ -8,15 +8,16 @@ import geopandas as gpd import ogr from skimage import transform as tf from scipy.spatial import ConvexHull from scipy.spatial import Voronoi from scipy.spatial import Voronoi, Delaunay, ConvexHull import shapely.geometry from shapely.geometry import Polygon, MultiPolygon, Point from shapely.affinity import scale from shapely import wkt from shapely import wkt, geometry from autocnet.utils import utils from autocnet.cg import cg from shapely.ops import cascaded_union, polygonize def two_point_extrapolate(x, xs, ys): Loading Loading @@ -570,7 +571,7 @@ def alpha_shape(points, alpha): 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)) triangles = list(cg.polygonize(m)) return cascaded_union(triangles) Loading Loading @@ -627,4 +628,3 @@ def rasterize_polygon(shape, vertices, dtype=bool): # Set all values inside polygon to one base_array[fill] = 1 return base_array
autocnet/cg/change_detection.py +146 −78 Original line number Diff line number Diff line Loading @@ -11,7 +11,11 @@ from sklearn.cluster import DBSCAN from sklearn.neighbors import NearestNeighbors from skimage.feature import blob_log, blob_doh from math import sqrt, atan2, pi from hoggorm.mat_corr_coeff import RVcoeff import math import scipy from scipy.spatial import cKDTree from plio.io.io_gdal import GeoDataset Loading @@ -19,11 +23,11 @@ from plio.io.io_gdal import GeoDataset from shapely import wkt from shapely.geometry import Point, MultiPoint, Polygon import richdem as rd import pandas as pd import geopandas as gpd from math import sqrt, atan2, pi import pysis from autocnet.utils.utils import bytescale Loading Loading @@ -469,6 +473,70 @@ def blob_detector(image1, image2, sub_solar_azimuth, image_func=image_diff_sq, changes = gpd.GeoDataFrame(geometry=polys) return changes, bdiff def rv_detector(im1, im2, search_size, pattern_size=None, threshold=.999): """ RV coefficient based change detection. Parameters ---------- image1 : GeoDataset Image representing the "before" state of the ROI, can be a 2D numpy array or plio GeoDataset. image2 : GeoDataset Image representing the "after" state of the ROI, can be a 2D numpy array or plio GeoDataset. search_size : int The size of the search space surrounding each pixel. pattern_size : int The size of the window used to calculate the RV score. This window slides through the search space. threshold : float The cutoff value for an RV value to be considered a change : pd.DataFrame A pandas dataframe containing a points of changed areas : np.ndarray A numpy array containing the RV values of each pixel. Note that the array is padded by NaN values for 1/2 window size on each size """ def get_window(arr, ulx, uly, size): return arr[ulx:ulx+size, uly:uly+size] if isinstance(im1, GeoDataset): im1 = im1.read_array() if isinstance(im2, GeoDataset): im2 = im2.read_array() if pattern_size is None: pattern_size = search_size if search_size < pattern_size: print("Pattern size must be <= search size. Setting pattern_size=search_size") search_size = pattern_size rv = np.empty(im1.shape) rv[:] = np.NaN for row in range(im1.shape[0] - search_size): for col in range(im1.shape[1] - search_size): best = -float("inf") # Windows are determined by ulx, uly, but rv corresponds to window's center center_x = col + (search_size//2) center_y = row + (search_size//2) for row_offset in range(search_size - pattern_size + 1): for col_offset in range(search_size - pattern_size + 1): pattern_uly = row + row_offset pattern_ulx = col + col_offset best = max(best, abs(RVcoeff([get_window(im1, row, col, pattern_size), get_window(im2, pattern_uly, pattern_ulx, pattern_size)])[0,1])) rv[center_y, center_x] = best # Get x/y coordinates of points with correlation <= threshold filtered_rv = np.asarray(np.where(rv<=threshold)).T change_geometries = gpd.GeoDataFrame(geometry=[Point(x[1],x[0]) for x in filtered_rv]) return change_geometries, rv def compute_depression(input_dem, scale_factor=1, curvature_percentile=75, return_polygon=True, alpha=0.5): """ Loading Loading @@ -550,7 +618,7 @@ def compute_depression(input_dem, scale_factor=1, curvature_percentile=75, retur dem = dem+demfilled if return_polygon: concave_hull = cg.alpha_shape(np.argwhere(dmask), alpha=alpha) concave_hull = cg.cg.alpha_shape(np.argwhere(dmask), alpha=alpha) return dem, concave_hull return dem, dmask Loading Loading @@ -584,12 +652,12 @@ def generate_dem(alpha=1.0, size=800, scales=[160,80,32,16,8,4,2,1], scale_facto """ topo=np.zeros((2,2))+random.rand(2,2)*(200/(2.**alpha)) topo=np.zeros((2,2))+np.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 = topo + np.random.rand(int(nn), int(nn))*(200/(nn**alpha)) topo = rd.rdarray(topo, no_data=0) Loading Loading @@ -710,6 +778,7 @@ def generate_boulder(dem, radius, height=None, x=None, y=None): ''' Generates a half dome with a given radius, at a given height, at a given x, y in 2D topology array Parameters ---------- dem : 2d array Loading Loading @@ -843,4 +912,3 @@ 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
environment.yml +1 −0 Original line number Diff line number Diff line Loading @@ -11,6 +11,7 @@ dependencies: - dill - geoalchemy2 - geopandas - hoggorm - imageio - ipykernel - jupyter Loading