ESTUDIO DE CASOS: TEMPERATURA DE SUPERFICIE E ISLA CALOR URBANA
G174: FOTOINTERPRETACIÓN Y TELEDETECCIÓN
Domingo F. Rasilla y Francisco Conde
2023-11-08
1 INTRODUCCIÓN
El término temperatura de la superficie terrestre (conocida como LST por sus siglas en inglés) es muy utilizado en diferentes estudios científicos. Constituye una medida del calentamiento de la superficie terrestre, que recibe y reemite los rayos del sol, y depende en gran medida de las características y propiedades de sus diferentes superficies. Una de las superficies que experimenta un calentamiento diferencial son las ciudades, que suelen experimentar temperaturas más altas que las zonas circundantes. La causa principal son la acumulación de estructuras, como edificios, aceras o asfaltos, que absorben más calor y lo liberan más lentamente, al revés que superficies naturales como bosques, ríos o lagos. En el caso de la temperatura del aire, el efecto de isla de calor suele ser más acusado por las noches, llegando a registrarse en momentos puntuales diferencias superiores a 5ºC. En el caso de la temperatura de la superficie terrestre, las mayores diferencias también se observan por la noche, mientras que las diferencias diurnas dependen en cierta medida del tipo de superficies que rodean la ciudad. Así, en zonas áridas, es habitual que la ciudad constituya en realidad una “isla de frescor”, debido al elevado calentamiento de las superficies rurales.
En la actualidad la única forma de obtener una detallada, global, precisa (1a2◦C) y frecuente medida de LST es mediante el uso de las bandas termales infrarrojas (TIR, ver 7.4.4) de las diversas misiones espaciales.
En condiciones de cielo despejado, la medición mediante las bandas TIR recogen al menos tres componentes principales. Estos componentes se originan tanto en la atmósfera como en la superficie:
Radiación emitida desde la superficie.
Radiación TIR atmosférica reflejada por la superficie en el campo de visión del sensor.
Radiación TIR emitida en el trayecto mismo.
Para estimar LST a partir del TIR requiere desacoplar la radiación emitida de la superficie de los dos componentes atmosféricos y ajustar por la propia atenuación atmosférica e incluso la misma emisividad de la superficie terrestre (LSE). El valor LSE es previamente desconocido, así el cálculo LST es matemáticamente no resoluble. Entonces, o se debe conocer previo a la estimación o debe ser estimada simultáneamente a partir de modelos, restricciones y suposiciones.
En la actualidad las distintas misiones de EO proveen a los usuarios productos ya corregidos y modelados para dar una alta precisión de los valores estimados (< 2◦C).
2 DATOS
2.1 Ficheros vectoriales:
Para este ejercicio se utilizarán dos ficheros con formato vectorial (shapefiles):
Un fichero con los límites administrativos del área de trabajo.
Un fichero con la delimitación de las unidades administrativas (distritos).
2.2 Raster data (Landsat)
3 ÁLGEBRA RASTER
El objetivo de esta sección es obtener los valores de temperatura de la superficie terrestre y delimitar la superficie urbana de nuestro ámbito de estudio para determinar si existe una relación entre ambos (la conocida como isla de calor urbana).
3.1 Cálculo del índice NDVI.
Tras una visualización y análisis básico, comienza el procedimiento para el cálculo de la temperatura de la superficie terrestre. El procedimiento requiere generar varias capas adicionales en formato ráster.
La primer de estas capas es el índice NDVI. Como ya se ha expuesto previamente, la vegetación viva se puede representar con el NIR y las bandas rojas, ya que la clorofila se refleja en la longitud de onda NIR, pero se absorbe en la longitud de onda roja. Recuérdese que la fórmula es
\[NDVI= \frac{NIR-Red}{NIR+Red}\]
\[NDVI (Landsat 8) = \frac{B5 – B4} {B5 + B4}\]
\[NDVI (Landsat 4-7) = \frac{B4 – B3} {B4 + B3}\]
El cálculo del índice NDVI se realizará mediante una función específica, que ya se mencionó en capítulos anteriores.
NDVIfun <- function(NIR, Red) {
NDVI <- (NIR - Red) / (NIR + Red)
return(NDVI)
}Esta función podría ser reutilizada posteriormente en otros análisis.
Para ello sólo habría que crear y guardar un nuevo script de R, que
luego podría ser activado mediante la función source():
# source('insert file name')A continuación, se calculará el índice NDVI y se representará gráficamente
ndvi <- NDVIfun(imagen$NIR, imagen$red)
plot(ndvi,
col = rev(terrain.colors(10)),
main = "Landsat-NDVI")
Mapa con el NDVI
La distribución de valores de esta capa puede analizarse mediante un histograma.
hist(ndvi,
breaks = 40,
main = "Histograma NDVI",
xlim = c(-.3,.8))
Distribución del índice NDVI
Una posibilidad para identificar las zonas con vegetación es reclasificar los valores de NDVI por encima de un determinado valor. En este caso, se considerará que un valor de 0.3 es adecuado para este propósito, pero convendría sustentar esta cifra u otro similar a través de un razonamiento sustentado en las fuentes consultadas por el investigador.
veg <- reclassify(ndvi,
cbind(-Inf, 0.3, NA))
plot(veg,
main = 'Posible cubierta vegetal')
Mapa con la posible vegetación del área de estudio
A continuación, se puede visualizar la localización y extensión de las superficies verdes sobre la imagen en color verdadero.
plotRGB(imagen_rgb,
axes = TRUE,
stretch = "lin",
main = "Imagen en color verdadero")
plot(veg,
add=TRUE,
legend=FALSE)
Superposición de la vegetación sobre la imagen en color verdadero
3.2 Cálculo de la temperatura de la superficie terrestre
Hay diferentes métodos para calcular la temperatura de la superficie terrestre. A continuación se usará el desarrollado originalmente por Artis y Carnahan (1982), que ha sido reformulado recientemente por Guha et al. 2018 y Avdan y Jovanovska (2016). Este método comprende varios pasos
Procedimiento de cálculo de la LST
3.2.1 Obtención de la radiancia espectral en la parte superior de la atmósfera TOA
La radiancia espectral en la parte superior de la atmósfera (TOA) es la luz reflejada en la Tierra vista desde el satélite medida en unidades de radiancia. Su cálculo se lleva a cabo a partir de los números digitales que conforman cada banda (DN) usando la siguiente fórmula:
\[\lambda= Grescale * QCAL + Brescale\]
En esta ecuación, Grescale y Brescale representan la ganancia (REFLECTANCE_MULT_BAND) y el sesgo (REFLECTANCE_ADD_BAND) de la imagen, siendo QCAL el número digital (DN).
El cálculo se puede realizar manualmente. Grescale y Brescale están
disponibles en el archivo _MTL.txt que se incluye en la
escena descargada. Este fichero puede abrirse con el bloc de notas,
extrayéndose manualmente los valores correspondientes a la banda 10.
Imagen del fichero *MTL
La conversión de la banda 10 original en valores de radiancia en lo alto de la atmósfera se puede realizar fácilmente:
TOA_manual <- 2.0E-5 * lsatmask$LC08_L1TP_203023_20190513_20190521_01_T1_B10mask -0.1
TOA_manualUna segunda posibilidad es utilizar la función
RStoolbox::MTL() dentro del paquete
RStoolbox
library(RStoolbox)
MTL <- list.files("./datos/raster/",
pattern="*_MTL.txt",
ignore.case=TRUE,
full.names = TRUE)
readMTL<-readMeta(MTL)
readMTLA partir del objeto readMTL se extraen los valores
correspondientes a la Banda 10 y se procede al cálculo de la radiancia
TOA para esa banda.
offsetandgain <- subset(readMTL$CALRAD,
rownames(readMTL$CALRAD) == "B10_dn")
TOA <- offsetandgain$gain * lsatmask$LC08_L1TP_203023_20190513_20190521_01_T1_B10maskmask + offsetandgain$offset3.2.2 Conversión de TOA a temperatura de brillo.
La temperatura de brillo es la radiación emitida desde la parte superior de la atmósfera hasta el satélite en unidades de temperatura de un cuerpo negro equivalente. La conversión de radiancia TOA en temperatura de brillo \(T_b\) requiere la siguiente ecuación:
\[T_b=\frac{K_2}{ln((K_1/\lambda)+1)}\]
K1 (774.8853) y K2 (1321.0789) son las constantes de calibración previas al lanzamiento proporcionadas por USGS handbook.
Estos valores también pueden ser extraidos del fichero
_MTL.txt
Calidata <- as.data.frame(readMTL$CALBT)
# subconjunto de las filas
Calidata <- subset(Calidata,
rownames(Calidata) %in% "B10_dn")
# subconjunto las columnas
K1<-subset(Calidata, select=(K1))
K2<-subset(Calidata, select=(K2))
# get just the values out
K1<-K1$K1
K2<-K2$K2
# esto también funcionaría para K1
K1<-readMTL$CALBT$K1[1]
# para K2
K2<-readMTL$CALBT$K2[1]
Brighttemp<-(K2 / log((K1 / TOA) + 1))3.2.3 Cálculo de la emisividad
La emisividad de cada pixel requiere calcular la fracción que representa la vegetación en cada pixel, y esto a su vez necesita de los valores del índice NDVI. La fracción de vegetación es la relación entre el área con vegetación proyectada verticalmente y la extensión total de cada pixel.La fracción de vegetación de cada píxel se calcula mediante la siguiente ecuación:
\[F_v= \left( \frac{NDVI - NDVI_{min}}{NDVI_{max}-NDVI_{min}} \right)^2\]
Aquí, \(NDVI_{min}\) es el valor mínimo de NDVI (0,2) donde los píxeles se consideran tierra desnuda y \(NDVI_{max}\) es el valor en el que los píxeles se consideran vegetación saludable (0,5).
facveg <- (ndvi-0.2/0.5-0.2)^2La emisividad es la relación entre la energía de radiación absorbida y la energía de radiación entrante total en comparación con un cuerpo negro (que absorbería todo), siendo una medida de absoptividad. La emisividad se calcularía mediante la siguiente ecuación.
\[\varepsilon = 0.004*F_v+0.986\]
emiss <- 0.004*facveg + 0.986Por último, la temperatura de la superficie terrestre (LST) se obtendía siguiendo la ecuación de Weng et al. 2004 (también resumido en Guja et al. (2018) y Avdan y Jovanovska (2016)):
\[LST= \frac{T_b}{1+(\lambda \varrho T_b / (p))ln\varepsilon}\]
donde:
\[p= h\frac{c}{\varrho}\]
donde
- \(h\) que es la constante de Plank \(6.626 × 10^-34 Js\)
Plank <- 6.626*10e-34- \(c\) que es la velocidad de la luz en el vacío \(2.998 × 10^8 m/seg\)
c <- 2.998*10e8- \(\varrho\) que es la constante de Boltzmann de \(1.38 × 10^-23 J/K\)
Boltzmann <- 1.38*10e-23Con lo que el valor de \[p\] sería
p <- Plank*(c/Boltzmann)- Por último, \(\lambda\) es la longitud de onda efectiva para Landsat 8 banda 10.
lambda <- 1.09e-5Obtención de la LST mediante la siguiente ecuación.
LST <- Brighttemp/(1 +(lambda*Brighttemp/p)*log(emiss)) - 273.15
LST
plot(LST)
Mapa con la LST
3.3 Calculo del área urbana
La estimación de la superficie ocupada por superficies edificadas o en desarrollo de construcción se puede alcanzar mediante el cálculo de varios índices espectrales.
Uno de los más usados es el NDBI o Índice de Diferencia Normalizada Edificada (Zha et al., 2003). El cálculo del índice NDBI requiere de las bandas de análisis del infrarrojo a través de las bandas SWIR y NIR:
\[NDBI= \frac{SWIR - NIR}{SWIR + NIR}\]
En el caso de Landsat 8
\[NDBI = \frac{(B6 - B5)} {(B6 + B5)}\]
En el caso de Landsat 5 TM
\[NDBI = \frac{(B5 - B4)} {(B5 + B4)}\]
Al igual que en índices análogos el intervalo de valores resultantes oscila entre -1 y 1, donde valores negativos indican presencia de zonas con vegetación; valores intermedios comienzan a determinar zonas desnudas, cultivos en crecimiento o zonas o en fase de construcción, mientrs que valores positivos elevados indicaría áreas edificadas o infraestructuras antrópicas.
Otro índice similar es el UI o Urban Index de Kawamura, que ofrece grandes similitudes, aunque incorporando la banda VNIR:
\[UI= \frac{SWIR - VNIR}{SWIR + VNIR}\]
El cálculo del índice NDBI se puede realizar por dos vías diferentes:
- Manualmente
NDBI <-((lsatmask$LC08_L1TP_203023_20190513_20190521_01_T1_B6maskmask-
lsatmask$LC08_L1TP_203023_20190513_20190521_01_T1_B5maskmask)/
(lsatmask$LC08_L1TP_203023_20190513_20190521_01_T1_B6maskmask+
lsatmask$LC08_L1TP_203023_20190513_20190521_01_T1_B5maskmask))
plot(NDBI)- Mediante una función:
NDBIfunexample <- NDVIfun(lsatmask$LC08_L1TP_203023_20190513_20190521_01_T1_B6maskmask,
lsatmask$LC08_L1TP_203023_20190513_20190521_01_T1_B5maskmask)
plot(NDBI)4 ANÁLISIS DE LA RELACIÓN ENTRE EL ÁREA URBANA Y LA TEMPERATURA DE LA SUPERFICIE TERRESTRE
4.1 Cuantificación de la relación
La relación entre ambos parámetros se puede visualizar rápidamente mediante una nube de puntos.
plot(values(NDBI),
values(LST))
Relación entre la temperatura de superficie y el índice NDBI
Sin embargo, dado que se dispone de numerosos puntos, se produce un
problema conocido como superposición de gráficos. Una alternativa para
solucionarlo es aislar un subconjunto aleatorio de los mismos píxeles de
ambas capas ráster. Para hacerlo, es necesario volver a crear un
RasterStack con ambos conjuntos
computeddata <- stack(LST,NDBI) # Creación de un rasterStack.
random <- sampleRandom(computeddata, 500, cells=TRUE) # Extracción de una muestra aleatoria de 500 casos.
plot(random[,2], random[,3]) # Verificación del objeto `random`. 
Relación entre la temperatura de superficie y el índice NDBI (selección aleatoria)
Para obtener una visualización óptima, se necesitan los siguientes paquetes.
library(plotly)
library(tidyverse)
library(htmlwidgets)
library(ggplot2)La visualización de esos datos requiere su transformación en un dataframe sobre el que se trabaja con ggplot.
randomdf <- as.data.frame(random) # Conversión en un dataframe
names(randomdf) <- c("cell", "LST", "NDBI") # Renombramos las columnas del dataframe
heat<-ggplot(randomdf, aes(x = Urban, y = Temperature)) +
geom_point(alpha=2, colour = "#51A0D5") +
labs(x = "LST",
y = "NDBI",
title = "Relación entre LST y NDBI") +
geom_smooth(method='lm', se=FALSE)+
theme_classic()+
theme(plot.title = element_text(hjust = 0.5))
ggplotly(heat) # Dibujo interactivoA pesar de lo anterior, se podría trazar todo el conjunto de datos recurriendo a un gráfico de hexágonos
computeddatadf <- as.data.frame(computeddata)
names(computeddatadf) <- c("Temperature", "Urban")
hexbins <- ggplot(computeddatadf, aes(x=Urban, y=Temperature) ) +
geom_hex(bins=100, na.rm=TRUE) +
labs(fill = "Número de casos según bin")+
geom_smooth(method='lm', se=FALSE, size=0.6)+
theme_bw()
ggplotly(hexbins)Además de su representación gráfica, conviene cuantificar la relación visual observada. Para ello se puede utilizar el coeficiente de correlación de Pearson.
stat.sig <- cor.test(computeddatadf$Temperature, computeddatadf$Urban,
use = "complete.obs", # Evitar que NA impidan su cálculo
method = c("pearson"))
stat.sigA partir de este procedimiento se obtienen una serie de estadísticos adicionales
- t: valor del estadístico del t-test, cuyo valor crítico se puede obtener así. Si el valor de t obtenido es mayor que el valor crítico se puede señalar que existe una relación estadísticamente significativa entre ambos fenómenos.
abs(qt(0.05/2, 198268))- df: número de grados de libertad.
length(computeddatadf$Temperature)
removena <- na.omit(computeddatadf) # Eliminamos algunos datos con NA
length(removena$Temperature)
length(removena$Urban)p-value: informa si se puede rechazar la hipótesis hula cuando p<0.05, lo que supone una relación estadísticamente significativa entre ambos conjuntos de datos (hay un 95% de posibilidades de que la relación sea real).
cor: coeficiente de correlación de Pearson.
Intervalo de confianza del 95 %: confianza al 95% de que el coeficiente de correlación de la población se sitúe dentro de este intervalo.
Dado que p-value es < 0.05 se podría sostener que la temperatura de la superficie terrestre aumenta a medida que aumenta/densifica la zona urbanizada (pero habría que explicar por qué).
4.2 Cálculo de la temperatura según distritos
Para este apartado es necesario cargar paquetes adicionales
library(dplyr)
library(sf)En primer lugar, se importará el fichero con los distritos, reproyectándolo al mismo CRS para ajustarlo al fichero con los límites del área de trabajo.
UK_LSOA <- st_read("./datos/vectorial/LSOA/Lower_Super_Output_Area_(LSOA)_IMD_2019__(OSGB1936).shp")
UK_LSOA <- UK_LSOA %>%
st_transform(., 32630)A continuación, se procederá a la extracción de los valores medios de
LST dentro de cada unidad espacial (distritos) mediante la herramienta
raster::extract().
LST_per_LSOA <- raster::extract(LST, # Objeto raster de entrada
manchester_LSOA, # Objeto vectorial de entrada
fun=mean, # función para resumir la información de los píxeles dentro de cada unidad espacial (municipios): promedio.
na.rm=TRUE, # Se ignora los valores NA.
df=TRUE) # Formato -dataframe- de los resultados. Una vez creado, el dataframe con los resultados (promedio de la LST
según cada municipio) será unido al dataframe original. Para ello se
requiere que ambos dataframes tengan una variable común idéntica, en
este caso el ID. En primer lugar, se añadirá el identificador del
dataframe original al nuevo dataframe con los valores de LST, lo que
permitirá realizar una unión left_join() entre ambos. Para
ello:
Se añaden los ID al nuevo dataframe.
LST_per_LSOA$FID <- manchester_LSOA$FIDSe procede a una Left join, es decir, a añadir los
ID.
manchester_LSOA_temp <- merge(x=manchester_LSOA,
y= LST_per_LSOA,
by="FID",
all.x=TRUE)
names(manchester_LSOA_temp) <- Añadir la temperatura media según polígonos al fichero vectorial
manchester_LSOA_temp <- manchester_LSOA %>%
left_join(.,
LST_per_LSOA,
by="FID")%>%
dplyr::rename(temp=layer)Si con los comandos anteriores se ha obtenido la temperatura promedio
en cada municipio, la cantidad de superficie construida en cada
municipio es desconocida. Para ello, consideraremos que cualquier valor
del índice NDBI por enicma de 0 significa que todo el pixel es urbano.
La herramienta raster::extract() volverá a ser utilizada
para obtener el % de superficie urbana dentro de cada unidad espacial.
Se define como urbano todo pixel con un NDBI > 0.
NDBI_urban<- NDBI > 0 A continuación, se calculará la suma de todos los píxeles urbanos (NDBI mayor que 0) dentro de cada unidad espacial.
NDBI_urban_per_LSOA <- raster::extract(NDBI_urban,
manchester_LSOA,
fun=sum,
na.rm=TRUE,
df=TRUE)También se calculan todos los píxeles (independientemente de si son urbanos o no) dentro de cada unidad espacial:
NDBI_per_LSOA_cells <- raster::extract(NDBI_urban,
manchester_LSOA,
cellnumbers=TRUE,
na.rm=TRUE,
df=TRUE)Recuento de los píxeles según unidades admnistrativas:
NDBI_per_LSOA2_cells <- NDBI_per_LSOA_cells %>%
count(ID)Adición del ID al área urbana
NDBI_urban_per_LSOA$FID <- manchester_LSOA$FIDAñadir el ID al dataframe con el número de píxeles según unidades administrativas:
NDBI_per_LSOA2_cells$FID <- manchester_LSOA$FIDFusión de ambos dataframes utilizando la variable
FID:
Urban_info_LSOA <- NDBI_urban_per_LSOA %>%
left_join(.,
NDBI_per_LSOA2_cells,
by="FID")Se elimina las variables que no se necesitan y se renombran las restantes
Urban_info_LSOA_core_needed <- Urban_info_LSOA %>%
dplyr::rename(urban_count=layer,
LSOA_cells=n) %>%
dplyr::select(urban_count,
LSOA_cells,
FID) %>%
dplyr::mutate(percent_urban=urban_count/LSOA_cells*100)Elaboración final del fichero vectorial con los datos de temepratura y % de superficie urbana según distrito.
manchester_LSOA_temp_urban <- manchester_LSOA_temp %>%
left_join(.,
Urban_info_LSOA_core_needed,
by="FID")4.3 Cartografía
Una vez disponemos del fichero vectorial con los datos de temperatura y % de área urbana en cada distrito, podemos incluir algunos nombres de lugares obtenidos de Open Street Map.
Places <- st_read(here::here("prac8_data",
"greater-manchester-latest-free.shp",
"gis_osm_places_free_1.shp")) %>%
st_transform(., 32630)
manchester_Places <- Places[manchester_boundary,]%>%
filter(fclass=="city")A continuación, se elaborará un mapa, usando el paquete
tmap. Recuerde agregar un título en Rmarkdown y incluya el
argumento fig.cap="caption here" en el encabezado del
fragmento de código. Este primer fragmento de código aumenta el tamaño
de la caja del mapa, para que se pueda incluir una flecha indicando el
norte que no se superponga al propio mapa (veáse https://www.jla-data.net/eng/adjusting-bounding-box-of-a-tmap-map/).
bbox_new <- st_bbox(manchester_LSOA_temp_urban) # "Bounding box" actual
yrange <- bbox_new$ymax - bbox_new$ymin # Rango de valores y
bbox_new[4] <- bbox_new[4] + (0.1 * yrange) # ymax - parte superior
bbox_new[2] <- bbox_new[2] - (0.1 * yrange) # ymin - parte inferiorCargamos el paquete necesario y establecemos el modo de visualización.
library(tmap)
tmap_mode("plot")Establecimiento del nuevo “bbox” y eliminación del bbox=bbox_new para ver la diferencia
tm1 <- tm_shape(manchester_LSOA_temp_urban, bbox = bbox_new) +
tm_polygons("temp",
palette="OrRd",
legend.hist=TRUE,
title="Temperatura")+
tm_shape(manchester_Places, bbox=bbox_new)+
tm_dots(size=0.1, col="white")+
tm_text(text="name", size=0.75, ymod=-0.5, col="white", fontface = "bold")+
#tm_legend(show=FALSE)+
tm_layout(frame=FALSE,
legend.outside=TRUE)+
tm_compass(type = "arrow", size=1, position = c("left", "top")) +
tm_scale_bar(position= c("left", "bottom"), breaks=c(0,2,4), text.size = .75)
#tm_credits("(a)", position=c(0,0.85), size=1.5)
tm1
Temperatura de la superficie terrestre según distritos
Este ejericicio también podría repetirse para el área urbana y
trazarse en una figura, tal vez con un mapa insertado como lo hicimos
anteriormente. Pero ahora también podemos trazar dos mapas de coropletas
uno encima del otro, mediante el paquete biscale(). Sin
embargo, ahora debemos usar ggplot en lugar de
tmap, puedes forzar tmap para hacerlo, pero la
característica aún está en desarrollo.
4.4 Mapa bivariado (optional)
Si bien los mapas bivariados son muy sugerentes, no dicen nada sobre los valores de datos reales, simplemente dividen los datos en clases según un determinado criterio (jenks, igual número de casos, etc). Es necesario, por tanto, establecer algunos criterios adicionales:
Un mapa bivariado central de LSOA dentro de Manchester
El bivariado central tendrá los límites del MSOA y algunos topónimos.
Dos gráficos adicionales mostrando la distribución de los datos.
library(biscale)
library(cowplot)
library(sysfonts)
library(extrafont)
library(showtext) # Fuentes adicionales
#font_add_google("Lato", regular.wt = 300, bold.wt = 700) # I like using Lato for data viz (and everything else...). Open sans is also great for web viewing.
showtext_auto()Creación de las clases
data <- bi_class(manchester_LSOA_temp_urban, x = temp, y = percent_urban, style = "jenks", dim = 3)Mapa con ggplot
map <- ggplot() +
geom_sf(data = data, mapping = aes(fill = bi_class), color=NA, lwd = 0.1, show.legend = FALSE) +
bi_scale_fill(pal = "DkViolet", dim = 3) +
geom_sf(data = manchester_MSOA, mapping = aes(fill=NA), color="black", alpha=0, show.legend = FALSE)+
geom_sf(data=manchester_Places, mapping=aes(fill=NA), color="white", show.legend = FALSE)+
geom_sf_text(data=manchester_Places, aes(label = name, hjust = 0.5, vjust = -0.5),
nudge_x = 0, nudge_y = 0,
fontface = "bold",
color = "white",
show.legend = FALSE,
inherit.aes = TRUE)+
labs(
title = "",
x="", y=""
) +
bi_theme()
legend <- bi_legend(pal = "DkViolet",
dim = 3,
xlab = "Temperatura ",
ylab = "% de superficie urbana",
size = 8)
credit<- ("Temperatura de superficie y área urbana, imagen del 13/05/2022")
Relación entre la temperatura de superficie y el índice NDBI según barrios
Combinación de mapa con la leyenda.
finalPlot <- ggdraw() +
draw_plot(map, 0, 0, 1, 1) +
draw_plot(legend, 0.1, 0.1, 0.2, 0.2)
#draw_text(credit, 0.68, 0.1, 0.2, 0.2, size=10)
finalPlotUna vez dibujado el gráfico principal, hay que hacer lo propio con los gráficos laterales.
urban_box<-ggplot(data, aes(x=bi_class, y=percent_urban, fill=bi_class)) +
geom_boxplot()+
scale_fill_manual(values=c("#CABED0", "#BC7C8F", "#806A8A", "#435786", "#AE3A4E", "#77324C", "#3F2949", "#3F2949"))+
labs(x="Clases bivariadas (temp, urban)",
y="% superficie urbana")+
theme_light()+
theme(legend.position="none") # Remove legend
temp_violin<-ggplot(data, aes(x=bi_class, y=temp, fill=bi_class))+
geom_violin()+
scale_fill_manual(values=c("#CABED0", "#BC7C8F", "#806A8A", "#435786", "#AE3A4E", "#77324C", "#3F2949", "#3F2949"))+
labs(x="",
y="Temperatura")+
guides(fill=guide_legend(title="Clase"))+
theme_light()+
theme(legend.position="none") # Remove legendÚnalos todos en dos pasos: haga las tramas secundarias y luego únalas a la trama principal.
side <- plot_grid(temp_violin, urban_box, labels=c("B","C"),label_size = 12, ncol=1)
all <- plot_grid(finalPlot, side, labels = c('A'), label_size = 12, ncol = 2, rel_widths = c(2, 1))
Ahora bien, esto parece un poco basura en la página. Por lo general,
es posible que quieras exportar el mapa a un .png para
usarlo en otro lugar, pero necesitarás cambiar su tamaño. Me gusta
trazar el mapa en la consola (escriba todo en la consola y luego
presione Intro), haga clic en exportar en la pestaña de trazados, cambie
la imagen y anote los números y luego introdúzcalos a continuación.
all
dev.copy(device = png, filename = here::here("prac8_images", "bivariate.png"),
width = 687,
height = 455)
dev.off()4.5 Consideraciones finales
Si desea explorar más a fondo este tipo de análisis, deberá considerar lo siguiente:
Otros métodos para extraer temperatura de datos Landsat (véase apartado siguiente)
Validación de su capa de temperatura (por ejemplo, datos de la estación meteorológica)
La fórmula utilizada para calcular la emisividad — hay muchas
El uso de datos satelitales sin procesar en lugar de eliminar los efectos de la atmósfera. En esta práctica solo hemos utilizado índices espaciales relativos (por ejemplo, NDVI). Sin embargo, si utilizara métodos alternativos, podría ser más apropiado utilizar datos de reflectancia de la superficie (también proporcionados por USGS).
…O ejecutar el análisis con diferentes datos y métodos, por ejemplo
Datos procedentes de MODIS (datos diarios de LST https://modis.gsfc.nasa.gov/data/dataprod/mod11.php)
Métodos: clasificación de usos de suelo supervisada o no supervisada (veáse https://rspatial.org/rs/5-supclassification.html). Se puede clasificar una imagen en varias clases de cobertura terrestre y explorar su relación con la temperatura.
4.6 References
Avdan, U. and Jovanovska, G., 2016. Algorithm for automated mapping of land surface temperature using LANDSAT 8 satellite data. Journal of Sensors, 2016.
Guha, S., Govil, H., Dey, A. and Gill, N., 2018. Analytical study of land surface temperature with NDVI and NDBI using Landsat 8 OLI and TIRS data in Florence and Naples city, Italy. European Journal of Remote Sensing, 51(1), pp.667-678.
Weng, Q., Lu, D. and Schubring, J., 2004. Estimation of land surface temperature–vegetation abundance relationship for urban heat island studies. Remote sensing of Environment, 89(4), pp.467-483.
Young, N.E., Anderson, R.S., Chignell, S.M., Vorster, A.G., Lawrence, R. and Evangelista, P.H., 2017. A survival guide to Landsat preprocessing. Ecology, 98(4), pp.920-932.
Zha, Y., Gao, J. and Ni, S., 2003. Use of normalized difference built-up index in automatically mapping urban areas from TM imagery. International journal of remote sensing, 24(3), pp.583-594.
5 CÁLCULO DE LA TEMPERATURA DE SUPERFICIE SEGÚN EL PAQUETE R LST
# library(raster)
library(LST)
set.seed(2)5.1 Datos de partida
Banda 10 (TIR 1)
Landsat_10 <- raster::raster(ncol=100, nrow=100)
raster::values(Landsat_10) = runif(10000, min=27791, max=30878)Banda 11 (TIR 2)
Landsat_11 <- raster::raster(ncol=100, nrow=100)
raster::values(Landsat_11) = runif(10000, min=25686, max=28069)Banda 4 (rojo)
red <- raster::raster(ncol=100, nrow=100)
raster::values(red) = runif(10000, min=0.1, max=0.4)Banda 5 Infrarrojo cercano (NIR)
NIR <- raster::raster(ncol=100, nrow=100)
raster::values(NIR) = runif(10000, min=0.1, max=0.6)5.2 Cálculos previos necesarios
Temperatura at sensor o de la temperatura de brillo
BTem <- BT(Landsat_10 = Landsat_10,
Landsat_11 = Landsat_11)NDVI
ndvi <- NDVI(Red = red, NIR = NIR)Emisividad de la superficie terrestre
# Método de Skokovic et al. 2014
e_skokovic <- E_Skokovic(red = red, # red band of remote sensing imagery
NDVI = ndvi, # NDVI calculated from remote sensing imagery
band = "band 11") # Especificación de la banda térmica de Landsat 8 (puede ser "band 10" o "band 11").
# Método de Sobrino et al. 2008
E_Sobrino(red = red, NDVI = ndvi)
# Método de Valor y Caselles 1996
E_Valor(NDVI = ndvi)
# Método de Van de Griend y Owe 1993
E_VandeGriend(NDVI = ndvi)
# Método de Yu et al. 2014
E_Yu(red = red, NDVI = ndvi, band = "band 11")Transmitancia atmosférica (cantidad de energía que atraviesa un cuerpo en la unidad de tiempo (potencia). Esta función calcula este parámetro utilizando la temperatura del aire (To, °C) y la humedad relativa (RH, %) de la fecha y del lugar en el que se tomó la imagen Landsat.
tau <- tau(To = 26,
RH = 42,
band = "band 11") # Especificación de la banda térmica de Lansat 8 usada (puede ser "band 10" o "band 11").Proporción de la vegetación (fracción de cubierta verde) extraída a partir de NDVI.
Pv(NDVI = ndvi,
minNDVI = 0.2, # Ref. Sobrino et al. 2004
maxNDVI = 0.5) # idem** Temperatura atmosférica media** (ºK)
Ta <- Ta(To = 26, # Temperatura del aire de la fecha y el lugar del que se tomó la imagen
mod = "Mid-latitude Winter Region") # Modelos a utilizar: "USA 1976 Standard", "Tropical Region", "Mid-latitude Summer Region" or "Mid-latitude Winter Region"5.3 Cálculo de la temperatura de superficie
Existen diferentes procedimientos para su cálculo
Mediante el algorítmo “mono window”.
BTemp <- raster::raster(ncol=100, nrow=100)
raster::values(BTemp) = runif(10000, min=298, max=305)
Ts_MWA <- MWA(BT = BTemp, # Raster* object, brightness temperature
tau = tau, # Atmospheric transmittance
E = e_skokovic, # Raster* object, Land Surface Emissivity calculated according to one of the methods
Ta = Ta) # Mean atmospheric temperature (K) of the date when Landsat passed over the study area
plot(Ts_MWA)Mediante el método de la ecuación de transferencia radiativa.
Ts_RTE <- RTE(TIR = TIR, # Landsat band 10 or 11
tau = 0.86, # Atmospheric transmittance
E = E, # Land Surface Emissivity calculated according to
dlrad = 2.17, # Downwelling radiance calculated from https://atmcorr.gsfc.nasa.gov/
ulrad = 1.30, # Upwelling radiance calculated from https://atmcorr.gsfc.nasa.gov/
band = "band 11") # A string specifying which Landsat 8 thermal band to use (can be "band 10" or "band 11")Mediante el algorítmo “single channel”.
Ts_SCA <- SCA(TIR = TIR, # Landsat band 10 or 11
BT = BT, # at-Sensor Temperature or brightness temperature
tau = 0.86,
E = E,
dlrad = 2.17,
ulrad = 1.30,
band = "band 11")Mediante el algorítmo “split-window”.
Ts_SWA <- SWA(TIR_10 = TIR_10, # Banda 10 Landsat
TIR_11 = TIR_11, # Banda 11 Landsat
tau_10 = 0.86, # Transmitancia atmosférica correpondiente a la banda 10 de Landsat.
tau_11 = 0.87, # Transmitancia atmosférica correpondiente a la banda 11 de Landsat.
E_10 = E_10, # Emisividad de la superficie (banda 10 de Landsat) según Skokovic et al. 2014 o Yu et al. 2014
E_11 = E_11) # Emisividad de la superficie (banda 11 de Landsat) según Skokovic et al. 2014 o Yu et al. 20146 BUILT UP INDEX CON SPATIAL ECHO
Hay muchos índices para el análisis de superficies construidas por el ser humano: los más comunes son el Normalized Difference Built-up Index (NDBI), Built-up Index (BU), Urban Index (UI), Index-based Built-up Index (IBI) y el Enhanced Built-up and Bareness Index (EBBI). Cada uno de ellos tiene una fórmula específic, y su propio método de cálculo. Las áreas urbanizadas y el suelo desnudo reflejan más SWIR que NIR, mientras que el agua no se refleja en el espectro infrarrojo. En el caso de la superficie verde, la reflexión de NIR es más alta que el espectro SWIR.
Un índice sencillo es el NDBI o Índice de Diferencia Normalizada Edificadam, que requiere de las bandas de análisis del infrarrojo a través de las bandas SWIR y NIR. El indicador se obtiene con la con la sencilla relación:
\[NDBI = \frac{SWIR-NIR} {SWIR+NIR}\] \[NDBI (Landsat 4-7) = \frac{B5 – B4} {B5 + B4}\] \[NDBI (Landsat 8) = \frac{B6 – B5} {B6 + B5}\] Al igual que en índices análogos (como el NDWI o el NDVI) el intervalo de valores resultantes oscila entre -1 y 1, donde aquellos valores de tendencia negativa indican presencia de zonas con vegetación. Valores intermedios comienzan a determinar zonas desnudas, cultivos en crecimiento o zonas o en fase de construcción a medida que adquieren valores de tendencia positivos elevados para indicar zonas territoriales con coberturas de suelo edificadas o infraestructuras antrópicas.
Este índice puede servir para una posterior clasificación supervisada, que ayudará a descubrir e identificar zonas de vegetación, zonas desnudas y zonas edificadas, localizando los núcleos urbanos e infraestructuras antrópicas como carreteras y autovías o aeropuertos.
Para obtener un mejor resultado se puede usar el Build-up Index (BI), que combina el NDBI y NDVI.
\[NDBI = NDBI - NDVI\]
Una variante es el Index based Built up Index (IBI), que hace uso de los conocidos índices SAVI y del MNDWI, según la siguiente expresión:
\[IBI = \frac{NDBI−(SAVI+MNDWI)/2} {NDBI + (SAVI+MNDWI)/2}\]
Finalmente, cabe resaltar el UI o Urban Index de Kawamura, que incorpora la banda VNIR:
\[UI = \frac{SWIR-VNIR} {SWIR+VNIR}\]
El paquete SpatialEcho incluye una función para calcular el Built-up index.
La sintaxis será la siguiente
built.index( green, red, nir, swir1, swir2, L = 0.5, method = c(“Bouhennache”, “Zha”, “Xu”) )
Argumentos: green Green band (0.53 - 0.59mm), landsat 5&7 band 3, OLI (landsat 8) band 3 red Red band (0.636 - 0.673mm), landsat 5&7 band 3, OLI (landsat 8) band 4 nir Near infrared band (0.851 - 0.879mm) landsat 5&7 band 4, OLI (landsat 8) band 5 swir1 Short-wave infrared band 1 (1.566 - 1.651mm), landsat 5&7 band 5, OLI (landsat 8) band 6 swir2 Short-wave infrared band 2 (2.11 - 2.29mm), landsat 5&7 band 7, OLI (landsat 8) band 7 L Factor L en el índice SAVI. method Método de cálculo. Las opciones son “Bouhennache”, “Zha”, “Xu”
Esta función calcula el índice siguiendo tres métodos diferentes:
Bouhennache es un nuevo método que utiliza una porción más grande del VIR/NIR siguiendo las bandas OLI (((b3+b4+b7)-b6)/3) / (((b3+b4+b7)+b6)/3)
Zha es el método de relación de banda original usando TM5 ndbi = (b5 - b4) / (b5 + b4)
Xu es una modificación para eliminar el ruido usando ETM+7 (ndbi-((savi-nndwi)/2) / (ndbi+((savi-nndwi)/2)
En general, el agua tiene los valores más altos donde se producirán áreas edificadas en la parte media de la distribución. Dado que el índice de Bouhennache et al (2018) explota una porción más grande del espectro visible (Vis) e infrarrojo, la vegetación se producirá con los valores más bajos y los estériles exhibirán valores mayores que la vegetación y valores más bajos que las áreas edificadas.
library(raster)
library(RStoolbox)
library(spatialEco)data(lsat)
lsatCálculo de la radiancia en lo alto de la atmósfera ((W/(m^2 * srad * μm)))
lsat <- radCor(lsat,
metaData = readMeta(system.file("external/landsat/LT52240631988227CUB02_MTL.txt", package="RStoolbox")),
method = "apref")
plotRGB(lsat,
r=3, g=2, b=1,
scale=1.0, stretch="lin")Cálculo según el método de Bouhennache et al. (2018; requiere las bandas green, red, swir1 and swir2)
bouh <- built.index(red = lsat[[3]],
green = lsat[[4]],
swir1 = lsat[[5]],
swir2 = lsat[[7]])
plotRGB(lsat, r=6,g=5,b=2, scale=1, stretch="lin")
plot(bouh,
legend=FALSE,
col=rev(terrain.colors(100, alpha=0.35)),
add=TRUE )Cálculo según Zha et al. (2003; requiere nir and swir1)
zha <- built.index(nir = lsat[[4]],
swir1 = lsat[[5]],
method = "Zha")
plotRGB(lsat,
r=6,g=5,b=2,
scale=1,
stretch="lin")
plot(zha, legend=FALSE, col=rev(terrain.colors(100, alpha=0.35)), add=TRUE )Cálculo segúnel método de Xu (2008), que es una modificación normalizada del me´todo de Zha (requiere las bandas green, red, nir and swir1)
xu <- built.index(green= lsat[[3]],
red = lsat[[3]],
nir = lsat[[4]],
swir1 = lsat[[5]],
method = "Xu")
plotRGB(lsat, r=6,g=5,b=2, scale=1, stretch="lin")
plot(xu, legend=FALSE, col=rev(terrain.colors(100, alpha=0.35)), add=TRUE )