magassagvizsgalat=0 #0-AVERAGE WATER SURFACE ELEVATION, 1- mOST FREQUENT VALUE OF ELEVATION (ON THE LAKE SURFACE)


to="SOMESHAPEFILENAME.shp" #INSERT HERE THE SHAPEFILE NAME, WHICH CONTAINS THE LAKE POLYGONS, YOU MAY NEED A PRJ PROJECTION FILE
import ogr, osr
driver=ogr.GetDriverByName('ESRI Shapefile') 
to1=driver.Open(to,1)
layer=to1.GetLayer() 
numFeatures=layer.GetFeatureCount()


import gdal
from gdalconst import *
raszter="SOMEGEOTIFF.tif" #YOUR TERRAIN MODEL IN GEOTIFF FORMAT (32 BIT FLOATING POINT)
dataset=gdal.Open(raszter,GA_ReadOnly)
tx=dataset.RasterXSize
ty=dataset.RasterYSize
cs=abs(dataset.GetGeoTransform()[1])
cx=dataset.GetGeoTransform()[3]
cy=dataset.GetGeoTransform()[0]
t=dataset.ReadAsArray()
ujt=t

bfx=cx
bfy=cy
tile=[]
i=0

for k in range(0,tx,100):
    for l in range (0,ty,100):
        bfx=cx-k*cs
        bfy=cy+l*cs
        jay=bfy+cs*100
        jax=bfx-cs*100
        tile.append([i,bfx,bfy,jax,jay])
        i+=1
print (tile)
for i in range (0,numFeatures): 
    print (i) 
    feat=layer.GetNextFeature()
    geom=feat.GetGeometryRef()
    #Vetületi transzformáció EOTRbe
    '''sourceSR = layer.GetSpatialRef()
    targetSR = osr.SpatialReference()
    targetSR.ImportFromEPSG(23700) # EOTR
    coordTrans = osr.CoordinateTransformation(sourceSR,targetSR)
    geom.Transform(coordTrans)
    ter=geom.GetArea() #Vízfelület m2-ben
    coordTrans2 = osr.CoordinateTransformation(targetSR,sourceSR)
    geom.Transform(coordTrans2)'''
    #--
    #A felület befoglalója
    yk=geom.GetEnvelope()[0]
    yn=geom.GetEnvelope()[1]
    xk=geom.GetEnvelope()[2]
    xn=geom.GetEnvelope()[3]
    print(xk)
    print(xn)
    print(yk)
    print(yn)
    for j in range (0,len(tile)): #Melyik tile-ra esik az objektum befoglaló négyszöge
        if (yk>=tile[j][2] and yk<=tile[j][4]):
            tn1=tile[j][0]
        if (yn>=tile[j][2] and yn<=tile[j][4]):
            tn2=tile[j][0]
        if xk>=tile[j][3] and xk<=tile[j][1]:
            tn3=tile[j][0]
        if xn>=tile[j][3] and xn<=tile[j][1]:
            tn4=tile[j][0]
    
    k1y=round((tile[tn1][2]-cy)/(cs)) #LEFTy
    k2y=round((tile[tn2][4]-cy)/(cs)) #RIGHTy
    k1x=round((tile[tn4][3]-cx)/(cs))*-1 #UPPERx
    k2x=round((tile[tn3][1]-cx)/(cs))*-1#LOWERx
    
    print (k1y)
    print (k2y)
    print (k1x)
    print (k2x)
    if k1y==k2y:
        k2y=k1y+100
    if k1x==k2x:
        k2x=k1x+100
    if (k2x>cx):
        k2x=ty
    if k2y>cy:
        k2y=tx
    if k1y!=0:
        k1y=k1y-100
    if k1x!=0:
        k1x=k1x-100 
    if k2y!=ty:
        k2y=k2y+100
    if k2x!=tx:
        k2x=k2x+100 
    magatlag_p=0
    darabp=0
    gy=[]     
    #NEW LAKE SURFACE HEIGHT    
    for j in range(k1x,k2x):
        for k in range(k1y,k2y):
            x=cx-cs*j-cs/2
            y=cy+k*cs+cs/2
            pg=ogr.Geometry(ogr.wkbPoint)
            pg.AddPoint(y,x)

            
            if geom.Contains(pg) and t[j][k]>0 and magassagvizsgalat==0: 
                magatlag_p=magatlag_p+t[j][k]
                darabp+=1
            if geom.Contains(pg) and t[j][k]>0 and magassagvizsgalat==1: 
                gy.append(round(t[j][k]))
    
        
           
    if magassagvizsgalat==1:
        dt = {x:gy.count(x) for x in gy}
        mag=max(dt,key=dt.get)
    if darabp!=0 and magassagvizsgalat==0:
        mag=round(magatlag_p/darabp,3)
    for j in range(k1y,k2y):
        
        for k in range(k1x,k2x):
            x=cx-cs*j-cs/2
            y=cy+k*cs+cs/2
            pg=ogr.Geometry(ogr.wkbPoint)
            pg.AddPoint(y,x)
            if geom.Contains(pg): 
                ujt[j][k]=mag
                
to1.Destroy()

fn=open("NEW.asc",'w') #NEW DIGITAL ELEVATION MODEL
for i in range(0,len(ujt)):
    if i==0:
        fn.write('ncols '+str(tx)+'\n'+'nrows '+str(ty)+'\n'+'xllcorner '+str(cy)+'\n'+'yllcorner '+str(cx-cs*ty)+'\n'+'cellsize '+str(cs)+'\n')    
    for j in range(0,len(ujt[i])):
        fn.write(str(ujt[i][j])+' ')
    fn.write('\n')
        
fn.close()