1. gee examples
2. gee Beginner's Cookbook

Defining the GEO-function

#################################### 
######### computing pixel offsets
####################################
#  Need to get pixel offsets from the upper
#  left corner for specific coordinates x,y
# index from 0 to n-1

def get_xy(coordinate, ori_x, ori_y, pixel_width, pixel_height):
    """
    get image x,y index by GIS coordinate 
    
    #########
    parameters:
    
    coordinate: pandas DataFrame, (Latitude - y, Longitude - x)
    
    ori_x: geotransform[0],  left-upper longitude in GEOTiff, origin x coordinate
    
    ori_y: geotransform[3],  left-upper latitude in GEOTiff, origin y coordinate
    
    pixel_width: geotransform[1], pixel width
    
    pixel_height: geotransform[5], pixel height (negative)
    
    #########
    return:
    
    coordinate_xy: pandas, (x, y)
    """
    x_offset = ((coordinate['lon'].to_numpy() - ori_x) / pixel_width).astype(int)
    y_offset = ((coordinate['lat'].to_numpy() - ori_y) / pixel_height).astype(int)
    coordinate_xy = coordinate.copy()
    coordinate_xy['x'] = x_offset
    coordinate_xy['y'] = y_offset
    
    return coordinate_xy[['lat', 'lon', 'x', 'y']]
def get_value(data, xy):
    
    """
    """
    
    xy_copy = xy.copy()
    value_list = []
    y, x = xy['y'].to_numpy(), xy['x'].to_numpy()
    for i in range(len(x)):
        #value = np.mean(data[y[i]-1:y[i]+2 , x[i]-1: x[i] + 2])
        value = data[y[i], x[i]]
        value_list.append(value)
    xy_copy['value'] = value_list
    return xy_copy[['lat', 'lon', 'x', 'y','value']]
def GEO_resample(data, x, y, kernel = 3):
    minus = kernel//2
    return np.mean(data[y - minus: x + kernel - minus , x - minus: x + kernel - minus ])

def make_raster(in_ds, fn, data, data_type, nodata=None):
    """Create a one-band GeoTIFF.
    in_ds - datasource to copy projection and geotransform from
    fn - path to the file to create
    data - NumPy array containing data to write
    data_type - output data type
    nodata - optional NoData value
    """
    driver = gdal.GetDriverByName('GTiff')
    out_ds = driver.Create(fn, in_ds.RasterXSize, in_ds.RasterYSize, 1, data_type)
    out_ds.SetProjection(in_ds.GetProjection())
    out_ds.SetGeoTransform(in_ds.GetGeoTransform())
    out_band = out_ds.GetRasterBand(1)
    
    if nodata is not None:
        out_band.SetNoDataValue(nodata)
    
    out_band.WriteArray(data)
    out_band.FlushCache()
    out_band.ComputeStatistics(False)
    
    return out_ds