Source code for raidnr.bwgrid.grid_transformations





# -----------------------------------------------------------------


# ██████╗  █████╗ ██╗██████╗     ███╗   ██╗██████╗ 
# ██╔══██╗██╔══██╗██║██╔══██╗    ████╗  ██║██╔══██╗
# ██████╔╝███████║██║██║  ██║    ██╔██╗ ██║██████╔╝
# ██╔══██╗██╔══██║██║██║  ██║    ██║╚██╗██║██╔══██╗
# ██║  ██║██║  ██║██║██████╔╝    ██║ ╚████║██║  ██║
# ╚═╝  ╚═╝╚═╝  ╚═╝╚═╝╚═════╝     ╚═╝  ╚═══╝╚═╝  ╚═╝
                                               

# -----------------------------------------------------------------


#-----------------------------------------------------

# Grid Transformations v.01
# 21.06.18
# Ioannis Toumpalidis

#-----------------------------------------------------
# Description
#-----------------------------------------------------

# Set of scripts to create and manipulate a geo lattice grid 
# and burn values to it. 
# Grid alligns with the what3word. (Not tested for the alligment, errors might vary due to projection).
# 
#-----------------------------------------------------


import numpy as np
import pandas as pd

import matplotlib.pyplot as plt

import geopandas as gpd
import os


import math

import requests
import json

from shapely.geometry import Point
from pyproj import Proj, transform
from shapely.geometry.linestring import LineString


    
import ogr
import gdal
from PIL import Image

from raidnr.config import TOPO_LAYERS
from multiprocessing import Pool


import gdal2tiles
import os


import osr
from osgeo import gdal_array
import gdal
import json

from scipy import interpolate

from numpy import trapz
from scipy.integrate import simps





def introLogo(name,date,author=""):

    print("______________________________________________________________________________________________")
    print("______________________________________________________________________________________________")
    
    print(
    """ 
    ::::::::  ::::::::  :::::::: 
    ::::::::  ::::::::  ::::::::
    ::::::::::::::::::::::::::::
    ::::::::::::::::::::::::::::
        ::::::::::::::::::::
    ::::::::::       :::::::::::
    ::::::::::       :::::::::::
    ::::::::::       :::::::::::  ######                                       #     # 
    ::::::::::       :::::::::::  #     # #####  #   # #####  ###### #    #    #  #  #  ####   ####  ##### 
        ::::::::::::::::::::      #     # #    #  # #  #    # #      ##   #    #  #  # #    # #    # #    #
    ::::::::::::::::::::::::::::  ######  #    #   #   #    # #####  # #  #    #  #  # #    # #    # #    #
    ::::::::::::::::::::::::::::  #     # #####    #   #    # #      #  # #    #  #  # #    # #    # #    #
    ::::::::  ::::::::  ::::::::  #     # #   #    #   #    # #      #   ##    #  #  # #    # #    # #    #
    ::::::::  ::::::::  ::::::::  ######  #    #   #   #####  ###### #    #     ## ##   ####   ####  #####


    """
    )


    print("______________________________________________________________________________________________")
    print("______________________________________________________________________________________________")
    

    print(""" 

            :----------------------------------------------------------------------:
            :                                                                      :
            :  Title of Project : %s"""%(name))

    if(author!=""):
        print("""
            :  Authored by : %s"""%(author))

    print("""
            :                                                       Date: %s
            :----------------------------------------------------------------------:
    """%(date))



    print("-----------------------------------------------------------------------------------------------")


[docs]def w3wGetGrid(long,lat): """ Return the what3words grid id Parameters ---------- long : float coordinates in WGS84. lat : float coordinates in WGS84. link to :func:`~grid_transformations.w3wGetCoordinates` """ url="https://api.what3words.com/v3/convert-to-3wa?coordinates=%s,%s&key=78O2F1VX"%(long,lat) response=requests.get(url) return json.loads(response.text)
[docs]def w3wGetCoordinates(w3w): """ Parameters ---------- w3w: string w3w in the format of "word1.word2.word3" List of strings containing the what 3 words for a specific area. """ url= ("https://api.what3words.com/v3/convert-to-coordinates?"+ "words=%s&key=78O2F1VX"%(w3w)) response=requests.get(url) coords=json.loads(response.text)['coordinates'] wgsCoords=(coords['lng'],coords['lat']) return wgsCoords
[docs]def wgs_to_bng(center): """ Project to BNG Parameters ---------- center: tupple format longitude, latitude """ inProj = Proj(init='epsg:4326') outProj = Proj(init='epsg:27700') x1,y1 = center x2,y2 = transform(inProj,outProj,x1,y1) return(x2,y2)
[docs]def bng_to_wgs(center): """ Transfrom coordinates from BNG to WGS Parameters ---------- center: tupple format longitude, latitude """ inProj = Proj(init='epsg:27700') outProj = Proj(init='epsg:4326') x1,y1 = center x2,y2 = transform(inProj,outProj,x1,y1) return(x2,y2)
def midpoint(p1, p2): return Point((p1[0]+p2[0])/2, (p1[1]+p2[1])/2)
[docs]def try_bounds(polygon): """ Extract bounding box for shapely polygon. Parameters ---------- polygon : shapely.geometry.Polygon """ try: return polygon.exterior.bounds except: return 0
[docs]def get_Bounds_FromShape(geopandasDf, copyDf=False): """ Get Bounding Box from a geopandas dataframe Parameters ---------- bui """ if(type(geopandasDf['geometry'].iloc[0])==LineString or type(geopandasDf['geometry'].iloc[0])==Point): geopandasDf_bounds=geopandasDf.bounds geopandasDf=geopandasDf.join(geopandasDf_bounds,how='outer') minx=geopandasDf['minx'] maxx=geopandasDf['maxx'] miny=geopandasDf['miny'] maxy=geopandasDf['maxy'] else: geopandasDf['bbox']=geopandasDf['geometry'].apply(lambda x: try_bounds(x)) geopandasDf=geopandasDf[geopandasDf['bbox']!=0] minx,miny,maxx,maxy=zip(*geopandasDf['bbox']) bound_minx=min(minx) bound_miny=min(miny) bound_maxx=max(maxx) bound_maxy=max(maxy) if copyDf: return geopandasDf,[bound_minx,bound_miny,bound_maxx,bound_maxy] else: return [bound_minx,bound_miny,bound_maxx,bound_maxy]
[docs]def get_UG_grid_bounds(geopandasDf,corner_stone=[531087.2503762188, 181782.11191701167],offset=3,corners=False): """ Return the bounding box of a geoDataFrame in coordinates, alligned with the Universal Grid. Parameters ---------- geopadnasDf : (geopandas.DataFrame) corner_stone (list) : the coordinates of the zeroing poing is the SE corner of the Zero cell. offset (int) : dimensions of the cell. corners (bool) : If True return the corners of the boundings box, else the centroid of the corner cells. """ bounds_minx,bounds_miny,bounds_maxx,bounds_maxy=get_Bounds_FromShape(geopandasDf) bounds_NW_corner=(bounds_minx,bounds_maxy) bounds_SE_corner=(bounds_maxx,bounds_miny) # grid_NW_index= get_UG_index(bounds_NW_corner,offset,corner_stone) # grid_SE_index= get_UG_index(bounds_SE_corner,offset,corner_stone) cpNW = get_UG_centroid_from_coords(bounds_NW_corner,corner_stone,offset) cpSE =get_UG_centroid_from_coords(bounds_SE_corner,corner_stone,offset) if (corners==True): print("returning corners for NW and SE corner cells ") NWcorner= returnSquareCoords(cpNW,offset) SECorner= returnSquareCoords(cpSE,offset) return(NWcorner,SECorner) else: print("returning centroid coordinates for NW and SE corner cells ") return (cpNW,cpSE)
def get_UG_index(NW_corner,dx,corner_stone=[531087.2503762188, 181782.11191701167]): ii=math.floor((NW_corner[0]-corner_stone[0])/dx) jj=math.floor((NW_corner[1]-corner_stone[1])/dx) return ii,jj def get_UG_centroid_from_index(i,j,corner_stone,dx): NWx = corner_stone[0]+(i*dx)+(dx*0.5) NWy = corner_stone[1]+(j*dx)+(dx*0.5) return NWx,NWy def get_UG_centroid_from_coords(x,c=[531087.2503762188, 181782.11191701167],offset=3): i,j=get_UG_index(x,offset,c) NWx,NWy= get_UG_centroid_from_index(i,j,c,offset) return NWx,NWy
[docs]def get_UG_centroid_from_wgs(coords): """ wrapper for the get_UG_centroid_from_coords function. Parameters ---------- coords : array-like Coords in WGS (lat, long) Return ------ coords_ : tupple Coord tupple in WGS (lat,long) """ print("NOTICE: coords in lat long format") lon = coords[1] lat = coords[0] coords_bng= wgs_to_bng((lon,lat)) resp=get_UG_centroid_from_coords(coords_bng) coords_temp=bng_to_wgs(resp) coords_ = (coords_temp[1],coords_temp[0]) return coords_
[docs]def draw_UG_from_points(coordsX,coordsY,offset,ax,color): """ Draw UG given two lists for x,y coordinates of points. Parameters ---------- coordsX: list coordsY: list offset : the width of the grid cells ax: ax of the plot color : color for the lines of the grid """ for i in range(len(coordsX)): ccp= gt.get_UG_centroid_from_coords((coordsX[i],coordsY[i]),corner_stone,offset) ax.scatter(ccp[0],ccp[1],c='red',s=0.4) small_sqX,small_sqY=returnSquareCoords(ccp,offset) ax.plot(small_sqX,small_sqY,color=color)
[docs]def return_grid_cp(startX,endX,startY,endY,offset): """ Generate and return the centroids of a grid given two sets of coordinates, that specify the bounding box of the grid. Parameters ---------- startX,endX : float horizontal coordinates startY,endY : float vertical coordinates offset : int width of the cells """ gridXX=np.arange(startX,endX,offset) gridYY= np.arange(startY,endY,-offset) g=np.meshgrid(gridXX,gridYY) cells= list(zip(*(x.flat for x in g))) x,y=zip(*cells) return x,y
[docs]def returnSquareCoords(centroid,width): """ Return list of coordinates (X,Y) for the corners of a square with center the given centroid. Parameters ---------- centroid: array [x,y] width: float width of the cell """ squareX=[centroid[0]-width*0.5,centroid[0]+width*0.5, centroid[0]+width*0.5,centroid[0]-width*0.5, centroid[0]-width*0.5] squareY=[centroid[1]+width*0.5,centroid[1]+width*0.5, centroid[1]-width*0.5,centroid[1]-width*0.5, centroid[1]+width*0.5] return squareX,squareY
[docs]def rasterize_shape(filename, bbox, tifName="", attribute=""): """ Burn shape to a raster defined by the input bbox. Parameters ---------- filename: str path to the .shp file bbox: list bounding box for the raster [NWCorner_x, SECorner_y, SECorner_x, NWCorner_y] tifName : str [optional] """ vector_fn = filename # Define pixel_size and NoData value of new raster pixel_size = 3 NoData_value = -9999 # Open the data source and read in the extent source_ds = ogr.Open(vector_fn) source_layer = source_ds.GetLayer() x_min,x_max,y_min,y_max= bbox[0], bbox[2], bbox[1], bbox[3] # Create the destination data source x_res = int((x_max - x_min) / pixel_size) y_res = int((y_max - y_min) / pixel_size) # print (x_res, y_res) val=255 OPTIONS='ALL_TOUCHED=TRUE' if(attribute!=""): OPTIONS+=",ATTRIBUTE="+attribute memoryPref=gdal.GDT_Byte val=0 else: memoryPref=gdal.GDT_Byte # if tifName!="": target_ds = gdal.GetDriverByName("GTiff").Create(tifName, x_res, y_res,memoryPref) else: target_ds = gdal.GetDriverByName(str('MEM')).Create('', x_res, y_res, memoryPref) # target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0,-pixel_size)) band = target_ds.GetRasterBand(1) band.SetNoDataValue(NoData_value) # Rasterize if(attribute!=""): gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[val], options=[OPTIONS]) else: gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[val],options=[OPTIONS]) # print("normal") # Read as array array = band.ReadAsArray() target_ds=None return array
[docs]def return_cost(topo_array, layerId, wx, wy, ksize=21): """ Return the sum of the cost at wx, wy positions, given a specifc kernel-size. Parameters ---------- topo_array: nd-array terrain array layerID: int index to access the desired layer wx : int or array-like raster x-coordinates on the terrain wy : int or array-like raster y-coordinates on the terrain ksize : int kernel size to sum the values. Must be an odd number Returns ------- array array of sum """ ksize = (ksize-1)//2 #take the distance focalPt-to-kernel border sum_cost=[] if(isinstance(wx,int) and isinstance(wy,int)): wx=np.array([wx]) wy=np.array([wy]) #iterate through focuses if(len(topo_array.shape)<3): to_topo_array = np.ndarray(shape=(topo_array.shape[0],topo_array.shape[1],1)) to_topo_array[:,:,0] = topo_array topo_array = to_topo_array.copy() layer_array = np.pad(topo_array[:,:,layerId],ksize,mode='constant') for focal_pt in range(len(wx)): try: # since we applied the pad ROI=layer_array[ int(wy[focal_pt]) : int(wy[focal_pt])+2*ksize+1, int(wx[focal_pt]) : int(wx[focal_pt])+2*ksize+1] sum_cost.append(sum(sum(ROI))) except Exception as e: print(e) sum_cost.append(1) # add padding # return np.array(sum_cost) return np.array(sum_cost)
[docs]def get_grid_from_directory(raster_path, raster_size=None): """ Construct cost surface from raster grids. Parameters ---------- raster_path: str folder path to rasters raster_size: str Path to the directory where rasters are. Returns ------- nd-array: terrain layered array containing the topography """ if not raster_size: raster_size = get_rasters_size(raster_path) multi_d_array = np.empty(shape=(raster_size[0],raster_size[1],len(TOPO_LAYERS))) for i,tl in enumerate(TOPO_LAYERS): try: rasterFile= os.path.join(raster_path,"Processed_"+tl+'.tif') ds = gdal.Open(rasterFile) myarray = np.array(ds.GetRasterBand(1).ReadAsArray()) multi_d_array[:,:,i]=myarray except: multi_d_array[:,:,i]=np.zeros(shape=raster_size) return multi_d_array
[docs]def get_rasters_size(path): """ """ assert os.path.exists(path), "Path {} does not exist".format(path) file_paths = [os.path.join(path,f) for f in os.listdir(path) if os.path.isfile(os.path.join(path,f)) and f.endswith(('.tif', '.tiff', '.jpg', '.png'))] #take all files shapes = [Image.open(fp).size for fp in file_paths] shapes = np.array(shapes).transpose() assert(len(np.unique(shapes[0])) == 1 and len(np.unique(shapes[1])) == 1), "Rasters in path come in different sizes" r = shapes.transpose()[0] return (np.flip(r)) #rasters have col,row
[docs]def coords_geo_to_raster(bbox, coords, pixels_size=3., to_int=True): """ Transfrom coordinates to raster row,col coords. Parameters ---------- bbox : array minx, miny, maxx, maxy coords : tuple In format (long, lat) """ if((coords[1]<=bbox[3]) and (coords[1]>=bbox[1]) and (coords[0]>=bbox[0]) and (coords[0]<=bbox[2])): raster_coord_j=(bbox[3]-coords[1])/pixels_size raster_coord_i=(coords[0]-bbox[0])/pixels_size if True: point= (int(raster_coord_j),int(raster_coord_i)) else: point= (raster_coord_j,raster_coord_i) return point else: raise ValueError("coord not in bbox")
[docs]def raster_to_coords(bbox,raster_coords,pixel_size=3): """ Transfrom coordinates from raster to coords. Parameters ---------- bbox : array minx,maxx,miny,maxy raster_coords : tuple In format (row,col) Returns ------- tbc """ latitude=bbox[3]-raster_coords[0]*pixel_size longitude= bbox[0]+raster_coords[1]*pixel_size point = (longitude,latitude) return point
[docs]def w3w_to_raster(w3w_list,bbox,fileName="",cost_surface="",w=20): """ Return the coresponding raster coordinates from a w3w location. Parameters ---------- w3w_list : list list of words eg. w3w_list=['apple','grid','row'] bbox : list list containing [minx,maxx,miny,maxy] """ coords=w3wGetCoordinates(w3w_list)['coordinates'] wgsCoords=(coords['lng'],coords['lat']) # print(wgsCoords) bngCoords=wgs_to_bng(wgsCoords) raster_coords=coords_geo_to_raster(bbox,bngCoords) print(raster_coords) if(fileName!=""): if(cost_surface!=""): return_plot_grid(cost_surface,raster_coords,w) else: print("Error: provide cost_surface numpy array") else: return raster_coords
[docs]def return_plot_grid(sumarray,raster_coords,fileName="",w=20): """ """ fig,ax = plt.subplots(figsize=(10,10)) ax.axis("equal") j,i= raster_coords if(w>100): w=100 print("maximum width of plot is 200 grids") subset=sumarray[j-w:j+w+1,i-w:i+w+1] gridX= np.arange(0.5,subset.shape[1]+0.5) gridY= np.arange(0.5,subset.shape[0]+0.5) ax.set_xticks(gridX) ax.set_yticks(gridY) ax.scatter(w,w,c='red',marker='s') ax.axes.xaxis.set_ticklabels([]) ax.axes.yaxis.set_ticklabels([]) ax.imshow(subset) if(fileName!=""): plt.grid() plt.savefig(file_name+".png") else: plt.grid() plt.show()
[docs]def distance_transform(world,movements,size=(3,3),backward=False,debug=False): """ Implementation of the distance transform algorithm. Parameters --------- world: np.array raster array movements : dict dictionary to save the movement of the cells backward : boolean perform backward scan, default is forward debug : boolean debug the process with plots for every itteration References --------- """ centerx = int(size[0]/2) distances = np.zeros(shape=size) for i in range(len(distances)): for j in range(len(distances[i])): distances[i,j]=math.sqrt(math.pow(int(size[0]/2)-i,2) + math.pow(int(size[1]/2)-j,2) ) if(backward): origin_col,destination_col,step_col =[(len(world)-centerx-1),centerx-1,-1] origin_row,destination_row,step_row =[len(world[i])-centerx-1,centerx-1,-1] fromThis,toThat = 1,0 distances = [[np.inf,np.inf,np.inf], [np.inf,np.inf,1], [math.sqrt(2),1,math.sqrt(2)]] fromThis,toThat = 1,1 else: origin_col,destination_col,step_col =[centerx,len(world)-centerx,1] origin_row,destination_row,step_row =[centerx,len(world[i])-centerx,1] fromThis,toThat = 0,1 distances = [ [math.sqrt(2), 1 , math.sqrt(2)], [1,np.inf,np.inf], [np.inf,np.inf,np.inf] ] fromThis,toThat = 1,1 # print(distances) # distances[centerx][centerx]+=np.inf for i in range(origin_col,destination_col,step_col): for j in range(origin_row,destination_row,step_row): if(world[i][j]!=0): subtest=world[i-centerx*fromThis:i+centerx*toThat+1,j-centerx*fromThis:j+centerx*toThat+1]+distances flat_indx=np.argmin(subtest) minInd_= np.unravel_index(flat_indx, subtest.shape) world_posi=(minInd_[0]-centerx) + i world_posj=(minInd_[1]-centerx) + j if(world[i][j]>subtest[minInd_]): world[i][j]=subtest[minInd_] movements[(i,j)]=(world_posi,world_posj) # if(debug): # fig,ax = plt.subplots(figsize=(10,10)) # ax.imshow(world) # ax.plot((j,world_posj),(i,world_posi)) # ax.imshow(world) # ax.scatter(j,i,c='red') # plt.show() # print(subtest[minInd_]) else: movements[(i,j)]=(i,j) return world,movements
[docs]def compute_distance_transform(target_layer, target=None): """ Distance transform algorithm implementation Parameters ---------- target_layer : nd.array An nd.array of shape (row_num,col_num). The raster layer that you need to compute the dt on. target : array-like A tupple specifying the location of the target. Default value is None Returns ------- distance_r : np.array An array of the same shape as the target_layer. Each element of the output array recieves as value, the distance to the closer element of the target_layer matrix. m : dictionary A dictionary holding a vector of the slope for each cell. """ distance_raster_zeros=np.zeros_like(target_layer) distance_raster = np.add(distance_raster_zeros, np.inf,casting="unsafe") if target!=None: distance_raster[target]=0 else: if(np.max(target_layer)>0): targetY,targetX=np.where(target_layer>0) distance_raster[targetY,targetX]=0 else: print("input target_layer does not have positive values") return movements={} distance_r,m=distance_transform(distance_raster,movements) distance_r,m=distance_transform(distance_raster,m,backward=True) return distance_r,m
[docs]def tile_layers (tile_path, rasterPath): """ """ filenames=[] layer_names=[] output_dirs=[] for filename in os.listdir(rasterPath): if filename.endswith(".tif"): # print(filename) filenames.append(os.path.join(rasterPath,filename)) layer_names.append(filename.split(".")[0]) output_dirs.append(tile_path) else: continue # print(filenames) args_ = zip(filenames,layer_names,output_dirs) print("Start Tilling...") pool=Pool(processes=4) M=pool.starmap(tileThisFile, args_) get_bbox_from_raster(filenames[0],tile_path,layer_names) print("Done")
def tileThisFile(file_name,layer_name,ouput_dir): ouput_dir_ = os.path.join(ouput_dir,layer_name) options = {'zoom': (12,16), 'resume': True,'s_srs':"EPSG:27700","webviewer":"none"} gdal2tiles.generate_tiles(file_name, ouput_dir_, **options) return "ok" def get_bbox_from_raster(file_name,ouput_dir,layer_names): src = gdal.Open(file_name) ulx, xres, xskew, uly, yskew, yres = src.GetGeoTransform() lrx = ulx + (src.RasterXSize * xres) lry = uly + (src.RasterYSize * yres) source = osr.SpatialReference() source.ImportFromEPSG(27700) # The target projection target = osr.SpatialReference() target.ImportFromEPSG(4326) # Create the transform - this can be used repeatedly transform = osr.CoordinateTransformation(source, target) # Transform the point. You can also create an ogr geometry and use the more generic `point.Transform()` NWcorner=transform.TransformPoint(ulx, uly) SEcorner=transform.TransformPoint(lrx,lry) bounds_data={"bounds":{"minx":NWcorner[0],"miny":SEcorner[1],"maxx":SEcorner[0],"maxy":NWcorner[1]},"raster_layers":layer_names} filePath = os.path.join(ouput_dir,'tile_bounds.json') with open(filePath, 'w') as f: json.dump(bounds_data, f)
[docs]def array_to_tiff(fileName,nparray,bbox): """ Burn np.array to geoTiff. Parameters ---------- References ---------- https://gis.stackexchange.com/questions/37238/writing-numpy-array-to-raster-file """ # For each pixel I know it's latitude and longitude. # As you'll see below you only really need the coordinates of # one corner, and the resolution of the file. xmin,ymin,xmax,ymax = bbox nrows,ncols = np.shape(nparray) xres = (xmax-xmin)/float(ncols) yres = (ymax-ymin)/float(nrows) geotransform=(xmin,xres,0,ymax,0, -yres) # That's (top left x, w-e pixel resolution, rotation (0 if North is up), # top left y, rotation (0 if North is up), n-s pixel resolution) # I don't know why rotation is in twice??? output_raster = gdal.GetDriverByName('GTiff').Create(fileName,ncols, nrows, 1 ,gdal.GDT_Float32) # Open the file output_raster.SetGeoTransform(geotransform) # Specify its coordinates srs = osr.SpatialReference() # Establish its coordinate encoding srs.ImportFromEPSG(27700) # This one specifies WGS84 lat long. output_raster.SetProjection( srs.ExportToWkt()) # Exports the coordinate system # to the file output_raster.GetRasterBand(1).WriteArray(nparray) # Writes my array to the raster output_raster.FlushCache()
[docs]def reproject_tiff(input_file,output_file,to_crs="4326"): """ Wrap for the GDAL.warp function Parameters ---------- input_file: str Path to input file .tif ouput_file : str Path to ouput .tif """ input_raster = gdal.Open(input_file) output_raster = output_file gdal.Warp(output_raster,input_raster,dstSRS='EPSG:'+to_crs)
[docs]def get_cut_and_fills(len_alignment,alignment_elevation,step=1): """ Parameters ---------- len_alignment : int Length of alignment alignment_elevation : array-like Elevation values of alignment. step : int Number of steps to calculate the cut/fills. Returns ------- cut,fill,tunnels,bridges : array-like Indicies where cut and fills are needed or tunnels and bridges. ydif : array-like """ ground_x = np.arange(0, len_alignment) f = interpolate.interp1d(ground_x, alignment_elevation) dif_ele=0.02*len(alignment_elevation) current_dif=abs(alignment_elevation[0]-alignment_elevation[-1]) dif_to_spent=abs(current_dif-dif_ele)/2 if alignment_elevation[0]>alignment_elevation[-1]: sign=1 else: sign=-1 origin=(0,alignment_elevation[0]-(sign*dif_to_spent)) end = (len(alignment_elevation),alignment_elevation[-1]+(sign*dif_to_spent)) # origin=(0,alignment_elevation[0]) # end=(len(alignment_elevation),alignment_elevation[-1]) fnew =interpolate.interp1d((origin[0],end[0]),(origin[1],end[1])) aligment_heights=fnew(ground_x) # set the steps to evaluate the ground surface line # step=1 # arange of x ground values x_dx = np.arange(0,len(ground_x),step) # height in the x values actual_heights=f(x_dx) aligment_heights = fnew(x_dx) # find difference ydif=actual_heights-aligment_heights # NOTICE #-------------------- # atm we save indicies, not coords...maybe we need to change that tunnels=[] bridges=[] cut=[] fill=[] # NOTICE #-------------------- # this needs to go to the config file tunnel_dist=40 bridge_gap=40 for i in range(len(ydif)): # cut if(ydif[i]>0): if(ydif[i]>tunnel_dist): tunnels.append(i) else: cut.append(i) # fill else: if(ydif[i]<-bridge_gap): bridges.append(i) else: fill.append(i) return cut,fill,tunnels,bridges,ydif,aligment_heights,actual_heights
[docs]def fill_zeros_elevation(elevation,plot=False): """ Fill elevation array averaging from the near values. Parameters ---------- elevation : array-like An array with values of elevation. plot : bool Bool for plotting the results (deafault is true) """ indexes_zero=np.where(elevation<1) if len(indexes_zero[0])>0: dif_index= np.subtract(indexes_zero[0][1:],indexes_zero[0][:-1]) break_points=np.where(dif_index>1) zero_segments=[] if(len(break_points[0])>0): start_point =0 for i in range(len(break_points[0])): break_index = indexes_zero[0][break_points[0][i]] range_=(indexes_zero[0][start_point],break_index) zero_segments.append(range_) start_point = break_points[0][i]+1 zero_segments.append((indexes_zero[0][break_points[0][-1]],indexes_zero[0][-1])) else: zero_segments.append((indexes_zero[0][0],indexes_zero[0][-1])) for i in zero_segments: if((i[0]!=0) and (i[1]!= len(elevation)-1)): prev_value= elevation[i[0]-1] next_value = elevation[i[1]+1] #we can also sample some noise here avg_elevation=(prev_value+next_value)/2 elevation[i[0]:i[1]+1]= avg_elevation elif (i[0]==0): if (i[1]+1>len(elevation)-1): elevation[0:i[1]]=elevation[i[1]] else: elevation[0:i[1]+1]=elevation[i[1]+1] elif(i[1]==len(elevation)-1): prev_value=elevation[i[0]-1] elevation[i[0]:i[1]+1]=prev_value else: print(i) if(plot): plt.plot(elevation) else: plt.savefig("elevation.png") return elevation