Add FAPAR

sat.py

@@ -1,7 +1,8 @@
 from PIL import Image
 import requests
 from sat7_pointer import *
-from numpy import np
+import numpy as np
+from tqdm import tqdm
 
 # Load Landsat 7 band 1, 2, 3 TIF images and create a composite image
 # from the three bands.
@@ -14,17 +15,17 @@ from numpy import np
 #
 # The composite image is saved as a PNG file.
 
-band1 = Image.open('LE07_L1TP_177025_20210723_20210818_02_T1_B1.TIF')
+band1 = Image.open("LE07_L1TP_177025_20210723_20210818_02_T1_B1.TIF")
-band2 = Image.open('LE07_L1TP_177025_20210723_20210818_02_T1_B2.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')
+band3 = Image.open("LE07_L1TP_177025_20210723_20210818_02_T1_B3.TIF")
 
-composite = Image.merge('RGB', (band3, band2, band1))
+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')
+mtl_data = load_metadata("LE07_L1TP_177025_20210723_20210818_02_T1_MTL.txt")
 
 # Fetch coordinates of a city.
 #
@@ -35,14 +36,18 @@ mtl_data = load_metadata('LE07_L1TP_177025_20210723_20210818_02_T1_MTL.txt')
 #
 # # City is Belgorod, Russia.
 
-url = 'http://nominatim.openstreetmap.org/search?q=Belgorod,+Russia&format=json'
+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)
+x0, y1 = lat_lot_to_pixel(
-x1, y0 = lat_lot_to_pixel(data[0]['boundingbox'][1], data[0]['boundingbox'][3], mtl_data)
+    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)
@@ -51,7 +56,7 @@ print(composite.size)
 
 cropped = composite.crop((x0, y0, x1, y1))
 
-cropped.save('cropped.png')
+cropped.save("cropped.png")
 
 ###############################################################################
 # LAB 2
@@ -60,7 +65,7 @@ cropped.save('cropped.png')
 # Load landsat band 4 TIF image.
 # Band 4 is near infrared.
 
-band4 = Image.open('LE07_L1TP_177025_20210723_20210818_02_T1_B4.TIF')
+band4 = Image.open("LE07_L1TP_177025_20210723_20210818_02_T1_B4.TIF")
 
 # Claculate the NDVI.
 #
@@ -116,7 +121,8 @@ def get_color(value):
     else:
         return (0, 0, 0)
 
-for x in range(ndvi.size[0]):
+
+for x in tqdm(range(ndvi.size[0]), desc="NDVI"):
     for y in range(ndvi.size[1]):
         r = red.getpixel((x, y))
         nir = ndvi.getpixel((x, y))
@@ -125,7 +131,7 @@ for x in range(ndvi.size[0]):
         else:
             result.putpixel((x, y), get_color((nir - r) / (nir + r)))
 
-result.save('ndvi.png')
+result.save("ndvi.png")
 
 
 # Calculate FAPAR (Fraction of Absorbed Photosynthetically Active Radiation)
@@ -133,28 +139,141 @@ result.save('ndvi.png')
 #
 # Bands 1, 3, 4 are used.
 
-# Coefficients of polynominal g_0
+solar_zenith_angle = np.radians(float(mtl_data["SUN_ELEVATION"]))
-# l0,1 l0,2 l0,3 l0,4 l0,5 l0,6
+sensor_zenith_angle = np.radians(0)
-# 0.27505 0.35511 −0.004 −0.322 0.299 −0.0131
+sun_sensor_relative_azimuth = np.radians(float(mtl_data["SUN_AZIMUTH"]))
-l_0 = np.array([0.27505, 0.35511, -0.004, -0.322, 0.299, -0.0131])
 
-# Coefficients of polynominal g_1
+gain = [float(mtl_data["RADIANCE_MULT_BAND_" + str(i)]) for i in [1, 3, 4]]
-# l1,1 l1,2 l1,3 l1,4 l1,5 l1,6 l1,7 l1,8 l1,9 l1,10 l1,11
+offset = [float(mtl_data["RADIANCE_ADD_BAND_" + str(i)]) for i in [1, 3, 4]]
-# −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
+dsol = float(mtl_data["EARTH_SUN_DISTANCE"])
-# 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]
+pic = [0.643, 0.80760, 0.89472]
+k = [0.76611, 0.63931, 0.81037]
+theta = [-0.10055, -0.06156, -0.03924]
 
-# Function to calculate ratio of polynominals 
+k = [0.63931,0.81037, 0.76611]
-def g_n(n, x, y):
+pic = [0.80760, 0.89472, 0.643]
-    if n == 0:
+theta = [-0.06156, -0.03924, -0.10055]
-        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
+E0 = [1969, 1551, 1044]
+
+cosg = np.cos(solar_zenith_angle) * np.cos(sensor_zenith_angle) + np.sin(
+    solar_zenith_angle
+) * np.sin(sensor_zenith_angle) * np.cos(sun_sensor_relative_azimuth)
+G = (
+    np.tan(solar_zenith_angle) ** 2
+    + np.tan(sensor_zenith_angle) ** 2
+    - 2
+    * np.tan(solar_zenith_angle)
+    * np.tan(sensor_zenith_angle)
+    * np.cos(sun_sensor_relative_azimuth)
+) ** 0.5
+
+
+polynoms = np.array(
+    [
+        [0.27505, 0.35511, -0.004, -0.322, 0.299, -0.0131, 0, 0, 0, 0, 0],
+        [-10.036, -0.019804, 0.55438, 0.14108, 12.494, 0, 0, 0, 0, 0, 1],
+        [
+            0.42720,
+            0.069884,
+            -0.33771,
+            0.24690,
+            -1.0821,
+            -0.30401,
+            -1.1024,
+            -1.2596,
+            -0.31949,
+            -1.4864,
+            0,
+        ],
+    ]
+)
+
+blue = band1.copy().crop((x0, y0, x1, y1))
+red = band3.copy().crop((x0, y0, x1, y1))
+nir = band4.copy().crop((x0, y0, x1, y1))
+result = cropped.copy()
+
+f1 = [
+    ((np.cos(solar_zenith_angle) * np.cos(sensor_zenith_angle)) ** (k[i] - 1))
+    / (np.cos(solar_zenith_angle) + np.cos(sensor_zenith_angle)) ** (1 - k[i])
+    for i in range(3)
+]
+f2 = [
+    (1 - theta[i] ** 2) / (1 + 2 * theta[i] * cosg + theta[i] ** 2) ** (3 / 2)
+    for i in range(3)
+]
+f3 = [1 + (1 - pic[i]) / (1 + G) for i in range(3)]
+F = [f1[i] * f2[i] * f3[i] for i in range(3)]
+
+def get_color_fapar(value, rho):
+    if (0 < rho[0] and rho[0] < 0.257752) \
+        and (0 < rho[1] and rho[1] < 0.48407) \
+        and (0 < rho[2] and rho[2] < 0.683928) \
+        and (rho[0] <= rho[2]) \
+        and (rho[2] >= 1.26826*rho[1]):
+            return get_color(value)
+    if (rho[0] <= 0) or (rho[1] <= 0) or (rho[2] <= 0):
+        return (0, 0, 0)
+    if (rho[0] >= 0.257752) or (rho[1] >= 0.48407) or (rho[2] >= 0.683928):
+        return (255, 255, 255)
+    if (0 < rho[0] and rho[0] < 0.257752) \
+        and (0 < rho[1] and rho[1] < 0.48407) \
+        and (0 < rho[2] and rho[2] < 0.683928) \
+        and (rho[0] >= rho[2]):
+            return (0, 0, 255)
+    if (0 < rho[0] and rho[0] < 0.257752) \
+        and (0 < rho[1] and rho[1] < 0.48407) \
+        and (0 < rho[2] and rho[2] < 0.683928) \
+        and (rho[0] <= rho[2]) \
+        and (1.25*rho[1] > rho[2]):
+            return (255, 150, 150)
+    if (rho[1] < 0) or (rho[2] < 0):
+        return (0, 0, 0)
+    if value < 0 or value > 1:
+        return (0, 0, 0)
+    
+    return (int(180.0 * (1 - value)), 255, 255)
+
+for x in tqdm(range(result.size[0]), desc="FAPAR"):
+    for y in range(result.size[1]):
+        bands = [blue.getpixel((x, y)), red.getpixel((x, y)), nir.getpixel((x, y))]
+
+        rho_i = [
+            (
+                (np.pi * (gain[i] * bands[i] + offset[i]) * dsol ** 2)
+                / (E0[i] * np.cos(sensor_zenith_angle))
+            )
+            / F[i]
+            for i in range(3)
+        ]
+
+        g1 = (
+            (polynoms[1, 0] * (rho_i[0] + polynoms[1, 1]) ** 2)
+            + (polynoms[1, 2] * (rho_i[1] + polynoms[1, 3]) ** 2)
+            + polynoms[1, 4] * rho_i[0] * rho_i[1]
+        ) / (
+            polynoms[1, 5] * (rho_i[0] + polynoms[1, 6]) ** 2
+            + polynoms[1, 7] * (rho_i[1] + polynoms[1, 8]) ** 2
+            + polynoms[1, 9] * rho_i[0] * rho_i[1]
+            + polynoms[1, 10]
+        )
+        g2 = (
+            (polynoms[2, 0] * (rho_i[0] + polynoms[2, 1]) ** 2)
+            + (polynoms[2, 2] * (rho_i[2] + polynoms[2, 3]) ** 2)
+            + polynoms[2, 4] * rho_i[0] * rho_i[2]
+        ) / (
+            polynoms[2, 5] * (rho_i[0] + polynoms[2, 6]) ** 2
+            + polynoms[2, 7] * (rho_i[2] + polynoms[2, 8]) ** 2
+            + polynoms[2, 9] * rho_i[0] * rho_i[2]
+            + polynoms[2, 10]
+        )
+        FAPAR = ((polynoms[0, 0] * g2) - polynoms[0, 1] * g1 - polynoms[0, 2]) / (
+            (polynoms[0, 3] - g1) ** 2 + (polynoms[0, 4] - g2) ** 2 + polynoms[0, 5]
+        )
+
+        result.putpixel((x, y), get_color_fapar(FAPAR, rho_i))
+
+result.save('fapar.png')