Shuichi/isak
1
0
Fork 0
forked from inex/isak
Code Issues Pull requests Projects Releases Wiki Activity

Compare commits

merge into: Shuichi:master
Branches Tags
Shuichi:master
inex:master
...
pull from: inex:master
Branches Tags
inex:master
Shuichi:master

4 commits
master ... master

Author SHA1 Message Date
Inex Code f2bc29dad0 Add etc tasks 2021-12-22 11:16:42 +03:00
Inex Code 7d7d68ad67 FAPAR to KML! 2021-12-17 22:32:04 +03:00
Inex Code b92c7be4dc Add FAPAR 2021-12-17 21:30:42 +03:00
Inex Code 3770c4dc2f some polynoms 2021-12-17 18:55:29 +03:00
8 changed files with 14721 additions and 12 deletions
Show stats Expand all files Collapse all files

45
etc/circle.py Normal file
View file

@ -0,0 +1,45 @@
class Shape:
def __init__(self, x, y):
self.x = x
self.y = y
def move(self, delta_x, delta_y):
self.x = self.x + delta_x
self.y = self.y + delta_y
class Square(Shape):
def __init__(self, side=1, x=0, y=0):
super().__init__(x, y)
self.side = side
class Circle(Shape):
pi = 3.14159
all_circles = []
def __init__(self, radius=1, x=0, y=0):
super().__init__(x, y)
self.radius = radius
self.all_circles.append(self)
@property
def radius(self):
return self._radius
@radius.setter
def radius(self, value):
if value < 0:
raise ValueError("Radius cannot be negative")
self._radius = value
@classmethod
def total_area(cls):
area = 0
for circle in cls.all_circles:
area += cls.circle_area(circle.radius)
return area
@staticmethod
def circle_area(radius):
return __class__.pi * radius * radius

12
etc/decorators.py Normal file
View file

@ -0,0 +1,12 @@
def wrap_with_html(fn):
def wrapper_func(*args):
print('<html>')
fn(*args)
print('</html>')
return wrapper_func
@wrap_with_html
def echo(param):
print(param)
echo('Hello')

33
etc/kindahtml.py Normal file
View file

@ -0,0 +1,33 @@
class HtmlElement:
tag = "html"
indent = " "
def __init__(self, content=None, **kwargs):
if content is None:
self.contents = []
else:
self.contents = [content]
self.attributes = kwargs
def append(self, new_content):
self.contents.append(new_content)
def render(self, cur_ind=""):
print(cur_ind + f"<{self.tag}>")
for content in self.contents:
try:
content.render(cur_ind + self.indent)
except AttributeError:
print(cur_ind + self.indent + content)
print(cur_ind + f"</{self.tag}>")
class Body(HtmlElement):
tag = "body"
class Paragraph(HtmlElement):
tag = "p"
para = Paragraph("hello")
body = Body(para)
body.append("world")
doc = HtmlElement(body)
doc.render()

11
etc/nasa.py Normal file
View file

@ -0,0 +1,11 @@
import requests
import json
api_key = "DEMO_KEY"
def get_weather():
url = f"https://api.nasa.gov/insight_weather/?api_key={api_key}&feedtype=json&ver=1.0"
response = requests.get(url)
return json.loads(response.text)
print(get_weather())

14404
fapar.kml Normal file
View file

File diff suppressed because it is too large Load diff

BIN
fapar.png Normal file
View file

Binary file not shown.
Side by side

After

Width:  |  Height:  |  Size: 689 KiB

214
sat.py
View file

@ -1,6 +1,9 @@
from PIL import Image
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.
@ -13,17 +16,17 @@ from sat7_pointer import *
#
# 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')
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))
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.
#
@ -34,14 +37,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)
x1, y0 = lat_lot_to_pixel(data[0]['boundingbox'][1], data[0]['boundingbox'][3], mtl_data)
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)
@ -50,7 +57,7 @@ print(composite.size)
cropped = composite.crop((x0, y0, x1, y1))
cropped.save('cropped.png')
cropped.save("cropped.png")
###############################################################################
# LAB 2
@ -59,7 +66,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.
#
@ -115,7 +122,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))
@ -124,6 +132,188 @@ 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)
# for each pixel in the cropped image.
#
# Bands 1, 3, 4 are used.
solar_zenith_angle = np.radians(float(mtl_data["SUN_ELEVATION"]))
sensor_zenith_angle = np.radians(0)
sun_sensor_relative_azimuth = np.radians(float(mtl_data["SUN_AZIMUTH"]))
gain = [float(mtl_data["RADIANCE_MULT_BAND_" + str(i)]) for i in [1, 3, 4]]
offset = [float(mtl_data["RADIANCE_ADD_BAND_" + str(i)]) for i in [1, 3, 4]]
dsol = float(mtl_data["EARTH_SUN_DISTANCE"])
pic = [0.643, 0.80760, 0.89472]
k = [0.76611, 0.63931, 0.81037]
theta = [-0.10055, -0.06156, -0.03924]
k = [0.63931,0.81037, 0.76611]
pic = [0.80760, 0.89472, 0.643]
theta = [-0.06156, -0.03924, -0.10055]
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')
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))

14
sat7_pointer/__init__.py
View file

@ -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])