gee_colocate_spectra.py

import ee
import math
import numpy as np
import datetime
import pandas as pd
#from atmospheric import Atmospheric

ee.Initialize()

# global variables using in functions
JRC_STARTTIME = ee.Number(ee.Date('1984-03-17').millis()) # start time of jrc availability
JRC_ENDTIME = ee.Number(ee.Date('2015-10-17').millis()) # end time of jrc availability

# helper function to unpack bit values
def getQABits(image, start, end, newName):
    #Compute the bits we need to extract.
    pattern = 0;
    for i in range(start,end+1):
       pattern += math.pow(2, i);

    # Return a single band image of the extracted QA bits, giving the band
    # a new name.
    return image.select([0], [newName])\
                  .bitwiseAnd(int(pattern))\
                  .rightShift(start);

#function to mask clouds in landsat imagery using the pixel_qa band
def maskClouds(image):
    qaClouds = getQABits(image.select('pixel_qa'),5,5,'qa').neq(1)
    qaShadows = getQABits(image.select('pixel_qa'),3,3,'qa').neq(1)
    mask = qaClouds.And(qaShadows)
    noClouds = image.updateMask(mask).set("system:time_start",image.get("system:time_start"))
    return noClouds

# function to mask land areas in Landsat imagery
def maskLand(image):
    # get the jrc mask if within jrc time period
    def jrcMask(image):
        date = ee.Date(image.get('system:time_start'))
        blank = ee.Image(0);
        jrc = ee.ImageCollection('JRC/GSW1_0/MonthlyHistory')
        waterClass = ee.Image(jrc.select('water').filterDate(date.advance(-30,'day'),date.advance(30,'day')).max())
        water = blank.where(waterClass.eq(2),1)
        return water

    # get internal landsat water mask if outside of jrc time period
    def qaMask(image):
        qa = self.getQABits(image.select('pixel_qa'),2,2,'qa');
        water = qa.eq(1);
        return water

    # get image date and select which water mask to use based on availability
    imgDate = ee.Number(image.get('system:time_start'))
    result = ee.Image(
                ee.Algorithms.If(imgDate.gt(JRC_STARTTIME).And(imgDate.lt(JRC_ENDTIME)),
                jrcMask(image),
                qaMask(image)
                ))

    return result


#  read in station location and water quality data
statLoc = pd.read_csv("/Users/kmarkert/Documents/Kel/research/mekong_tss/mrcWQData/statID_local.csv",low_memory=False)
wqtable = pd.read_csv("/Users/kmarkert/Documents/Kel/research/mekong_tss/mrcWQData/mrc_wq_dataTable.csv",low_memory=False)

# extract out the geographic location for all stations
statNames = np.array(statLoc.StatID)
stat_X = np.array(statLoc.X_Deg)
stat_Y = np.array(statLoc.Y_Deg)

# extract station id, time, and tss values
statWq = np.array(wqtable.StatID)
wqTime = np.array(wqtable.SDate)
wqTss = np.array(wqtable.TSS_mgL)

# From: https:#landsat.usgs.gov/what-are-best-spectral-bands-use-my-study
wavelengths = [0.485, 0.56, 0.66, 0.83, 1.65, 2.22]
# landsat band names for LT5 and LE7 (will be different for LC8)
bands = ['B1','B2','B3','B4','B5','B7']

# define empty lists to populate with values from ee requests
b1 = []
b2 = []
b3 = []
b4 = []
b5 = []
b7 = []
scnid = []
statid = []
locx = []
locy = []
tss = []
date = []
acqtime = []
distance = []

dem = ee.Image('USGS/SRTMGL1_003')# Shuttle Radar Topography mission covers *most* of the Earth

senCode = {'LT4':'LT04','LT5':'LT05','LE7':'LE07'}

# sensor 
lt4 = ee.ImageCollection('LANDSAT/LT04/C01/T1_SR')
lt5 = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR')
le7 = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR')


jrc = ee.ImageCollection('JRC/GSW1_0/MonthlyHistory')
blank = ee.Image(1)

wetSeason = np.arange(6,12)

collection = ee.ImageCollection(lt4.merge(lt5).merge(le7)).map(maskClouds)

for i in range(statWq.size):
    
    timeOut = wqTime[i].split('/')
#    if int(timeOut[2]) > 20:
#        timeOut[2] = '19' + timeOut[2]
#    else:
#        timeOut[2] = '20' + timeOut[2]
        
    try:
    
        scnTime = datetime.datetime(int(timeOut[2]),int(timeOut[1]),int(timeOut[0]))
    
        idx = np.where(statWq[i].capitalize() == statNames)
    
        geom = ee.Geometry.Point([stat_X[idx][0], stat_Y[idx][0]])
    
        station = ee.Feature(geom, {'label': 'WaterQualityStation'})
    
        statPoints = ee.FeatureCollection([station])

        dt = ee.Date('{0}'.format(scnTime.strftime('%Y-%m-%d')))
        
        waterClass = ee.Image(jrc.select('water').filterDate(dt.advance(-15,'day'),dt.advance(30,'day')).max());
        water = blank.updateMask(waterClass.neq(2));
            
        dist = blank.cumulativeCost(water,1000).rename('dist')
        
        if scnTime.month in wetSeason:
            offset = 2 # use two day window in wet season to overcome cloud issues
        else:
            offset = 1 

        inImg = ee.Image(collection.filterDate(dt.advance(-offset,'day'),dt.advance(offset,'day')).filterBounds(geom).first())

        acqt = scnTime.strftime('%Y-%m-%d')

        lsMekong = maskLand(inImg).select(bands).addBands(dist)

        result = lsMekong.reduceRegion(ee.Reducer.mean(), statPoints, 30)

    except:
        continue

    try:
        tssval = wqTss[i]

        if np.isnan(tssval) == False:
            
            data = result.getInfo()
            #print data

            b1.append(data['B1'])
            b2.append(data['B2'])
            b3.append(data['B3'])
            b4.append(data['B4'])
            b5.append(data['B5'])
            b7.append(data['B7'])
            distance.append(data['dist'])
            acqtime.append(acqt)

            statid.append(statWq[i].capitalize())
            date.append(wqTime[i])
            locx.append(stat_X[idx][0])
            locy.append(stat_Y[idx][0])

            tss.append(tssval)

            print('B1:{0}\tB2:{1}\tB3:{2}\tB4:{3}\tB5:{4}\tB7:{5}\tTSS:{6}\tELV:{7}\tTime:{8}\tOZN:{9}\tWV:{10}\tAOT:{11}\tSZ:{12}'.format(
                    len(b1),len(b2),len(b3),len(b4),len(b5),len(b7),len(tss),len(elv),len(acqtime),len(ozn),len(wv),len(aod),len(sz)))
            
            
            data = {'B1':b1,'B2':b2,'B3':b3,'B4':b4,'B5':b5,'B7':b7,'SceneID':scnid,'StatID':statid,
                    'lat':locy,'lon':locx,'Date':date,'TSS_mgL':tss,'Elv':elv,'AcqTime':acqtime,
                    'Ozone':ozn,'WaterVapor':wv,'SolarZenith':sz,'AOT':aod,'Distance':distance}
            
            df = pd.DataFrame(data)
            df.to_csv('gee_spectra.csv')
            
            df = pd.DataFrame(data)
            
            df.to_csv('/Users/kmarkert/Documents/Kel/research/mekong_tss/ls_sr_tss_coincidence_data.csv')

        else:
            continue

    except:
        pass


#