160 lines
4.7 KiB
Python
160 lines
4.7 KiB
Python
from PIL import Image
|
||
import requests
|
||
from sat7_pointer import *
|
||
from numpy import np
|
||
|
||
# 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')
|
||
|
||
|
||
# Calculate FAPAR (Fraction of Absorbed Photosynthetically Active Radiation)
|
||
# for each pixel in the cropped image.
|
||
#
|
||
# Bands 1, 3, 4 are used.
|
||
|
||
# Coefficients of polynominal g_0
|
||
# l0,1 l0,2 l0,3 l0,4 l0,5 l0,6
|
||
# 0.27505 0.35511 −0.004 −0.322 0.299 −0.0131
|
||
l_0 = np.array([0.27505, 0.35511, -0.004, -0.322, 0.299, -0.0131])
|
||
|
||
# Coefficients of polynominal g_1
|
||
# l1,1 l1,2 l1,3 l1,4 l1,5 l1,6 l1,7 l1,8 l1,9 l1,10 l1,11
|
||
# −10.036 −0.019804 0.55438 0.14108 12.494 0 0 0 0 0 1.0
|
||
l_1 = np.array([-10.036, -0.019804, 0.55438, 0.14108, 12.494, 0, 0, 0, 0, 0, 1.0])
|
||
|
||
# Coefficients of polynominal g_2
|
||
# l2,1 l2,2 l2,3 l2,4 l2,5 l2,6 l2,7 l2,8 l2,9 l2,10
|
||
# 0.42720 0.069884 −0.33771 0.24690 −1.0821 −0.30401 −1.1024 −1.2596 −0.31949 −1.4864
|
||
l_2 = np.array([0.42720, 0.069884, -0.33771, 0.24690, -1.0821, -0.30401, -1.1024, -1.2596, -0.31949, -1.4864])
|
||
|
||
l = [l_0, l_1, l_2]
|
||
|
||
# Function to calculate ratio of polynominals
|
||
def g_n(n, x, y):
|
||
if n == 0:
|
||
return l_0[0] * y - l_0[1] * x + l_0[2] / (l_0[3] - x)**2 + (l_0[4] - y)**2 + l_0[5]
|
||
if n == 1:
|
||
return l_1[0]
|
||
|
||
l[n][0] * (x + l[n][1])**2 + l[n][2] * (y + l[n][3])**2 + l[n][4] * x * y / l[n][5] * (x + l[n][6])**2 + l |