Image processing using thermal data¶

Notebook prepared by Ben Maathuis, Willem Nieuwenhuis and Bas Retsios. ITC-University of Twente, Enschede. The Netherlands

Within this notebook you are going to retrieve from thermal image data:

  • the land surface temperature using emissivity derived from the ASTER emissivity data base (The ASTER Global Emissivity Database - https://emissivity.jpl.nasa.gov/aster-ged, contained wihin the Lamdsat product) and calculated yourself using a NDVI approach
  • apply a thermal image to derive the evapotranspiration using the DATTUTDUT model

Thermal Processing Landsat - Urban Heat Island Amsterdam¶

In [1]:
# import libraries
import os
import sys
import ilwis
import numpy as np
#import numpy.ma as ma
import netCDF4 as nc
from math import *
import matplotlib.pyplot as plt
import warnings

Retrieve thermal image data: Copernicus Global Land Collection¶

In openEO, Earth observation data sources are organised in collections. To access Landsat Thermal image data, using the Copernicus Global Land Data Service at VITO (url="openeo.vito.be"), we have used the collection_id = "LANDSAT8-9_L2". Landsat 8-9 Level 2 collection includes both Landsat 8 and the most recently launched Landsat 9 satellites (provided by NASA/USGS), both carrying the Operational Land Imager (OLI) and the Thermal Infrared Sensor (TIRS) instruments, with 9 optical and 2 thermal bands. These two sensors provide seasonal coverage of the global landmass. Landsat 8-9 Level 2 data from the most recently released collection 2, provides atmospherically corrected Surface Reflectance and Surface Brightness Temperature products generated from Collection 2 Level-1 scenes that have been processed to Tier 1 or Tier 2. Collection 2 level 2 data are available since February 2013 for Landsat 8 and since January 2022 for Landsat 9 and new data are continously added usually within 24 hours after Level 1 scenes are available.

A small data slice, given a selected temporal and spatial range, of the collection was downloaded. Details are:

  • temporal_extent = ["2023-07-07", "2023-07-08"] # Sunny day
  • Specify Region of Interest in Longitude-Latitude coordinates (degree.decimals)

Here region of interest specified over Amsterdam:

  • spatial_extent = {"west": 4.84, "south": 52.28, "east": 5.01, "north": 52.44}

The dataset was downloaded and is available in your sample data folder, file format is NetCDF

In [2]:
work_dir = os.getcwd()+r'/T_ET'

print("current dir is: %s" % (os.getcwd()))
print("current working directory is:",work_dir) 

if os.path.isdir(work_dir):
    print("Folder exists")
else:
    print("Folder doesn't exists")
    os.makedirs(work_dir, exist_ok=True)

current dir is: d:\jupyter\notebook_scripts\ilwispy_tutorial\tutorial_release
current working directory is: d:\jupyter\notebook_scripts\ilwispy_tutorial\tutorial_release/T_ET
Folder exists
In [3]:
#set the working directory for ILWISPy
ilwis.setWorkingCatalog(work_dir)
print(work_dir)

d:\jupyter\notebook_scripts\ilwispy_tutorial\tutorial_release/T_ET

Load the dataset and review the metadata¶

In [4]:
#the sample data is expected within a specific T_ET folder within the notebook folder  
fn = work_dir+'/L8_9_AMS.nc'
ds = nc.Dataset(fn)
In [5]:
print(ds)

<class 'netCDF4.Dataset'>
root group (NETCDF4_CLASSIC data model, file format HDF5):
    Conventions: CF-1.9
    institution: openEO platform
    dimensions(sizes): t(1), y(605), x(403)
    variables(dimensions): int32 t(t), float64 x(x), float64 y(y), |S1 crs(), float32 B01(t, y, x), float32 B02(t, y, x), float32 B03(t, y, x), float32 B04(t, y, x), float32 B05(t, y, x), float32 B06(t, y, x), float32 B07(t, y, x), float32 B10(t, y, x), float32 BQA(t, y, x), float32 QA_RADSAT(t, y, x), float32 SR_QA_AEROSOL(t, y, x), float32 ST_QA(t, y, x), float32 ST_TRAD(t, y, x), float32 ST_URAD(t, y, x), float32 ST_DRAD(t, y, x), float32 ST_ATRAN(t, y, x), float32 ST_EMIS(t, y, x), float32 ST_EMSD(t, y, x), float32 ST_CDIST(t, y, x), float32 dataMask(t, y, x)
    groups:
In [6]:
for dim in ds.dimensions.values():
    print(dim)

"<class 'netCDF4.Dimension'>" (unlimited): name = 't', size = 1
"<class 'netCDF4.Dimension'>": name = 'y', size = 605
"<class 'netCDF4.Dimension'>": name = 'x', size = 403
In [7]:
for var in ds.variables.values():
    print(var)

<class 'netCDF4.Variable'>
int32 t(t)
    standard_name: t
    long_name: t
    units: days since 1990-01-01
    axis: T
unlimited dimensions: t
current shape = (1,)
filling on, default _FillValue of -2147483647 used
<class 'netCDF4.Variable'>
float64 x(x)
    standard_name: projection_x_coordinate
    long_name: x coordinate of projection
    units: m
unlimited dimensions: 
current shape = (403,)
filling on, default _FillValue of 9.969209968386869e+36 used
<class 'netCDF4.Variable'>
float64 y(y)
    standard_name: projection_y_coordinate
    long_name: y coordinate of projection
    units: m
unlimited dimensions: 
current shape = (605,)
filling on, default _FillValue of 9.969209968386869e+36 used
<class 'netCDF4.Variable'>
|S1 crs()
    crs_wkt: PROJCS["WGS 84 / UTM zone 31N", GEOGCS["WGS 84", DATUM["World Geodetic System 1984", SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], UNIT["degree", 0.017453292519943295], AXIS["Geodetic longitude", EAST], AXIS["Geodetic latitude", NORTH], AUTHORITY["EPSG","4326"]], PROJECTION["Transverse_Mercator", AUTHORITY["EPSG","9807"]], PARAMETER["central_meridian", 3.0], PARAMETER["latitude_of_origin", 0.0], PARAMETER["scale_factor", 0.9996], PARAMETER["false_easting", 500000.0], PARAMETER["false_northing", 0.0], UNIT["m", 1.0], AXIS["Easting", EAST], AXIS["Northing", NORTH], AUTHORITY["EPSG","32631"]]
    spatial_ref: PROJCS["WGS 84 / UTM zone 31N", GEOGCS["WGS 84", DATUM["World Geodetic System 1984", SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]], PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]], UNIT["degree", 0.017453292519943295], AXIS["Geodetic longitude", EAST], AXIS["Geodetic latitude", NORTH], AUTHORITY["EPSG","4326"]], PROJECTION["Transverse_Mercator", AUTHORITY["EPSG","9807"]], PARAMETER["central_meridian", 3.0], PARAMETER["latitude_of_origin", 0.0], PARAMETER["scale_factor", 0.9996], PARAMETER["false_easting", 500000.0], PARAMETER["false_northing", 0.0], UNIT["m", 1.0], AXIS["Easting", EAST], AXIS["Northing", NORTH], AUTHORITY["EPSG","32631"]]
unlimited dimensions: 
current shape = ()
filling on, default _FillValue of  used
<class 'netCDF4.Variable'>
float32 B01(t, y, x)
    long_name: B01
    units: 
    _FillValue: 0.0
    grid_mapping: crs
unlimited dimensions: t
current shape = (1, 605, 403)
filling on
<class 'netCDF4.Variable'>
float32 B02(t, y, x)
    long_name: B02
    units: 
    _FillValue: 0.0
    grid_mapping: crs
unlimited dimensions: t
current shape = (1, 605, 403)
filling on
<class 'netCDF4.Variable'>
float32 B03(t, y, x)
    long_name: B03
    units: 
    _FillValue: 0.0
    grid_mapping: crs
unlimited dimensions: t
current shape = (1, 605, 403)
filling on
<class 'netCDF4.Variable'>
float32 B04(t, y, x)
    long_name: B04
    units: 
    _FillValue: 0.0
    grid_mapping: crs
unlimited dimensions: t
current shape = (1, 605, 403)
filling on
<class 'netCDF4.Variable'>
float32 B05(t, y, x)
    long_name: B05
    units: 
    _FillValue: 0.0
    grid_mapping: crs
unlimited dimensions: t
current shape = (1, 605, 403)
filling on
<class 'netCDF4.Variable'>
float32 B06(t, y, x)
    long_name: B06
    units: 
    _FillValue: 0.0
    grid_mapping: crs
unlimited dimensions: t
current shape = (1, 605, 403)
filling on
<class 'netCDF4.Variable'>
float32 B07(t, y, x)
    long_name: B07
    units: 
    _FillValue: 0.0
    grid_mapping: crs
unlimited dimensions: t
current shape = (1, 605, 403)
filling on
<class 'netCDF4.Variable'>
float32 B10(t, y, x)
    long_name: B10
    units: 
    _FillValue: 0.0
    grid_mapping: crs
unlimited dimensions: t
current shape = (1, 605, 403)
filling on
<class 'netCDF4.Variable'>
float32 BQA(t, y, x)
    long_name: BQA
    units: 
    _FillValue: 0.0
    grid_mapping: crs
unlimited dimensions: t
current shape = (1, 605, 403)
filling on
<class 'netCDF4.Variable'>
float32 QA_RADSAT(t, y, x)
    long_name: QA_RADSAT
    units: 
    _FillValue: 0.0
    grid_mapping: crs
unlimited dimensions: t
current shape = (1, 605, 403)
filling on
<class 'netCDF4.Variable'>
float32 SR_QA_AEROSOL(t, y, x)
    long_name: SR_QA_AEROSOL
    units: 
    _FillValue: 0.0
    grid_mapping: crs
unlimited dimensions: t
current shape = (1, 605, 403)
filling on
<class 'netCDF4.Variable'>
float32 ST_QA(t, y, x)
    long_name: ST_QA
    units: 
    _FillValue: 0.0
    grid_mapping: crs
unlimited dimensions: t
current shape = (1, 605, 403)
filling on
<class 'netCDF4.Variable'>
float32 ST_TRAD(t, y, x)
    long_name: ST_TRAD
    units: 
    _FillValue: 0.0
    grid_mapping: crs
unlimited dimensions: t
current shape = (1, 605, 403)
filling on
<class 'netCDF4.Variable'>
float32 ST_URAD(t, y, x)
    long_name: ST_URAD
    units: 
    _FillValue: 0.0
    grid_mapping: crs
unlimited dimensions: t
current shape = (1, 605, 403)
filling on
<class 'netCDF4.Variable'>
float32 ST_DRAD(t, y, x)
    long_name: ST_DRAD
    units: 
    _FillValue: 0.0
    grid_mapping: crs
unlimited dimensions: t
current shape = (1, 605, 403)
filling on
<class 'netCDF4.Variable'>
float32 ST_ATRAN(t, y, x)
    long_name: ST_ATRAN
    units: 
    _FillValue: 0.0
    grid_mapping: crs
unlimited dimensions: t
current shape = (1, 605, 403)
filling on
<class 'netCDF4.Variable'>
float32 ST_EMIS(t, y, x)
    long_name: ST_EMIS
    units: 
    _FillValue: 0.0
    grid_mapping: crs
unlimited dimensions: t
current shape = (1, 605, 403)
filling on
<class 'netCDF4.Variable'>
float32 ST_EMSD(t, y, x)
    long_name: ST_EMSD
    units: 
    _FillValue: 0.0
    grid_mapping: crs
unlimited dimensions: t
current shape = (1, 605, 403)
filling on
<class 'netCDF4.Variable'>
float32 ST_CDIST(t, y, x)
    long_name: ST_CDIST
    units: 
    _FillValue: 0.0
    grid_mapping: crs
unlimited dimensions: t
current shape = (1, 605, 403)
filling on
<class 'netCDF4.Variable'>
float32 dataMask(t, y, x)
    long_name: dataMask
    units: 
    _FillValue: 0.0
    grid_mapping: crs
unlimited dimensions: t
current shape = (1, 605, 403)
filling on

Import Thermal Channel (band 10)¶

In [8]:
data = ds['B10'][0,:,:]
In [9]:
print(type(data))
print(data.shape)

<class 'numpy.ma.MaskedArray'>
(605, 403)
In [10]:
plt.figure(figsize=(5,7))
plt.imshow(data, cmap = 'jet', interpolation = 'None') 
plt.title('Temperature_Kelvin')
plt.colorbar(shrink=0.4)
#plt.grid(False)
plt.show()
No description has been provided for this image

Data processing using ILWISPy¶

Now continue with data processing using ILWISpy

In [11]:
#define projection
ysize,xsize=np.shape(data)
grf = ilwis.GeoReference('code=georef:type=corners,csy=proj4:+proj=utm +proj=utm +zone=31 +north +ellps=WGS84 +datum=WGS84 +units=m, envelope=625050 5811900 637140 5793750, gridsize=' + str(xsize) + ' ' + str(ysize) + ',cornerofcorners=yes')
In [12]:
#create empty ilwis raster 
rc = ilwis.RasterCoverage()
defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(200, 350, 0))
rc.setDataDef(defNumr)
rc.setSize(ilwis.Size(xsize, ysize))
rc.setGeoReference(grf)
In [13]:
#assign the values to the newly created raster
rc.array2raster(np.array(data).flatten())
print(rc.size())

Size(403, 605, 1)
In [14]:
# Temperature from Kelvin to Celsius
Bt_C = ilwis.do('mapcalc', '(@1-273.15)', rc) # transform to degree Celsius
In [15]:
Bt_C.store('Bt_C.mpr') # check the image created with iliws386

Calculate LST using Emissivity information contained in "ST_EMIS"¶

In [16]:
#See: https://emissivity.jpl.nasa.gov/aster-ged
# read the emissivity values from the layer 'ST_EMIS'
data = ds['ST_EMIS'][0,:,:]
In [17]:
plt.figure(figsize=(5,7))
plt.imshow(data, cmap = 'jet', interpolation = 'None')
plt.title('Aster derived Emissivity')
plt.colorbar(shrink=0.4)
#plt.grid(False)
plt.show()
No description has been provided for this image
In [18]:
#define projection
ysize,xsize=np.shape(data)
grf = ilwis.GeoReference('code=georef:type=corners,csy=proj4:+proj=utm +proj=utm +zone=31 +north +ellps=WGS84 +datum=WGS84 +units=m, envelope=625050 5811900 637140 5793750, gridsize=' + str(xsize) + ' ' + str(ysize) + ',cornerofcorners=yes')
In [19]:
print(np.shape(data))

(605, 403)
In [20]:
#create empty ilwis raster 
rc = ilwis.RasterCoverage()
defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 1, 0))
rc.setDataDef(defNumr)
rc.setSize(ilwis.Size(xsize, ysize))
rc.setGeoReference(grf)
In [21]:
rc.array2raster(np.array(data).flatten())
print(rc.size())

Size(403, 605, 1)
In [22]:
rc.store('st_emis.mpr')

The formula below (brightness temperature /(1+(0.00115brightness temperature/1.4388)log10(emissivity)), is used to correct for emissivity in brightness temperature measurements, particularly in remote sensing applications. It adjusts the measured brightness temperature (which is a measure of radiance expressed as a temperature) to account for the fact that real-world surfaces don't perfectly emit radiation like a black body. The formula essentially estimates the true physical temperature of the surface based on its brightness temperature and emissivity.

Here's a breakdown of the formula:

  • Brightness Temperature (Tb): This is the measured temperature of the object, expressed as the temperature of a black body that would emit the same intensity of radiation.
  • Emissivity (ε): This is a measure of how well a surface emits thermal radiation compared to a black body. It ranges from 0 to 1, with 1 being a perfect black body.
    1. 00115: This is a constant, related to the physical properties of the atmosphere and the wavelength of microwave radiation.
    1. 4388: This is another constant, derived from the Wien's displacement law and represents a scaling factor related to the wavelength of the radiation.
  • log10(ε): The logarithm of emissivity is used to linearize the relationship between emissivity and brightness temperature. The formula essentially calculates a correction factor (the denominator) that is applied to the brightness temperature to obtain a more accurate estimate of the true physical temperature. This correction is important because surfaces with lower emissivities will have lower brightness temperatures than their actual physical temperature. In essence, the formula is used to estimate the surface temperature by taking into account the emissivity of the surface, which varies depending on the material and its condition.

Additional reading: https://www.researchgate.net/publication/296414003_Algorithm_for_Automated_Mapping_of_Land_Surface_Temperature_Using_LANDSAT_8_Satellite_Data

In [23]:
#calculate LST
LST1 = ilwis.do('mapcalc', '(@1/(1+(0.00115*@1/1.4388)*log10(@2)))', Bt_C, rc)
In [24]:
LST1.store('LST1.mpr')

Calculate LST using emissivity derived from NDVI¶

The formula applied below ε = 0.004 * Pv + 0.986 is used to calculate land surface emissivity (ε) from the proportion of vegetation (Pv) derived from the Normalized Difference Vegetation Index (NDVI) (see: https://www.uv.es/ucg/articulos/2005/Publications_2004_10.pdf). This formula is commonly used in remote sensing applications, particularly when estimating land surface temperature from satellite imagery like Landsat. The 0.986 value is a standard emissivity value for vegetation, while the 0.004 represents the emissivity variation due to vegetation.

Here's a more detailed breakdown:

  • NDVI and Proportion of Vegetation (Pv): The Normalized Difference Vegetation Index (NDVI) is a widely used indicator of vegetation health and density. It's calculated from the red and near-infrared spectral bands of satellite imagery. The proportion of vegetation (Pv) is then derived from the NDVI values, representing the fraction of vegetation within a given pixel.

  • Emissivity (ε): Emissivity is a measure of how efficiently a surface emits thermal radiation. It's a crucial parameter for accurate land surface temperature (LST) calculations.

  • The Formula: The formula ε = 0.004 * Pv + 0.986 relates the proportion of vegetation (Pv) to the land surface emissivity (ε). The coefficient 0.004 is a factor that accounts for the emissivity variation due to vegetation. The value 0.986 is a standard emissivity value for vegetation. This means that as the proportion of vegetation increases (higher NDVI values), the emissivity of the surface will also increase.

In [25]:
# import bands 4 and 5 (red and nir)
b4data = ds['B04'][0,:,:]
b5data = ds['B05'][0,:,:]
In [26]:
#create empty ilwis raster 
b4 = ilwis.RasterCoverage()
defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 100000, 0))
b4.setDataDef(defNumr)
b4.setSize(ilwis.Size(xsize, ysize))
b4.setGeoReference(grf)
In [27]:
#create empty ilwis raster 
b5 = ilwis.RasterCoverage()
defNumr = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(0, 1000000, 0))
b5.setDataDef(defNumr)
b5.setSize(ilwis.Size(xsize, ysize))
b5.setGeoReference(grf)
In [28]:
b4.array2raster(np.array(b4data).flatten())
print(b4.size())

Size(403, 605, 1)
In [29]:
b5.array2raster(np.array(b5data).flatten())
print(b5.size())

Size(403, 605, 1)
In [30]:
b4.store('b4.mpr')
b5.store('b5.mpr')
In [31]:
#calculate ndvi
ndvi = ilwis.do('mapcalc', '(@2-@1)/(@2+@1)', b4, b5) 
print(ndvi.size())

Size(403, 605, 1)
In [32]:
ndvi.store('ndvi.mpr')
In [33]:
stats = ndvi.statistics(ilwis.PropertySets.pHISTOGRAM)
In [34]:
ndv_min= (stats[ilwis.PropertySets.pMIN]) # minimum value on the map
ndv_max = (stats[ilwis.PropertySets.pMAX]) # maximum value on the map
In [35]:
print(ndv_min)
print(ndv_max)

-1.0
1.0
In [36]:
pv = ilwis.do('mapcalc', 'sqrt(@1-@2)/(@3-@2)', ndvi, ndv_min, ndv_max)
In [37]:
pv.store('pv.mpr')
In [38]:
emis_ndv = ilwis.do('mapcalc', '(0.004*@1+0.986)', pv)
In [39]:
emis_ndv.store('emis_ndv.mpr')
In [40]:
#calculate LST
LST2 = ilwis.do('mapcalc', '(@1/(1+(0.00115*@1/1.4388)*log10(@2)))', Bt_C, emis_ndv)
In [41]:
LST2.store('LST2.mpr')

Display the results using matplotlib¶

In [42]:
LST1_2np = np.fromiter(iter(LST1), np.float64, LST1.size().linearSize()) 
LST1_2np = LST1_2np.reshape((LST1.size().ysize, LST1.size().xsize))
In [43]:
LST2_2np = np.fromiter(iter(LST2), np.float64, LST2.size().linearSize()) 
LST2_2np = LST2_2np.reshape((LST2.size().ysize, LST2.size().xsize))
In [44]:
#create a plot
%matplotlib inline 
fig1 = plt.figure(figsize=(12, 8))
plt.subplot(1, 2, 1)
plt.imshow(LST1_2np, cmap = 'jet', interpolation = 'None', vmin=25, vmax=45)
plt.title('LST-Method 1')
plt.colorbar(shrink=0.4)

plt.subplot(1, 2, 2)
plt.imshow(LST2_2np, cmap = 'jet', interpolation = 'None', vmin=25, vmax=45)
plt.title('LST_Method 2')
plt.colorbar(shrink=0.4);
No description has been provided for this image

Dattutdut model: ET-calculation¶

Retrieve ET from thermal images, Aster level 1b thermal channel is used here as input image from 18 July 2004, acquired at 11:00 UTC. Barrax (ESA study site) - Spain (ast_l1b_00307182004110047_08122004112525).

In the DATTUTDUT model, the minimum temperature (Tmin) is typically determined as the 0.5% quantile of the coldest pixels in a land surface temperature (LST) map. This means that the coldest 0.5% of the pixels are identified, and their temperature is taken as the minimum temperature value. This value, along with the maximum temperature (Tmax), is used to scale pixels within the scene, with hotter pixels representing areas with less evapotranspiration and colder pixels representing areas with higher evapotranspiration.

  • Land Surface Temperature (LST) Map: The DATTUTDUT model relies on an LST map as its primary input.
  • Quantile Calculation: The model identifies the coldest pixels in the LST map and calculates the 0.5% quantile.
  • Minimum Temperature (Tmin): The temperature value corresponding to the 0.5% quantile is designated as the minimum temperature (Tmin) for the scene.
  • Scaling Pixels: Tmin and Tmax (the hottest pixel temperature) are then used to scale the temperature values of all other pixels within the scene.
  • Evapotranspiration Estimation: This scaling process helps in estimating evapotranspiration, with hotter pixels indicating lower evapotranspiration and colder pixels indicating higher evapotranspiration.

Additional reading: Chapter2 - https://research.utwente.nl/files/279986011/timmermans.pdf

In [45]:
#Retrieve input ASTER L1b thermal image - brightness temperature in Kelvin
Input_file = work_dir+r'/trad200.mpr'
print(Input_file)

# open a raster coverage object:
raster = ilwis.RasterCoverage(Input_file)

# probe its size
print(raster.size().xsize)
print(raster.size().ysize)

d:\jupyter\notebook_scripts\ilwispy_tutorial\tutorial_release/T_ET/trad200.mpr
64
64
In [46]:
#retrieve some image statistics
pThermal = ilwis.PropertySets.pMIN|ilwis.PropertySets.pMAX|ilwis.PropertySets.pDISTANCE|ilwis.PropertySets.pDELTA|ilwis.PropertySets.pNETTOCOUNT|ilwis.PropertySets.pCOUNT|ilwis.PropertySets.pSUM|ilwis.PropertySets.pMEAN|ilwis.PropertySets.pMEDIAN|ilwis.PropertySets.pSTDEV|ilwis.PropertySets.pHISTOGRAM
Image_stat = raster.statistics(pThermal)
print("Minimum:", (Image_stat.prop(ilwis.PropertySets.pMIN)))
print("Maximum:", (Image_stat.prop(ilwis.PropertySets.pMAX)))
print("Mean:", (Image_stat.prop(ilwis.PropertySets.pMEAN)))
print("Median:", (Image_stat.prop(ilwis.PropertySets.pMEDIAN)))
print("STD:", (Image_stat.prop(ilwis.PropertySets.pSTDEV)))
print("Count:", (Image_stat.prop(ilwis.PropertySets.pCOUNT)))

Minimum: 295.76635900212904
Maximum: 321.76386778598567
Mean: 311.0757558384233
Median: 313.1600490182734
STD: 6.633416144483858
Count: 4096.0
In [47]:
offset = 0.5
tmin_offset = Image_stat.calcStretchRange(offset)[0] #offset and first value (left - lower range of histogram)
print('Tmin + offset = ', tmin_offset)
tmax = Image_stat.prop(ilwis.PropertySets.pMAX)
print('Tmax = ',tmax)

Tmin + offset =  297.06623444132185
Tmax =  321.76386778598567
In [48]:
# values based on image acquisition date / time and brightness temperature recorded
# note Tmin has offset of 0.5 % from minimum temperature

Tmin = tmin_offset 
Tmax = tmax 

DOY = 200
UTC = 11.008
In [49]:
# display the image
#read raster as np-array and transform - reshape from 1D to 2D
raster_im = np.fromiter(iter(raster), np.float64, raster.size().linearSize())
raster_im = raster_im.reshape((raster.size().ysize, raster.size().xsize ))

#plot the result
plt.figure(figsize = (6,6))
plt.imshow(raster_im, interpolation='none', cmap ='terrain')
plt.title('Thermal band Aster-L1b')
plt.axis('off')
plt.colorbar(shrink=0.6)
plt.show()
No description has been provided for this image
In [50]:
# inspect image georeference and coordinate system

grf = raster.geoReference()
print (raster.geoReference().name())
csy = grf.coordinateSystem()
print('coordinate system name =', csy.name())
print()
print(csy.toProj4())
print()
print(csy.toWKT())

barrax4les.grf
coordinate system name = A20040718.csy

+proj=utm +zone=30 +a=6378137.000000000000 +b=6356752.314245179296

PROJCS["A20040718.csy",GEOCS["A20040718.csy",DATUM[" WGS 84",[DWGS84],ELLIPSOID["WGS 84",6378137.000000000000,298.257223563000],PRIMEM["Greenwich",0, AUTHORITY["EPSG",8901"]]],PROJECTION["Transverse_Mercator"],PARAMETER[,Yes]PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],PARAMETER["scale",0.9996],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-3],UNIT["meter",1.0]]
In [51]:
#retreive from UTM projection the coordinates (per pixel) and transform these into latitude and longitude

csyLL = ilwis.CoordinateSystem('epsg:4326') #define latlon to be used as output

sizex = grf.size().xsize
sizey = grf.size().ysize

latitudes = np.empty((sizey, sizex))
longitudes = np.empty((sizey, sizex))

for x in range(sizex):
    for y in range(sizey):
        crdUtm = grf.pixel2Coord(ilwis.Pixel(x,y))#read input coordinates
        crdLL = csyLL.coord2coord(csy, crdUtm)#transform coordinates to LatLon
        lon = crdLL.x
        lat = crdLL.y
        latitudes[y,x] = lat
        longitudes[y,x] = lon
In [52]:
#create latitude and longitude maps for the region
lonDOY = ilwis.RasterCoverage()
defll = ilwis.DataDefinition(ilwis.NumericDomain('code=value'), ilwis.NumericRange(-180, 180, 0))
lonDOY.setDataDef(defll)
lonDOY.setSize(ilwis.Size(sizex, sizey))
lonDOY.setGeoReference(grf)
lonDOY.array2raster(longitudes.flatten())

latDOY = ilwis.RasterCoverage()
latDOY.setDataDef(defll)
latDOY.setSize(ilwis.Size(sizex, sizey))
latDOY.setGeoReference(grf)
latDOY.array2raster(latitudes.flatten())
In [53]:
#utm projection has been transformed into lat / lon  coordinate maps
lonDOY.store('lon'+str(DOY)+'.mpr')
latDOY.store('lat'+str(DOY)+'.mpr')
In [54]:
#da = day angle
da = 0.01721420632*(DOY-1)
print(da)

#et = equation of time
et = (0.000075+0.001868*cos(da)-0.032077*sin(da)-0.014615*cos(2*da)-0.04089*sin(2*da))*229.18
print(et)

om200ega = ilwis.do('mapcalc','(12-((11.008)+@1/15+@2/60))*15', lonDOY, et)

3.4256270576800003
-6.198818682649018
In [55]:
#declination
de = (0.006918-0.399912*cos(da)+0.070257*sin(da)-0.006758*cos(2*da)+0.000907*sin(2*da)-0.002697*cos(3*da)+0.00148*sin(3*da))*57.29577951
print(de)

co200zen = ilwis.do('mapcalc','sin(@1/57.29577951)*sin(@2/57.29577951)+cos(@1/57.29577951)*cos(@2/57.29577951)*cos(@3/57.29577951)', de, latDOY, om200ega)

21.00300380507951
In [56]:
#eccentricity correction factor (sun orbit around the earth ) 
eo = 1.00011+0.034221*cos(da)+0.00128*sin(da)+0.000719*cos(2*da)+0.000077*sin(2*da)
print(eo)

ir200toa = ilwis.do('mapcalc', '1367*@1*@2', eo, co200zen)

0.9675489495248685
In [57]:
# Calculation of the net radiation albedo map - instantaneous
r200n = ilwis.do('mapcalc','(1-(0.05+(@1-@2)/(@3-@2)*0.2))*0.7*@4+0.8*0.000000056697*pow(@2,4)-1.0*0.000000056697*pow(@1,4)', raster, Tmin, Tmax, ir200toa)
In [58]:
# Calculation of the soil heat flux map - instantaneous
g200s = ilwis.do('mapcalc','(0.05+(@1-@2)/(@3-@2)*0.4)*@4', raster, Tmin, Tmax, r200n)
In [59]:
# Calculation of turbulent fluxes
TRange = Tmax-Tmin
rg = ilwis.do('mapcalc','@1-@2',r200n, g200s)

h200 = ilwis.do('mapcalc','iff(@1<@2,0,iff(@1>@3,@4,@4*(@1-@2)/@5))', raster, Tmin, Tmax, rg, TRange)
In [60]:
le200 = ilwis.do('mapcalc','@1 - @2 - @3', r200n, g200s, h200)
In [61]:
# Daily values:
lambda200 = ilwis.do('mapcalc','@1/(@1+@2)', le200, h200)
In [62]:
# Calculation of daily totals (assuming n/N equal to 1):
om200egasr = ilwis.do('mapcalc','acos(-1*tan(@1/57.29577951)*tan(@2/57.29577951))', latDOY, de)
In [63]:
ir200toaday = ilwis.do ('mapcalc','24/3.141592*1367*0.0036*@1*cos(@2/57.29577951)*cos(@3/57.29577951)*(sin(@4)-@4*cos(@4))', eo, latDOY, de, om200egasr)
In [64]:
qn200day = ilwis.do ('mapcalc', '(0.25+0.5*1)*(1-1.1*(0.05+(@1-@2)/@3*0.2))*@4-110*(0.7)/11.5741', raster, Tmin, TRange, ir200toaday)
In [65]:
et200day = ilwis.do ('mapcalc', '@1*@2/2.47', lambda200, qn200day)

#read raster as np-array and transform - reshape from 1D to 2D
et200day_im = np.fromiter(iter(et200day), np.float64, et200day.size().linearSize())
et200day_im = et200day_im.reshape((et200day.size().ysize, et200day.size().xsize ))

#plot the result
plt.figure(figsize = (6,6))
plt.imshow(et200day_im, interpolation='none', cmap ='jet')
plt.title('ET (mm) according to DATTUTDUT')
plt.axis('off')
plt.colorbar(shrink=0.6)
plt.show()
No description has been provided for this image
In [66]:
et200day.store('et'+str(DOY)+'sum_mm.mpr')