import random
from osgeo import ogr,gdal

InputImage = r'C:UsersUsuarioDesktopQGIS-PoolOrtofotosPNOA_MA_OF_ETRS89_HU31_h50_0823.tif'
Shapefile = r'C:UsersUsuarioDesktopQGIS-PoolGitPool_shapefilePool.shp'
RasterFormat = 'GTiff'
PixelRes = 0.5
VectorFormat = 'ESRI Shapefile'
folder = 'C:\Users\Usuario\Desktop\QGIS-Pool\Git\images\'
tilemapping = {}

# Open datasets
Raster = gdal.Open(InputImage, gdal.GA_ReadOnly)
Projection = Raster.GetProjectionRef()

VectorDriver = ogr.GetDriverByName(VectorFormat)
VectorDataset = VectorDriver.Open(Shapefile, 0) # 0=Read-only, 1=Read-Write
layer = VectorDataset.GetLayer()
FeatureCount = layer.GetFeatureCount()
print("Feature Count:",FeatureCount)

# Iterate through the shapefile features
Count = 0
for feature in layer:
    Count += 1
    #print("Processing feature "+str(Count)+" of "+str(FeatureCount)+"...")

    geom = feature.GetGeometryRef()
    xmid, xmid, ymid, ymid = geom.GetEnvelope()
    xmid = xmid + random.randint(-40, 40)
    ymid = ymid + random.randint(-40, 40)
    xmin, ymin = xmid - 50, ymid - 50
    xmax, ymax = xmid + 50, ymid + 50

    # Create raster
    OutTileName = folder + str(Count) + '.tif'
    OutTile = gdal.Translate(OutTileName, Raster, width=224, height=224, projWin = [xmin, ymax, xmax, ymin])
    OutTile = None # Close dataset

    centers = []
    for feature1 in layer:
        geom1 = feature1.GetGeometryRef()
        x, x, y, y = geom1.GetEnvelope()
        if (xmin < x) & (x < xmax) & (y > ymin) & (y < ymax):
            x = (x - xmin) / (xmax - xmin)
            y = (y - ymin) / (ymax - ymin)
            centers.append((x,y))


    if len(centers) > 0:
        tilemapping[Count] = centers


# Close datasets
Raster = None
VectorDataset.Destroy()
print("Done.")