<jupyter_start><jupyter_text>Function definitions<jupyter_code>## Cloud masking steps for Sentinel-2 images ##

# Define cloud mask parameters
CLOUD_FILTER = 80 # maximum percent of cloud cover in the image
CLD_PRB_THRESH = 50 # probability of clouds at pixel
NIR_DRK_THRESH = 0.15 # NIR reflectance; values below threshold are considered clouds
CLD_PRJ_DIST = 1 # maximum distance to search for cloud shadows from cloud edges
BUFFER = 50 # dilation distance from edge of cloud objects

## Creates an image collection from specified parameters of cloud probability
def get_s2_sr_cld_col(start_date, end_date):
    # Import and filter S2 SR.
    s2_sr_col = (ee.ImageCollection('COPERNICUS/S2_SR')
        .filterDate(start_date, end_date)
        .filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE', CLOUD_FILTER)))

    # Import and filter s2cloudless.
    s2_cloudless_col = (ee.ImageCollection('COPERNICUS/S2_CLOUD_PROBABILITY')
        .filterDate(start_date, end_date))

    # Join the filtered s2cloudless collection to the SR collection by the 'system:index' property.
    return ee.ImageCollection(ee.Join.saveFirst('s2cloudless').apply(**{
        'primary': s2_sr_col,
        'secondary': s2_cloudless_col,
        'condition': ee.Filter.equals(**{
            'leftField': 'system:index',
            'rightField': 'system:index'
        })
    }))

## Cloud component
## Adds the s2cloudless probability layers and derives a cloud mask
def add_cloud_bands(img):
    # Get s2cloudless image, subset the probability band.
    cld_prb = ee.Image(img.get('s2cloudless')).select('probability')

    # Condition s2cloudless by the probability threshold value.
    is_cloud = cld_prb.gt(CLD_PRB_THRESH).rename('clouds')

    # Add the cloud probability layer and cloud mask as image bands.
    return img.addBands(ee.Image([cld_prb, is_cloud]))

## Cloud shadow component
## Adds dark pixels, cloud projection, and identified shadows as bands to the image collection
def add_shadow_bands(img):
    # Identify water pixels from the SCL band.
    not_water = img.select('SCL').neq(6)

    # Identify dark NIR pixels that are not water (potential cloud shadow pixels).
    SR_BAND_SCALE = 1e4
    dark_pixels = img.select('B8').lt(NIR_DRK_THRESH*SR_BAND_SCALE).multiply(not_water).rename('dark_pixels')

    # Determine the direction to project cloud shadow from clouds (assumes UTM projection).
    shadow_azimuth = ee.Number(90).subtract(ee.Number(img.get('MEAN_SOLAR_AZIMUTH_ANGLE')));

    # Project shadows from clouds for the distance specified by the CLD_PRJ_DIST input.
    cld_proj = (img.select('clouds').directionalDistanceTransform(shadow_azimuth, CLD_PRJ_DIST*10)
        .reproject(**{'crs': img.select(0).projection(), 'scale': 100})
        .select('distance')
        .mask()
        .rename('cloud_transform'))

    # Identify the intersection of dark pixels with cloud shadow projection.
    shadows = cld_proj.multiply(dark_pixels).rename('shadows')

    # Add dark pixels, cloud projection, and identified shadows as image bands.
    return img.addBands(ee.Image([dark_pixels, cld_proj, shadows]))

## Assembles the cloud-shadow mask to produce a final masking of the images
def add_cld_shdw_mask(img):
    # Add cloud component bands.
    img_cloud = add_cloud_bands(img)

    # Add cloud shadow component bands.
    img_cloud_shadow = add_shadow_bands(img_cloud)

    # Combine cloud and shadow mask, set cloud and shadow as value 1, else 0.
    is_cld_shdw = img_cloud_shadow.select('clouds').add(img_cloud_shadow.select('shadows')).gt(0)

    # Remove small cloud-shadow patches and dilate remaining pixels by BUFFER input.
    # 20 m scale is for speed, and assumes clouds don't require 10 m precision.
    is_cld_shdw = (is_cld_shdw.focal_min(2).focal_max(BUFFER*2/20)
        .reproject(**{'crs': img.select([0]).projection(), 'scale': 20})
        .rename('cloudmask'))

    # Add the final cloud-shadow mask to the image.
    return img.addBands(is_cld_shdw)
## Converts an EarthEngine image collection into a Pandas array for a given point
## with the specified bands.
## The buffer attribute should be equal to half the spatial resolution of the final product.
## The bands should be given as a list, even for single bands.
def ee_to_df(ee_arr, lon, lat, buffer, int_limit, bands, start_date, end_date):
    # Converts columns to numeric values
    def to_numeric(dataframe, band):
        dataframe[band] = pd.to_numeric(dataframe[band], errors='coerce')
        return dataframe

    # Transform the client-side data to a dataframe
    poi = ee.Geometry.Point(lon, lat)
    arr = ee_arr.select(bands).getRegion(poi, buffer).getInfo()
    df = pd.DataFrame(arr)
    headers = df.iloc[0]
    df = pd.DataFrame(df.values[1:], columns=headers)

    # Applies the to_numeric function and fills NaN rows with interpolated values 
    for band in bands:
        df = to_numeric(df, band)
        if int_limit > 0:
            df[band].interpolate(method='linear', limit=int_limit, limit_direction='both',
                                 inplace=True)
        df.drop_duplicates(keep='first') # remove duplicates

    # Creates an index date column and drops unnecessary date, time, and coordinate columns
    df['Date'] = pd.to_datetime(df['time'], unit='ms')
    df['Date'] = df['Date'].dt.date
    df.set_index('Date', inplace=True)
    df.drop(['id', 'time', 'longitude', 'latitude'], axis=1, inplace=True)
    
    # Drop duplicate entries from the index and reindex the dataset to daily timesteps
    df = df[~df.index.duplicated()]
    df = df.reindex(pd.date_range(start=start_date, end=end_date, freq='D'))
    df.index.name = 'Date'
    
    return df

## Creates a new column that counts the number of days between precipitation events
def dry_days(df):
    dry_streak = []
    counter = 0
    for day, value in enumerate(df['total_precipitation']):
        if value > 0:
            dry_streak.append(0)
            counter = 0
        else:
            counter += 1
            dry_streak.append(counter)
    df['Dry_days'] = dry_streak
    df['Dry_days'] = df['Dry_days'].astype('float')
    return df

## Calculates and adds the vegetation index profiles from the S2 data.
## The indices are interpolated linearly 30 days in both directions.
## NDWI is based on https://doi.org/10.1016/S0034-4257(96)00067-3
def add_vegetation_index(df, name, band_1, band_2):
    df[name] = (df[band_1] - df[band_2]) / (df[band_1] + df[band_2])
    df[name].interpolate(method='linear', limit=30, limit_direction='both', inplace=True)
    return df
## Compiles a dataframes of all remotely sensed EE variables at a given point
def get_rs_data(lon, lat, point_id):
    # MODIS
    lst = ee_to_df(land_t, lon, lat, 5, 5, ['LST_Day_1km'], start_date, end_date)
    lst = lst * 0.02 - 273.15 # scaling factor and Kelvin to Celcius conversion 

    # Sentinel-1
    sigma = ee_to_df(s1, lon, lat, 5, 0, ['VV', 'VH', 'angle'], start_date, end_date)

    # Sentinel-2
    s2_vi = ee_to_df(s2, lon, lat, 5, 0, ['B4', 'B8', 'B8A', 'B11'], start_date, end_date)
    s2_vi = add_vegetation_index(s2_vi, 'NDVI', 'B4', 'B8')
    s2_vi = add_vegetation_index(s2_vi, 'NDWI', 'B8A', 'B11')
    s2_vi.drop(['B4', 'B8', 'B8A', 'B11'], axis=1, inplace=True)

    # ERA-5
    era5_weather = ee_to_df(era5, lon, lat, 5, 0, ['mean_2m_air_temperature', 'total_precipitation'],
                 start_date, end_date)
    era5_weather = dry_days(era5_weather)
    era5_weather['mean_2m_air_temperature'] = era5_weather['mean_2m_air_temperature'] - 273.15

    # Merge dataframes
    dfs = [lst, sigma, s2_vi, era5_weather]
    df = reduce(lambda  left,right: pd.merge(left, right, on=['Date'], how='outer'), dfs)
    df['ID'] = point_id
    
    return df


## Creates a dataframe from ISMN data 
## Using the ISMN reader https://github.com/TUW-GEO/ismn
def get_sm_data(network, depth_start, depth_end, snow_depth=False):
    # Retrieves the name of each station in the network
    stations = {}
    for n, station in enumerate(sm_data[network]):
        stations[n] = station.metadata['station'][1]
    
    # Creates a new dataframe with the station IDs and coordinates
    grid = sm_data.collection.grid
    gpis, lon, lat = grid.get_grid_points()
    df_coords = pd.DataFrame(index=pd.Index(gpis, name='ID'), data={'longitude': lon,
                                                                    'latitude': lat})
    df_coords = df_coords.rename(index=stations)

    # Writes soil moisture at 0-5cm depth to a dataframe
    sm_df = pd.DataFrame()
    for network, station, sensor in sm_data.collection.iter_sensors(
        variable='soil_moisture', depth=[depth_start, depth_end]):
        
        data = sensor.read_data()
        sensor_df = pd.DataFrame(data)
        sensor_df = sensor_df[sensor_df['soil_moisture_flag'] == 'G'] # filters out bad readings
        sensor_df = sensor_df.resample('D').mean() # resamples to daily values
        sensor_df['Date'] = sensor_df.index.date
        sensor_df['ID'] = station.metadata['station'][1]
        sm_df = sm_df.append(sensor_df)
    
    # Writes snow depth to a dataframe if specified and filters out days with snow cover
    if snow_depth:
        snow_df = pd.DataFrame()
        for network, station, sensor in sm_data.collection.iter_sensors(
            variable='snow_depth', depth=[0., 0.]):
            
            data = sensor.read_data()
            sensor_df = pd.DataFrame(data)
            sensor_df.loc[sensor_df['snow_depth']<=0, 'snow'] = 'yes'
            sensor_df = sensor_df.resample('D').count() # counts the entries with snow cover
            sensor_df['Date'] = sensor_df.index.date
            sensor_df['ID'] = station.metadata['station'][1]
            sensor_df = sensor_df[['ID', 'Date', 'snow']]
            snow_df = snow_df.append(sensor_df)
        try:
            df = pd.merge(sm_df, snow_df, how='left', left_on=['Date', 'ID'], right_on=['Date', 'ID'])
            df = df[df['snow']==0] # any day with snow cover entries is discarded
        except KeyError:
            print('No snow depth variable found, set snow_depth=False and run again')
    else:
        df = sm_df
        
    df.rename({'Date_x' : 'Date'}, axis=1, inplace=True)
    df.set_index('Date', inplace=True)
    df = df[['ID', 'soil_moisture']]
    print('ISMN data successfully translated to dataframe')
    return df, df_coords


## Queries the Earth Engine data for each soil moisture probe (lat, lon coordinates) in the
## dataset and merges the EE data with the sm data
def merge_sm_rs_data(soil_moisture_df, coordinate_df):
    df = pd.DataFrame()
    print('Extracting data from Earth Engine...')
    i_start, i_end = 1, len(coordinate_df)
    for station, coords in coordinate_df.iterrows():
        print('Station {} of {}'.format(i_start, i_end))
        df_sm_temp = soil_moisture_df[soil_moisture_df['ID'] == station]
        df_rs_temp = (get_rs_data(coords['longitude'], coords['latitude'], station))
        df_rs_temp = df_rs_temp.reindex(pd.date_range(start=start_date, end=end_date, freq='D'))
        df_temp = pd.merge(df_rs_temp, df_sm_temp, how='left', left_index=True, right_index=True)
        df_temp.rename({'ID_x' : 'ID'}, axis=1, inplace=True)
        df_temp.drop('ID_y', axis=1, inplace=True)
        df = df.append(df_temp)
        i_start += 1
    df.index.name = 'Date'
    return df<jupyter_output><empty_output><jupyter_text>Data specifications<jupyter_code>#!# Important to verify that all datasets contain data between the two periods #!#
#!# The current timeframe between start and end are the lower and upper limits, #!#
#!# defined by the start of Sentinel-2 SR data and end of ERA-5
start_date, end_date = '2017-03-28', '2020-07-09' 

# Sentinel-1 GRD C-band SAR data
# https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S1_GRD
# Values are expressed in decibel on a logarithmic scale
s1 = (ee.ImageCollection("COPERNICUS/S1_GRD")
      .filterDate(ee.Date(start_date), ee.Date(end_date))
      .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
      .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VH'))
      .filter(ee.Filter.eq('instrumentMode', 'IW')))

# Sentinel-2 Level-2A surface reflectance
# https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2_SR
s2_sr_cld_col_eval = get_s2_sr_cld_col(start_date, end_date)
s2 = s2_sr_cld_col_eval.map(add_cld_shdw_mask)

# MODIS daily land surface temperature
# https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MOD11A1
# Temperature is given in Kelvin with a scaling factor of 0.02
land_t = (ee.ImageCollection("MODIS/006/MOD11A1")
          .filterDate(start_date, end_date)
          .select('LST_Day_1km'))

# ERA-5 daily aggregate re-analysis data
# https://developers.google.com/earth-engine/datasets/catalog/ECMWF_ERA5_DAILY
# Maximum temperature (air) uses Kelvin
# Precipitation is in meters
era5 = (ee.ImageCollection('ECMWF/ERA5/DAILY')
        .filterDate(start_date, end_date)
        .select(['mean_2m_air_temperature', 'total_precipitation']))
# Soil moisture data retrieved from ISMN
# Retrieved from https://ismn.geo.tuwien.ac.at/en/
root = r'C:\Users\Nick\Documents\_Master_Thesis\Soil_Moisture\ISMN'
sm_data = ISMN_Interface(os.path.join(root, 'REMEDHUS.zip'), parallel=True)<jupyter_output><empty_output><jupyter_text>Dataframe compiler<jupyter_code>## Compile the final merged dataframe of RS and SM data

# Specify the probe network name
# For multiple networks, these must be iterated through either through a script or by saving
# a .csv for each network and then later merged
network = 'SNOTEL'

# Creates a soil moisture dataframe from the ISMN data
# Specify the depth intervals (in meters) of the moisture measurements, e.g. 0.0, 0.05 for soil
# moisture between 0-5cm
# Set snow_depth=True to filter out days with snow cover if the data contains a snow depth variable
df_sm, df_coords = get_sm_data(network, 0.0, 0.05, snow_depth=True)

# Merge the soil moisture and remote sensing data for each station
df = merge_sm_rs_data(df_sm, df_coords)

# Export the dataframe to a local disk
# Specify the out_dir with your local path
out_dir = r'C:\Users\Nick\Documents\_Master_Thesis\Soil_Moisture\ISMN'
df.to_csv(os.path.join(out_dir, 'SM_Data_{}.csv'.format(network)))
print('Done!')<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup