# -----------------------------------------------------------------
# ██████╗ █████╗ ██╗██████╗ ███╗ ██╗██████╗
# ██╔══██╗██╔══██╗██║██╔══██╗ ████╗ ██║██╔══██╗
# ██████╔╝███████║██║██║ ██║ ██╔██╗ ██║██████╔╝
# ██╔══██╗██╔══██║██║██║ ██║ ██║╚██╗██║██╔══██╗
# ██║ ██║██║ ██║██║██████╔╝ ██║ ╚████║██║ ██║
# ╚═╝ ╚═╝╚═╝ ╚═╝╚═╝╚═════╝ ╚═╝ ╚═══╝╚═╝ ╚═╝
# -----------------------------------------------------------------
#-----------------------------------------------------
# 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 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