Compare commits
4 commits
| Author | SHA1 | Date | |
|---|---|---|---|
 Inex Code
|
f2bc29dad0 | ||
|
Inex Code
|
7d7d68ad67 | ||
|
Inex Code
|
b92c7be4dc | ||
|
Inex Code
|
3770c4dc2f |
45
etc/circle.py
Normal file
45
etc/circle.py
Normal 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
12
etc/decorators.py
Normal 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
33
etc/kindahtml.py
Normal 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
11
etc/nasa.py
Normal 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())
|
||||||
214
sat.py
214
sat.py
|
|
@ -1,6 +1,9 @@
|
||||||
from PIL import Image
|
from PIL import Image
|
||||||
import requests
|
import requests
|
||||||
from sat7_pointer import *
|
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
|
# Load Landsat 7 band 1, 2, 3 TIF images and create a composite image
|
||||||
# from the three bands.
|
# from the three bands.
|
||||||
|
|
@ -13,17 +16,17 @@ from sat7_pointer import *
|
||||||
#
|
#
|
||||||
# The composite image is saved as a PNG file.
|
# 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.
|
# Load corner coordinates of the image.
|
||||||
#
|
#
|
||||||
# The coordinates are stored in a MTL text file.
|
# 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.
|
# 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.
|
# # 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)
|
response = requests.get(url)
|
||||||
data = response.json()
|
data = response.json()
|
||||||
|
|
||||||
# Convert bounding box coordinates to image coordinates.
|
# 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(x0, y0)
|
||||||
print(x1, y1)
|
print(x1, y1)
|
||||||
|
|
@ -50,7 +57,7 @@ print(composite.size)
|
||||||
|
|
||||||
cropped = composite.crop((x0, y0, x1, y1))
|
cropped = composite.crop((x0, y0, x1, y1))
|
||||||
|
|
||||||
cropped.save('cropped.png')
|
cropped.save("cropped.png")
|
||||||
|
|
||||||
###############################################################################
|
###############################################################################
|
||||||
# LAB 2
|
# LAB 2
|
||||||
|
|
@ -59,7 +66,7 @@ cropped.save('cropped.png')
|
||||||
# Load landsat band 4 TIF image.
|
# Load landsat band 4 TIF image.
|
||||||
# Band 4 is near infrared.
|
# 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.
|
# Claculate the NDVI.
|
||||||
#
|
#
|
||||||
|
|
@ -115,7 +122,8 @@ def get_color(value):
|
||||||
else:
|
else:
|
||||||
return (0, 0, 0)
|
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]):
|
for y in range(ndvi.size[1]):
|
||||||
r = red.getpixel((x, y))
|
r = red.getpixel((x, y))
|
||||||
nir = ndvi.getpixel((x, y))
|
nir = ndvi.getpixel((x, y))
|
||||||
|
|
@ -124,6 +132,188 @@ for x in range(ndvi.size[0]):
|
||||||
else:
|
else:
|
||||||
result.putpixel((x, y), get_color((nir - r) / (nir + r)))
|
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))
|
||||||
|
|
@ -51,6 +51,20 @@ def load_metadata(mtl_file):
|
||||||
|
|
||||||
return mtl_data
|
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):
|
def lat_lot_to_pixel(lat, lon, mtl_data):
|
||||||
target_vector = np.array([lat, lon, 0])
|
target_vector = np.array([lat, lon, 0])
|
||||||
|
|
||||||
|
|
|
||||||
Loading…
Reference in a new issue