Μελέτη Θερμικής Αστικής Νησίδας: Υπολογισμός θερμοκρασίας εδάφους (LST) μέσω Landsat-8 και R

Μελέτη Θερμικής Αστικής Νησίδας: Υπολογισμός θερμοκρασίας εδάφους (LST) μέσω Landsat-8 και RPhilip FayadBlockedUnblockFollowFollowingMar 28Οδηγός για δορυφορικές εικόνες Landsat-8Ο ακόλουθος οδηγός βασίζεται στις οδηγίες του Αμερικάνικου οργανισμού γεωλογικής επισκόπησης (USGS) και γίνεται χρήση των καναλιών 10 και 11 του θερμικού αισθητήρα (TIRS) του δορυφόρου Landsat-8.

Σε ποσοτικές αναλύσεις, δηλαδή εκεί όπου απαιτούνται μαθηματικές και στατιστικές αναλύσεις, συνίσταται η χρήση του καναλιού 10 διότι το κανάλι 11 είναι πολύ πιο επιρρεπές σε εσφαλμένες μετρήσεις ακτινών (Stray light).

Βήμα 1: Προετοιμασία δεδομένωνΑρχικά, δημιουργούμε έναν δωρεάν λογαριασμό ώστε να αποκτήσουμε προσβαση σε δορυφορικές εικόνες Landsat-8 απο ιστοσελίδα του Αμερικάνικου οργανισμού γεωλογικής επισκόπησης USGS (earthexplorer.

usgs.

gov) ως εξής.

Στην πρώτη καρτέλα (Search criteria) ορίζουμε τις παραμέτρους αναζήτησης, δηλαδή την περιοχή και την χρονική περίοδο ενδιαφέροντος.

Στην δεύτερη καρτέλα (Data sets) επιλέγουμε τα επιθυμητά δορυφορικά δεδομένα για κατέβασμα.

Σε αυτή τη περίπτωση επιλέγουμε δεδομένα τύπου Landsat Sentinel 8 OLI/TIRS C1 Level-1 απο τον φάκελο Landsat Collection 1 Level-1.

Οι δορυφορικές εικόνες Landsat-8 level 1C έχουν διορθωθεί γεωμετρικά αλλά δεν έχουν υποστεί ατμοσφαιρικές διορθώσεις (σε αντίθεση με τις εικόνες Landsat-8 level-2).

Για τον λόγο αυτό στο τελικό αποτέλεσμα υπάρχει η επίδραση της ατμόσφαιρας και διακρίνονται σύννεφα.

Μετά την αναζήτηση και το κατέβασμα των επιθυμητών εικόνων, ανοίγουμε το λογισμικό R και προχωρούμε στο επόμενο βήμα.

Περιβάλλον RΓια τους σκοπούς της ανάλυσης και υπολογισμού των τιμών Θερμότητας Εδάφους (LST) θα χρησιμοποιήσουμε μόνο τα κανάλια (μπάντες) 10 και 11 κάθε εικόνας όπου σε κάθε εικονοστοιχείο (pixel) τους περιέχονται ψηφιακές τιμές (DN) αποκρίσεων του αισθητήρα.

Οι τιμές αυτές πρέπει να μετατραπούν σε τιμές ακτινοβολίας (Radiance), ανακλαστικότητας (Reflectance) και εν τέλει σε τιμές LST (Land Surface Tempeature).

Χρησιμοποιούνται οι εξής πιο κάτω μαθηματικοί τύποι και όλοι οι υπολογισμοί γίνονται στο περιβάλλον της γλώσσας προγραμματισμού R.

Για την εισαγωγή και διαχείριση των δορυφορικών εικόνων στο περιβάλλον R, δείτε το βίντεο ή το PDF πιο κάτω.

Remote sensing using R: Viewing satellite images from the Rstudio consoleΑπαιτείται η εγκατάσταση των πακετων “Raster” και “rgdal”.

Βήμα 2: Μετατροπή Ψηφιακών Τιμών (Digital Number — DN) σε Ακτινοβολία (TOA Radiance) και ανακλαστικότητα (TOA Reflectance)Πηγή: https://www.

usgs.

gov/land-resources/nli/landsat/using-usgs-landsat-level-1-data-productΓια την επίλυση της εξίσωσης θα χρειαστούμε τις τιμές τών RADIANCE_MULT_BAND_x και RADIANCE_ADD_BAND_x από το αρχείο των μεταδεδομένων (metadata, _MTL), το οποίο δίδεται με κάθε εικόνα Landsat.

Οι τιμές αυτές μπορεί διαφέρουν απο εικόνα σε εικόνα, έτσι απαιτείται έλεγχος σε κάθε εικόνα.

Ανοίγουμε το αρχείο των μεταδεδομένων σε περιβάλλον λογισμικού κειμένου (π.

χ.

WordPad) και εντοπίζουμε τις τιμές.

Ο υπολογισμός σε περιβάλλον R γίνεται ως εξής.

#Values from MetafileRADIANCE_MULT_BAND_10 <- 3.

3420E-04RADIANCE_MULT_BAND_11 <- 3.

3420E-04RADIANCE_ADD_BAND_10 <- 0.

10000RADIANCE_ADD_BAND_11 <- 0.

10000#Load raster package and load band 10 & 11 into R (navigate to your image directory first)library(raster)band_10 <- raster("band10.

tif") #change image name accordinglyband_11 <- raster("band11.

tif") #change image name accordingly#Calculate TOA from DN:toa_band10 <- calc(band_10, fun=function(x){RADIANCE_MULT_BAND_10 * x + RADIANCE_ADD_BAND_10})toa_band11 <- calc(band_11, fun=function(x){RADIANCE_MULT_BAND_11 * x + RADIANCE_ADD_BAND_11})Σύμφωνα με την USGS, επόμενο βήμα είναι η διόρθωση γωνίας Ηλίου.

Ωστόσο αυτό δεν εφαρμόζεται στις περιπτώσεις των καναλιών 10 και 11, έτσι πηγαίνουμε στο επόμενο βήμα χωρίς εφαρμογή του πιο κάτω τύπου.

Πηγή: https://www.

usgs.

gov/land-resources/nli/landsat/using-usgs-landsat-level-1-data-productΓια την επίλυση της εξίσωσης θα χρειαστούμε τις τιμές τών REFLECTANCE_MULT_BAND_x και REFLECTANCE_ADD_BAND_x από το αρχείο των μεταδεδομένων (metadata), το οποίο δίδεται με κάθε εικόνα Landsat.

Αυτό το βήμα ΔΕΝ απαιτείται διότι σε αυτον τον οδηγό αξιοποιούνται τα κανάλια 10 και 11, στα οποία δεν υπάρχει η πληροφορία της ανακλαστικότητας στα μεταδεδομένα.

Βήμα 3: Μετατροπή τιμών ανακλαστικότητας σε τιμές θερμοκρασίας εδάφους (LST)Εν τέλει, η μετατροπή σε τιμές Θερμοκρασίας Εδάφους (LST) γίνεται ακολουθώντας τον εξής τύπο.

Πηγή: https://www.

usgs.

gov/land-resources/nli/landsat/using-usgs-landsat-level-1-data-productΓια την επίλυση της εξίσωσης θα χρειαστούμε τις τιμές τών K1_CONSTANT_BAND_10, K2_CONSTANT_BAND_10, K1_CONSTANT_BAND_11 και K2_CONSTANT_BAND_11 από το αρχείο των μεταδεδομένων (metadata), το οποίο δίδεται με κάθε εικόνα Landsat.

Ο υπολογισμός σε περιβάλλον R γίνεται ως εξής.

#Values from MetafileK1_CONSTANT_BAND_10 <- 774.

8853K1_CONSTANT_BAND_11 <- 480.

8883K2_CONSTANT_BAND_10 <- 1321.

0789K2_CONSTANT_BAND_11 <- 1201.

1442#Calculate LST in Kelvin for Band 10 and Band 11temp10_kelvin <- calc(toa_band10, fun=function(x){K2_CONSTANT_BAND_10/log(K1_CONSTANT_BAND_10/x + 1)})temp11_kelvin <- calc(toa_band11, fun=function(x){K2_CONSTANT_BAND_11/log(K1_CONSTANT_BAND_11/x + 1)})#Convert Kelvin to Celsius for Band 10 and 11temp10_celsius <- calc(temp10_kelvin, fun=function(x){x – 273.

15})temp11_celsius <- calc(temp11_kelvin, fun=function(x){x – 273.

15})#Export raster imageswriteRaster(temp10_celsius, "temp10_c.

tif")writeRaster(temp11_celsius, "temp11_c.

tif")*Για την μετατροπή των τιμών απο κλιμακα Kelvin σε κλίμακα Κελσίου, απαιτείται η αφαίρεση 273.

15 μονάδων ( function(x){x-273.

15} ).

Συνήθως και στις πλείστες μελέτες, για να υπολογιστεί το LST απαιτείται ο υπολογισμός των δεικτών NDVI και emissivity όπως φαίνεται στον οδηγό της ESA (2017), ωστόσο σε αυτή την μελέτη γίνεται η χρήση μόνο των καναλιών 10 και 11 του Θερμικού υπέρυθρου αισθητήρα (TIRS) του δορυφόρου Landsat-8 και δεν υπολογίζονται NDVI και Emissivity.

Βήμα 4: Απόδοση στο περιβάλλον QGISΣε αυτό το σημείο αποδίδουμε τα αποτελέσματα, ταξινομώντας τις τιμές της τελικής εικόνας (του καναλιού 10), στο λογισμικό QGIS.

Ένα παράδειγμα φαίνεται πιο κάτω.

Το τελικό αποτέλεσμα, τιμές σε κλίμακα κελσίου.

Τα σύννεφα αποτελούν “θόρυβο” και οι ακραίες τιμές απεικονίζονται με άσπρο χρώμα.

Για καλύτερα οπτικά αποτελέσματα και τα δυο αυτά στοιχεία πρέπει να μασκαριστούν και να αφαιρεθούν.

Ελπίζω το παρών άρθρο να ήταν βοηθητικό.

Για απορίες/εισηγήσεις/προτάσεις μην διστάσετε να επικοινωνήσετε μαζί μου μέσω της προσωπικής μου ιστοσελίδας: https://philipkfayad.

wixsite.

com/homeΠηγές/ΒιβλιογραφίαUSGS.

Using the USGS Landsat Level-1 Data Product.

(Accessed: 24/03/2019)Martin Šiklar (2016).

Calculation of Land Surface Temperature (LST) from Landsat-8 using R.

https://www.

gis-blog.

com/calculation-of-land-surface-temperature-lst-from-landsat-8-using-r/.

(Accessed: 24/03/2019).

EarthExplorer Platform.

https://earthexplorer.

usgs.

gov.

(Accessed: 24/03/2019)Department of the Interior U.

S.

Geological Survey, USGS (2018).

LANDSAT 8 (L8) DATA USERS HANDBOOK.

Επιπλέον ΥλικόΥπολογισμός LST χρησιμοποιώντας το ArcGIS: https://www.

youtube.

com/watch?v=uDQo2a5e7dMΥπολογισμός LST χρησιμοποιώντας το ENVI: https://gis.

stackexchange.

com/questions/137752/how-i-calculate-land-surface-temperature-using-landsat-8-in-envi-5-1Κανάλια/Μπάντες δορυφόρων Landsat: https://landsat.

usgs.

gov/what-are-band-designations-landsat-satellitesrLandsat8 — an R interface to Landsat8: http://terradue.

github.

io/rLandsat8/index.

html.

. More details

Leave a Reply