Melakukan image processing dari data input landsat berupa raster file dengan tipe tiff image, melakukan manipulasi pixel dengan operasi matematika (telah tersedia), kemudian plot pixel yang baru menjadi tiff image yang baru.
Program sudah ada, tetapi hasilnya masih belum sesuai, hanya perlu mengkoreksi dan menambahkan beberapa baris program jika diperlukan.
Pada dasarnya program sangat simple, untuk yang terbiasa koding beberapa jam saya yakin bisa selesai
import gdal
import sys, os
import numpy as np
def readGeo(rast):
"""
Function readGeo reads raster and mask files and makes Numpy array with
the data restricted by the mask.
Inputs:
- rast - raster in GDAL readable format
- mask - raster mask in GDAL readable format. Data should be 8bit integer,
typicaly 0 (nodata) and 1 (data).
"""
# raster processing
ds = gdal.Open(rast, GA_Update)
gtransf = ds.GetGeoTransform()
prj = ds.GetProjection()
rast_in = gdal.Dataset.ReadAsArray(ds).astype(np.float32)
ds = None
return rast_in, gtransf, prj
def outRast(rast_out, gtransf, prj, outFile):
driver = gdal.GetDriverByName("Gtiff")
ds = driver.Create(outFile, rast_out.shape[1],rast_out.shape[0], 1, gdal.GDT_Float32)
ds.SetProjection(prj)
ds.SetGeoTransform(gtransf)
ds.GetRasterBand(1).WriteArray(rast_out)
# close the dataset to flush it to disk
ds = None
def sensorRadiance(band, gain, offset):
L_lambda = gain * band + offset
return L_lambda
def brightTemperature(L_lambda, K1, K2):
T_B = K2/(np.log(K1/L_lambda + 1))
return T_B
#-------------------------------------------------------------------------------
# Jimenez-Munoz & Sobrino metoda
def calcGamma(L_lambda, T_B):
c1 = 1.19104 * 10**(8)
c2 = 14387.7
lambd = 11.457
gamma = 1/(c2*L_lambda/T_B**2 * (lambd**4/c1 * L_lambda + 1/lambd))
return gamma
def calcDelta(gamma, L_lambda, T_B):
delta = -1 * gamma * L_lambda + T_B
return delta
def calcPsi(humidity):
psi1 = 0.14714 * humidity**2 - 0.15583 * humidity + 1.1234
psi2 = (-1.1836) * humidity**2 - 0.37607 * humidity - 0.52894
psi3 = (-0.04554) * humidity**2 + 1.8719 * humidity - 0.39071
return psi1, psi2, psi3
def surfTempJMS(L_lambda, emis, gamma, delta, psi1, psi2, psi3):
Ts = gamma * (1/emis * (psi1 * L_lambda + psi2) + psi3) + delta - 273.16
return Ts
#-------------------------------------------------------------------------------
# vrstvy
thermal_band, gtransf, prj = readGeo('B10.TIF')
emis, gtransf, prj = readGeo('B10.TIF')
# mezivypocty
humidity = 1
Gain = 0.055376
Offset = 1.18
K1 = 607.76
K2 = 1260.56
L_lambda = sensorRadiance(thermal_band, Gain, Offset)
T_B = brightTemperature(L_lambda, K1, K2)
gamma = calcGamma(L_lambda, T_B)
delta = calcDelta(gamma, L_lambda, T_B)
psi1, psi2, psi3 = calcPsi(humidity)
# teplota povrchu
Ts_JMS = surfTempJMS(L_lambda, emis, gamma, delta, psi1, psi2, psi3)
# vystupy
outRast(Ts_JMS, gtransf, prj, Surface_temperature)
Loading ...
