painesville

import ee
import numpy as np
import matplotlib.pyplot as plt
import shapefile
import pandas as pd


#ee.Authenticate(force=True)
# Token generated with all permissions 29 May 2024, DMK
ee.Authenticate()
ee.Initialize(project = 'ee-fortschthomas52')

sf = shapefile.Reader("/Users/hydro3/Documents/USGSSites/shape_files/ohio/painesville.shp")
shapes = sf.shapes()
points = shapes[0].points
aoi = ee.Geometry.Polygon(points)

start = pd.date_range(start= '2017-12-05' , end='2017-12-22' , 
              freq='5d')
end = pd.date_range(start='2017-12-09' , end='2017-12-26' , 
              freq='5d')
dates = pd.DataFrame ({'start': start ,  'end': end})

def pulldata(startDate, endDate):
    # Define source data
    image = ee.ImageCollection('COPERNICUS/S2_SR_HARMONIZED')\
        .filterDate(startDate, endDate)\
        .select('B2', 'B3', 'B4', 'B8', 'B11')
         #.filterMetadata('CLOUDY_PIXEL_PERCENTAGE', 'less_than', 10)\


# CRS is not the same.
# proj = image.first().select('B2').projection() # EPSG:32656, UTM zone 56N (Siberia?)
# proj = balule.projection() # EPSG:4326

# Export arrays
# https://gist.github.com/jdbcode/f4d56d72f7fc5beeaa3859999b1f5c3d
# https://gist.github.com/jdbcode/f4d56d72f7fc5beeaa3859999b1f5c3d?permalink_comment_id=3355627#gistcomment-3355627
    mosaic = image.median().reproject(crs='EPSG:32617', scale=10) # This allows us to set the resolution.
    band_arrs = mosaic.sampleRectangle(region=aoi, defaultValue = 0)

# Get individual band arrays.
    band_arr_b2 = band_arrs.get('B2')   # Blue
    band_arr_b3 = band_arrs.get('B3')   # Green
    band_arr_b4 = band_arrs.get('B4')   # Red
    band_arr_b8 = band_arrs.get('B8')   # NIR
    band_arr_b11 = band_arrs.get('B11') # SWIR

# Transfer the arrays from server to client and cast as np array.
    b2 = np.array(band_arr_b2.getInfo())   # b2  Blue
    b3 = np.array(band_arr_b3.getInfo())   # b3  Green
    b4 = np.array(band_arr_b4.getInfo())   # b4  Red
    b8 = np.array(band_arr_b8.getInfo())   # b8  NIR
    b11 = np.array(band_arr_b11.getInfo()) # b11 SWIR
    if(np.max(b2)!=0):
        
    #np_arr_b4 = np.expand_dims(b4, 2)
    #np_arr_b3 = np.expand_dims(b3, 2)
    #np_arr_b2 = np.expand_dims(b2, 2)

    #rgb_img = np.concatenate((np_arr_b4, np_arr_b3, np_arr_b2), 2)

    #rgb_img = (255*((rgb_img)/3000)).astype('uint8')
    #plt.imshow(rgb_img)
    #plt.show()

# Normalized Difference Water Index (NDWI) 
# NDWI = ( G - NIR ) / ( G + NIR )
        ndwiG = (b3-b8)/(b3+b8) # Gao
# NDWI = ( NIR - SWIR ) / ( NIR + SWIR )
        ndwiM = (b8-b11)/(b8+b11) # McFeeters
# NDWI = ( G - SWIR ) / ( G + SWIR )
        mndwi = (b3-b11)/(b3+b11) # Modified NDWI
        water = ndwiG > -0.02

        TSS1 = np.NAN
        Secchi1 = np.NAN
        TSS2 = np.NAN
        Secchi2 = np.NAN
        TSS3 = np.NAN
        Secchi3 = np.NAN
        TSS4 = np.NAN
        Ratio = np.NAN



        TSS1 = (b3+b4)/2
        TSS1 = TSS1 * water
    #plt.imshow(TSS1 , cmap= "summer")
    #plt.colorbar(orientation="vertical")
    #plt.show
        TSS1 = np.sum(TSS1) / np.sum(TSS1>0)
    #print(TSS1)

        Secchi1 = (b2/b4)
        Secchi1 = Secchi1 * water
        Secchi1 = np.sum(Secchi1)/ np.sum(Secchi1>0)

        TSS2 = (b3/b4)
        TSS2 = TSS2 * water
        TSS2 = np.sum(TSS2) / np.sum(TSS2>0)

        Secchi2 = (b4/b3)
        Secchi2 = Secchi2 * water
        Secchi2 = np.sum(Secchi2)/ np.sum(Secchi2>0)

        TSS3 = (b8/b3 , b8/b4)
        TSS3 = TSS3 * water
        TSS3 = np.sum(TSS3) / np.sum(TSS3>0)

        Secchi3 = (b4/b2)+b2
        Secchi3 = Secchi3 * water
        Secchi3 = np.sum(Secchi3)/ np.sum(Secchi3>0)

        TSS4 = (b4/b3)+b8
        TSS4 = TSS4 * water
        TSS4 = np.sum(TSS4) / np.sum(TSS4>0)

        Ratio = (ndwiG/ndwiM)
        Ratio = Ratio * water
        Ratio = np.sum(Ratio)/ np.sum(Ratio>0) 


        f = open("Painesville5day.txt", "a")
        f.write(str(startDate) + ", " + str(endDate) + ", " + 
                str(TSS1) + ", " + str(Secchi1) + ", " + 
                str(TSS2) + ", " + str(Secchi2) + ", " + 
                str(TSS3) + ", " + str(Secchi3) + ", " + 
                str(TSS4) + "," + str(Ratio) + '\n') 
            
        f.close()

for i in range(len(dates)):
    pulldata(dates["start"][i], dates["end"][i])