# -*- coding: utf-8 -*-
import tifffile as tiff
import math
import numpy
from itertools import chain
img = tiff.imread('tiffilePath')
listImg = list(chain.from_iterable(img))
#Filter out negative numbers (optional)
x = numpy.array([num for num in listImg if num >= 0])
def percentile(data, perc: int):
    size = len(data)
    return sorted(data)[int(math.ceil((size * perc) / 100)) - 1]
    
#5th percentile
print("5th percentile："+str(percentile(x, 5)))
#95th percentile
print("95th percentile："+str(percentile(x, 95)))

# -*- coding: utf-8 -*-
import arcpy
import copy
from arcpy import env
from arcpy.sa import *

# Check out the ArcGIS Spatial Analyst extension license
arcpy.CheckOutExtension("Spatial")
env.workspace = "D:\\NonSourceInundation\\result"
env.overwriteOutput = True
# original_para
dem="D:\\NonSourceInundation\\dem.tif"

# Set environment settings
cellXSzie = float(arcpy.GetRasterProperties_management(dem, "CELLSIZEX")[0])
cellYSzie = float(arcpy.GetRasterProperties_management(dem, "CELLSIZEY")[0])
maxdem=float(arcpy.GetRasterProperties_management(dem,"MAXIMUM")[0])
mindem=float(arcpy.GetRasterProperties_management(dem,"MINIMUM")[0])

def getValueSum(raster):
    sum=0
    rstArray = arcpy.RasterToNumPyArray(raster)
    rows, cols = rstArray.shape
    for rowNum in xrange(rows):
        for colNum in xrange(cols):
            v=rstArray.item(rowNum, colNum)
            sum=sum+v
    return sum

#  Calculate volume and extent of inundation based on flood levels
def getVandRaster(cellxsize,cellysize,h):
    inRaster = Raster(dem)
    deepRaster = Con(inRaster<=h, h-inRaster, 0)
    deepRaster2 = SetNull(deepRaster == 0,deepRaster)
    Vraster=deepRaster*cellxsize*cellysize
    Vraster2=Con(IsNull(Vraster),0,Vraster)
    V=getValueSum(Vraster2)
    # Set the flood amount (m3)
    if V>=1921911.658:
        deepRaster2.save("Level"+str(h)+"m_V"+str(V)+"m3.tif")
        print "Level"+str(h)+"m_V"+str(V)+"m3.tif"
        return 1
    else:
        return 0

waterLevel=mindem
# Rate of water level rise (m)
step=0.01
while waterLevel<maxdem:
    if getVandRaster(cellXSzie,cellYSzie,waterLevel)==1:
        break
    else:
        waterLevel=waterLevel+step

# -*- coding: utf-8 -*-
import arcpy
from arcpy import env
from arcpy.sa import *
import copy
arcpy.CheckOutExtension("Spatial")
env.overwriteOutput = True

def getLocationOfPointInGrid(rasterPath,pointPath,worksapceTempPath):
    env.workspace = worksapceTempPath
    cellXSzie = arcpy.GetRasterProperties_management(rasterPath, "CELLSIZEX")
    cellYSzie = arcpy.GetRasterProperties_management(rasterPath, "CELLSIZEY")
    leftPosition = arcpy.GetRasterProperties_management(rasterPath, "LEFT")
    bottomPosition = arcpy.GetRasterProperties_management(rasterPath, "BOTTOM")
    # Converting raster to array
    rstArray = arcpy.RasterToNumPyArray(arcpy.Raster(rasterPath))
    # Number of rows and columns of raster
    rows, cols = rstArray.shape
    # Create row and column number grid
    rowcell = {}
    colcell = {}
    for rowNum in xrange(rows):
        for colNum in xrange(cols):
            rowcell[(rowNum, colNum)] = [rowNum]
            colcell[(rowNum, colNum)] = [colNum]
    for rowNum in xrange(rows):
        for colNum in xrange(cols):
            rstArray[(rowNum, colNum)] = rowcell[(rowNum, colNum)][0]
    result = arcpy.NumPyArrayToRaster(rstArray, arcpy.Point(leftPosition[0], bottomPosition[0]), float(cellXSzie[0]),float(cellYSzie[0]))
    result.save("row.tif")
    dsc = arcpy.Describe(pointPath)
    coord_sys = dsc.spatialReference
    arcpy.DefineProjection_management("row.tif", coord_sys)
    for rowNum in xrange(rows):
        for colNum in xrange(cols):
            rstArray[(rowNum, colNum)] = colcell[(rowNum, colNum)][0]
    result = arcpy.NumPyArrayToRaster(rstArray, arcpy.Point(leftPosition[0], bottomPosition[0]), float(cellXSzie[0]),float(cellYSzie[0]))
    result.save("col.tif")
    dsc = arcpy.Describe(pointPath)
    coord_sys = dsc.spatialReference
    arcpy.DefineProjection_management("col.tif", coord_sys)
    # Get the row and column of the point
    inPointFeatures = pointPath
    inRasterList = [["row.tif", "row"], ["col.tif", "col"]]
    # Execute ExtractValuesTo Points
    ExtractMultiValuesToPoints(inPointFeatures, inRasterList, "")
    cursor = arcpy.SearchCursor(inPointFeatures)
    rowN = ''
    colN = ''
    for rowcol in cursor:
        rowN = int(rowcol.row)
        colN = int(rowcol.col)
    return [rowN, colN]

def  getfloodcell(waterLevel,sourcerow,sourcecol,elevationCell,rows,cols):
    elevationCellCopy = copy.deepcopy(elevationCell)
    floodbufferArray=[]
    floodbufferArray .append(elevationCellCopy[sourcerow,sourcecol])
    while(len(floodbufferArray)!=0):
        rowfirst=floodbufferArray[0][0]
        colfirst=floodbufferArray[0][1]
        elevationCellCopy[(rowfirst,colfirst)][3] = True
        floodbufferArray.pop(0)
        for row in range(rowfirst - 1, rowfirst + 2):
            for col in range(colfirst - 1, colfirst + 2):
                if(row<rows and col<cols and row>=0 and col>=0):
                    elevationinfo=elevationCellCopy[(row,col)][2]
                    if elevationinfo<=waterLevel and elevationCellCopy[(row,col)][3]==False:
                        elevationCellCopy[(row,col)][3] = True
                        floodbufferArray.append((row, col))
    return elevationCellCopy

# original_para
worksapce_path = "D:\\sourceFlood\\temp"
dem = "D:\\sourceFlood\\DEM.tif"
seed = "D:\\sourceFlood\\seed.shp"

# Convert elevation raster to array
elevationCell={}
cellXSzie = arcpy.GetRasterProperties_management(dem, "CELLSIZEX")
cellYSzie = arcpy.GetRasterProperties_management(dem, "CELLSIZEY")
leftPosition = arcpy.GetRasterProperties_management(dem, "LEFT")
topPosition = arcpy.GetRasterProperties_management(dem, "TOP")
bottomPosition = arcpy.GetRasterProperties_management(dem, "BOTTOM")
maxdem=arcpy.GetRasterProperties_management(dem,"MAXIMUM")
mindem=arcpy.GetRasterProperties_management(dem,"MINIMUM")
rstArray = arcpy.RasterToNumPyArray(arcpy.Raster(dem))

# Number of rows and columns of dem
rows, cols = rstArray.shape
print rows,cols
for rowNum in xrange(rows):
    for colNum in xrange(cols):
        # Set elevation raster pixel, initialize to no flooding
        elevationCell[(rowNum, colNum)] = [rowNum,colNum,rstArray.item(rowNum, colNum), False]

RowCol=getLocationOfPointInGrid(dem,seed,worksapce_path)
# Setting the inundation level 100
floodelevationCell=getfloodcell(100,RowCol[0],RowCol[1],elevationCell,rows,cols)
# Convert array to  raster
for rowNum in xrange(rows):
    for colNum in xrange(cols):
        rstArray[(rowNum, colNum)] = floodelevationCell[(rowNum, colNum)][3]
result = arcpy.NumPyArrayToRaster(rstArray, arcpy.Point(leftPosition[0], bottomPosition[0]), float(cellXSzie[0]),float(cellYSzie[0]))
result.save("D:\\sourceFlood\\result.tif")

# -*- coding: utf-8 -*-

import math
import arcpy
import os
from arcpy import sa
from arcpy.sa import *
from time import strftime
arcpy.CheckOutExtension("Spatial")
arcpy.env.overwriteOutput = True
print("Program start time:"+strftime("%m/%d/%Y %H:%M"))

#------------------------The parameters ---------------------------------
# The path to the parent folder of the unzipped folder of landsat 8 img
Inputpath="E://waterL8data"                             
# Projection file
projectFilePath="E://waterL8data//AEA_WGS_1984.prj"     
# Study area shp
Zone="E://waterL8data//province.shp"                  
# Temp data and results data path
arcpy.env.workspace="E://waterL8data//test"      
# Water extraction final result name (grid value of 1 for water, 2 for non-water)       
finalwaterName="finalwater.tif"                         
#------------------------------------------------------------------------------
# Get all files and folders
allFilesOrFoldersList = os.listdir(Inputpath)
for FileOrFolder in allFilesOrFoldersList:
    FileOrFolder_path = os.path.join(Inputpath, FileOrFolder)
    # Is folder
    if os.path.isdir(FileOrFolder_path):
        FolderPath=FileOrFolder_path
        allFilesList = os.listdir(FolderPath)
        for file in allFilesList:
            fileFullPath=os.path.join(FolderPath, file)
            # No.3 band
            if file[-6:] == "B3.TIF":
                green = arcpy.Raster(fileFullPath)
                # No.6 band
            elif file[-6:] == "B6.TIF":
                Mir = arcpy.Raster(fileFullPath)
                # NDWI
                MNDVI=(green-Mir)*1.0/(green+Mir)
                # NDWI correction
                MNDVIcor=Con((MNDVI>0.05) & (Mir<7200),1,2)
                # Remove the black edge
                outputPath=FileOrFolder+".tif"
                arcpy.CopyRaster_management(MNDVIcor, outputPath, "", 0, "")
                print(FolderPath+"CopyRaster success")

# Mosaic rasters
rasters=[]
for raster in arcpy.ListRasters("*.tif"):
    rasters.append(raster)
arcpy.MosaicToNewRaster_management(rasters, str(arcpy.env.workspace), "water.tif", "","8_BIT_UNSIGNED", "30", "1", "MINIMUM","FIRST")
print("Mosaic success")

water= "water.tif"
arcpy.ProjectRaster_management(water, "water_rep.tif",projectFilePath, "NEAREST", "30","#", "#", "#")
print("ProjectRaster success")

ExtractByMask("water.tif",Zone).save(finalwaterName)
print("ExtractByMask success")

print("program end time:"+strftime("%m/%d/%Y %H:%M"))

# -*- coding: utf-8 -*-
import arcpy
from arcpy import env
from arcpy.sa import *
arcpy.CheckOutExtension("Spatial")
# Workspace
env.workspace = ""
env.overwriteOutput = True
# Dem data path
dem="XXX\\xx.tif"
cellXSzie = float(arcpy.GetRasterProperties_management(dem, "CELLSIZEX")[0])
cellYSzie = float(arcpy.GetRasterProperties_management(dem, "CELLSIZEY")[0])
maxdem=float(arcpy.GetRasterProperties_management(dem,"MAXIMUM")[0])
mindem=float(arcpy.GetRasterProperties_management(dem,"MINIMUM")[0])

def getValueSum(raster):
    sum=0
    rstArray = arcpy.RasterToNumPyArray(raster)
    rows, cols = rstArray.shape
    for rowNum in xrange(rows):
        for colNum in xrange(cols):
            s=rstArray.item(rowNum, colNum)
            sum=sum+s
    return sum

# Calculation of the volume and extent of inundation based on water levels
def getVandRaster(cellxsize,cellysize,h,V_threshold):
    inRaster = Raster(dem)
    # Water depth raster
    deepRaster = Con(inRaster<=h, h-inRaster, 0)
    deepRaster2 = SetNull(deepRaster == 0,deepRaster)
    # Volume of water
    Vraster=deepRaster*cellxsize*cellysize
    Vraster2=Con(IsNull(Vraster),0,Vraster)
    # Sum
    V=getValueSum(Vraster2)
    # Threshold of volume
    if V>=V_threshold:
        deepRaster2.save("Level"+str(h)+"m_Volume"+str(V)+"m3.tif")
        return 1
    else:
        return 0

h=mindem
step=0.1
while h<maxdem:
    if getVandRaster(cellXSzie,cellYSzie,h)==1:
        break
    else:
        h=h+step

# -*- coding: utf-8 -*-
from sentinelsat.sentinel import SentinelAPI
from collections import OrderedDict

# Your Copernicus Open Access Hub Account
# Register on https://scihub.copernicus.eu/dhus/#/self-registration
name=''
password=''
api = SentinelAPI(name,password, 'https://scihub.copernicus.eu/dhus/'
xmin,ymin,xmax,ymax=120,34,121,35
roi="POLYGON(("+str(xmin)+" "+str(ymin)+","+str(xmax)+" "+str(ymin)+","+str(xmax)+" "+str(ymax)+","+str(xmin)+" "+str(ymax)+","+str(xmin)+" "+str(ymin)+"))"
# roi = 'POLYGON((DEFINE ROI))' # http://geojson.io/#map=10/42.4012/117.2488

start_date = '20190524'
end_date = '20190530'

# S1/2
platformname='Sentinel-2'
# S1：SLC(Single Look Complex)、GRD(Ground Range Detected)、OCN
# product_type = 'SLC'
# S2：S2MSI1C、S2MSI2A(S2-L2A)、S2MSI2Ap
product_type = 'S2MSI2A'
# Min, Max percentage
cloud_cover = (0, 50)
# Img save path
download_path="E:\\"

# Ordered dictionary
success_products = OrderedDict()
products = OrderedDict()

results  = api.query(area=roi,date=(start_date, end_date),platformname=platformname,producttype=product_type, cloudcoverpercentage=cloud_cover)
products.update(results)

total = len(products)
print(" {} data finded ".format(total))

def download_one(api, product, product_info):
    # download
    api.download(product, directory_path=download_path)
    # save info
    success_products[product] = product_info
    print('\t[success:] {}/{}'.format(len(success_products), total))

# Trigger offline product first
try:
    print("---Trigger offline product---")
    cnt = 1
    for product in products:
        product_odata = api.get_product_odata(product)
        if not product_odata['Online']:
            print("Trigger offline data {}： {}".format(cnt, product_odata['date']))
            # Triggering retrieval from long term archive
            api.download(product, download_path)
            cnt += 1
except:
    print("Sorry, requests for retrieval from LTA exceed user quota (No more than 20 times)")
    pass

# Start downloading
print("---Download Start---")
while len(success_products)!=total:
        for product in products:
            product_odata = api.get_product_odata(product)
            # Online
            if product_odata['Online']:
                print('Data online，can be downloaded：{} {}'.format(product_odata['date'], product_odata['title']) )
                download_one(api, product, products[product])
            else:
                print("[Data offline，please wait：] {}".format(product_odata['date'] ))

# -*- coding: utf-8 -*-
# Python version 3.X
import landsatxplore.api
from landsatxplore.earthexplorer import EarthExplorer

def request_Landsat(username,password,product,bbox,start_date,end_date,cloud_max):
    api = landsatxplore.api.API(username, password)
    scenes = api.search(
        dataset=product,
        # latitude=lat,
        # longitude=lon,
        bbox=bbox,
        start_date=start_date,
        end_date=end_date,
        max_cloud_cover=cloud_max)
    print('{} data queried.'.format(len(scenes)))
    api.logout()
    return scenes

def download_landsat(username,password,Landsat_name,output_dir):
    Earth_Down = EarthExplorer(username, password)
    for scene in Landsat_name:
        ID = scene['entityId']
        print('Downloading %s '% ID)
        Earth_Down.download(scene_id=ID, output_dir=output_dir)
    Earth_Down.logout()


if __name__ == '__main__':
    # Your USGS account
    username = ''
    password = ''
    # Select data according to sensor type
    product = 'LANDSAT_TM_C1'    # Landsat 4、5
    # product = 'LANDSAT_ETM_C1'    # Landsat 7
    # product = 'LANDSAT_8_C1'    # Landsat 8
    # lat =
    # lon =
    # The Image Coverage (xmin、ymin、xmax、ymax)
    bbox = [34, 119, 35, 120]
    start_date='1999-01-01'
    end_date='2000-01-01'
    # Maximum cloud coverage
    cloud_max = 10
    # Save path
    output_dir = 'E:\\LandsatImg\\'
    Landsat_name = request_Landsat(username,password,product,bbox,start_date,end_date,cloud_max)
    download_landsat(username,password,Landsat_name,output_dir)