from PIL import Image
import requests
from sat7_pointer import *

# Load Landsat 7 band 1, 2, 3 TIF images and create a composite image
# from the three bands.
#
# The composite image is a greyscale image with a red, green and blue
# channel.
#
# The red channel is the red channel of band 3, the green channel is the
# green channel of band 2 and the blue channel is the blue channel of band 1.
#
# The composite image is saved as a PNG file.

band1 = Image.open('LE07_L1TP_177025_20210723_20210818_02_T1_B1.TIF')
band2 = Image.open('LE07_L1TP_177025_20210723_20210818_02_T1_B2.TIF')
band3 = Image.open('LE07_L1TP_177025_20210723_20210818_02_T1_B3.TIF')

composite = Image.merge('RGB', (band3, band2, band1))

# Load corner coordinates of the image.
#
# The coordinates are stored in a MTL text file.

mtl_data = load_metadata('LE07_L1TP_177025_20210723_20210818_02_T1_MTL.txt')

# Fetch coordinates of a city.
#
# The coordinates of a city are fetched from the OpenStreetMap API.
#
# The coordinates are used to calculate the corner coordinates of the
# composite image.
#
# # City is Belgorod, Russia.

url = 'http://nominatim.openstreetmap.org/search?q=Belgorod,+Russia&format=json'
response = requests.get(url)
data = response.json()

# Convert bounding box coordinates to image coordinates.

x0, y1 = lat_lot_to_pixel(data[0]['boundingbox'][0], data[0]['boundingbox'][2], mtl_data)
x1, y0 = lat_lot_to_pixel(data[0]['boundingbox'][1], data[0]['boundingbox'][3], mtl_data)

print(x0, y0)
print(x1, y1)
print(composite.size)
# Crop the image to the bounding box coordinates.

cropped = composite.crop((x0, y0, x1, y1))

cropped.save('cropped.png')

###############################################################################
# LAB 2
###############################################################################

# Load landsat band 4 TIF image.
# Band 4 is near infrared.

band4 = Image.open('LE07_L1TP_177025_20210723_20210818_02_T1_B4.TIF')

# Claculate the NDVI.
#
# The NDVI is calculated by dividing
# difference of the red and near infrared band 4 by the sum of the
# near infrared and near infrared band 4.
#
# The NDVI is calculated for each pixel in the cropped image.
#
# The NDVI is saved as a PNG image.

ndvi = band4.copy().crop((x0, y0, x1, y1))
red = band3.copy().crop((x0, y0, x1, y1))
result = cropped.copy()

print(red.getpixel((0, 0)))
print(ndvi.getpixel((0, 0)))

# Colors for the NDVI.
def get_color(value):
    if value < 0.033:
        return (255, 255, 255)
    elif value < 0.066:
        return (196, 184, 168)
    elif value < 0.1:
        return (180, 150, 108)
    elif value < 0.133:
        return (164, 130, 76)
    elif value < 0.166:
        return (148, 114, 60)
    elif value < 0.2:
        return (124, 158, 44)
    elif value < 0.25:
        return (148, 182, 20)
    elif value < 0.3:
        return (116, 170, 4)
    elif value < 0.35:
        return (100, 162, 4)
    elif value < 0.4:
        return (84, 150, 4)
    elif value < 0.45:
        return (60, 134, 4)
    elif value < 0.5:
        return (28, 114, 4)
    elif value < 0.6:
        return (4, 96, 4)
    elif value < 0.7:
        return (4, 74, 4)
    elif value < 0.8:
        return (4, 56, 4)
    elif value < 0.9:
        return (4, 40, 4)
    else:
        return (0, 0, 0)

for x in range(ndvi.size[0]):
    for y in range(ndvi.size[1]):
        r = red.getpixel((x, y))
        nir = ndvi.getpixel((x, y))
        if nir + r == 0:
            result.putpixel((x, y), get_color(0))
        else:
            result.putpixel((x, y), get_color((nir - r) / (nir + r)))

result.save('ndvi.png')