FAPAR to KML!

sat.py

@@ -3,6 +3,7 @@ import requests
 from sat7_pointer import *
 import numpy as np
 from tqdm import tqdm
+from lxml import etree
 
 # Load Landsat 7 band 1, 2, 3 TIF images and create a composite image
 # from the three bands.
@@ -277,3 +278,42 @@ for x in tqdm(range(result.size[0]), desc="FAPAR"):
         result.putpixel((x, y), get_color_fapar(FAPAR, rho_i))
 
 result.save('fapar.png')
+
+
+root = etree.Element("kml")
+doc = etree.SubElement(root, "Document")
+
+start_x = 550
+start_y = 900
+width = 30
+height = 30
+result = result.crop((start_x, start_y, start_x+width, start_y+height))
+
+for x in tqdm(range(result.size[0]), desc="KML"):
+    for y in range(result.size[1]):
+        color = result.getpixel((x, y))
+        if color == (0, 0, 0):
+            continue
+        coord = pixel_to_lat_lot(x+x0+start_x, y+y0+start_y, mtl_data)
+        neighbour_r = pixel_to_lat_lot(x+1+x0+start_x, y+y0+start_y, mtl_data)
+        neighbour_d = pixel_to_lat_lot(x+x0+start_x, y+1+y0+start_y, mtl_data)
+        neighbour_rd = pixel_to_lat_lot(x+1+x0+start_x, y+1+y0+start_y, mtl_data)
+        placemark = etree.SubElement(doc, "Placemark")
+        etree.SubElement(placemark, "name").text = f"{x}_{y}"
+        etree.SubElement(placemark, "styleUrl").text = f"#style_{x}_{y}"
+        polygon = etree.SubElement(placemark, "Polygon")
+        outer = etree.SubElement(polygon, "outerBoundaryIs")
+        linearring = etree.SubElement(outer, "LinearRing")
+        coordinates = etree.SubElement(linearring, "coordinates")
+        coordinates.text = f"{coord[0]},{coord[1]},0 {neighbour_r[0]},{neighbour_r[1]},0 {neighbour_rd[0]},{neighbour_rd[1]},0 {neighbour_d[0]},{neighbour_d[1]},0 {coord[0]},{coord[1]},0"
+        color = result.getpixel((x, y))
+        color = (color[2], color[1], color[0])
+        color = '%02x%02x%02x' % color
+        color = "64" + color
+        style = etree.SubElement(doc, "Style")
+        style.set("id", f"style_{x}_{y}")
+        polystyle = etree.SubElement(style, "PolyStyle")
+        etree.SubElement(polystyle, "color").text = color
+
+with open("fapar.kml", "wb") as f:
+    f.write(etree.tostring(root, pretty_print=True))

sat7_pointer/__init__.py

@@ -51,6 +51,20 @@ def load_metadata(mtl_file):
 
     return mtl_data
 
+def pixel_to_lat_lot(x, y, mtl_data):
+    k = (x - 0.5) / mtl_data['width']
+    r = (y - 0.5) / mtl_data['height']
+
+    y = k*(mtl_data['vector_a'] + mtl_data['mu'] * (mtl_data['vector_a'] - mtl_data['vector_s'])) \
+        + r * (mtl_data['vector_c'] + mtl_data['nu'] * (mtl_data['vector_c'] - mtl_data['vector_s']))
+    B = np.array([mtl_data['vector_a'], mtl_data['vector_c'], mtl_data['vector_s']-y]).T
+    ags = np.linalg.solve(B, y.T)
+    al = ags[0]
+    gm = ags[1]
+    ecef = mtl_data['ecef_corners'][0, :] + al*mtl_data['vector_a'] + gm*mtl_data['vector_c']
+    lat, lon, alt = pyproj.transform(mtl_data['ecef'], mtl_data['lla'], ecef[0], ecef[1], ecef[2], radians=False)
+    return [lat, lon]
+
 def lat_lot_to_pixel(lat, lon, mtl_data):
     target_vector = np.array([lat, lon, 0])