PROGRAMACIÓN EN TIGS CON R - TP FINAL - TEMÁTICA: DATOS RASTER
Caterina Tibaldo Stralla
2024-07-24
En el presente trabajo se utilizará una imagen satelital, a la que se le harán diferentes tipos de análisis, para determinar las cualidades superficiales del suelo. Como área de estudio, se escogió una zona exclusivamente agrícola del Departamento Paraná, Provincia de Entre Ríos. Al momento de proyectar obras viales, sobre caminos productivos, es necesario que se caracterize la zona en torno a los distintos tipos de producción, la existencia de humedales o arroyos, y se identifique la cantidad de vegetación sobre la traza del camino, entre otros.
En primera medida, para relizar el presente trabajo práctico, se descargó una imagen satelital que incluye al área de estudio escogida. A través de la página de USGS, se obtuvó una imagen Landat 8. Luego, se cargaron las siguientes librerías, para comenzar a procesar los datos.
library(raster)## Warning: package 'raster' was built under R version 4.1.3## Loading required package: sp## Warning: package 'sp' was built under R version 4.1.3library(sf)## Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 7.2.1; sf_use_s2() is TRUEA continuación, se irán detallando los pasos que se fueron ejecutando, junto a los distintos análisis y conclusiones:
1) Se cargan las bandas azul, verde y roja de la imagen Landat para analizarlas
Banda Azul
b2<- raster ('C:/Users/admin/Documents/Caterina/Maestria TIG/R/Final/Imagen/LC08_L2SP_226082_20240714_20240719_02_T1_SR_B2.tif')Banda Verde
b3<- raster ('C:/Users/admin/Documents/Caterina/Maestria TIG/R/Final/Imagen/LC08_L2SP_226082_20240714_20240719_02_T1_SR_B3.tif')Banda Roja
b4<- raster ('C:/Users/admin/Documents/Caterina/Maestria TIG/R/Final/Imagen/LC08_L2SP_226082_20240714_20240719_02_T1_SR_B4.tif')2)- Se corrobora información de las mismas
b2## class : RasterLayer
## dimensions : 7691, 7631, 58690021 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 146085, 375015, -3630315, -3399585 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=21 +datum=WGS84 +units=m +no_defs
## source : LC08_L2SP_226082_20240714_20240719_02_T1_SR_B2.TIF
## names : LC08_L2SP_226082_20240714_20240719_02_T1_SR_B2b3## class : RasterLayer
## dimensions : 7691, 7631, 58690021 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 146085, 375015, -3630315, -3399585 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=21 +datum=WGS84 +units=m +no_defs
## source : LC08_L2SP_226082_20240714_20240719_02_T1_SR_B3.TIF
## names : LC08_L2SP_226082_20240714_20240719_02_T1_SR_B3b4## class : RasterLayer
## dimensions : 7691, 7631, 58690021 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 146085, 375015, -3630315, -3399585 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=21 +datum=WGS84 +units=m +no_defs
## source : LC08_L2SP_226082_20240714_20240719_02_T1_SR_B4.TIF
## names : LC08_L2SP_226082_20240714_20240719_02_T1_SR_B4Se puede observar que todas las bandas, son de clase ráster, poseen la misma dimensión, cantidad de celdas y resolución, y su sistema de coordenadas de referencia se corresponde con la EPSG 32621 (WGS 84 - UTM ZONA 21N). En el caso de que se desee realizar un estudio puntual de las características de la imagen, se pueden ir ejecutando diferentes códigos que resuelven la solicitud.
ncell(b2)## [1] 58690021dim(b2)## [1] 7691 7631 1nlayers(b2)## [1] 1Y en caso, de que se desee comparar solamente la extensión, dimensión, proyección y resolución de diferentes bandas, se puede ejecutar el siguiente código.
compareRaster(b2, b3, b4)## [1] TRUE3) Se procede a graficar las bandas
plot(b2,
zlim=c(0,15000),
main= 'Banda azul en ND')
plot(b3,
zlim=c(0,15000),
main= 'Banda verde en ND')
plot(b4,
zlim=c(0,15000),
main= 'Banda roja en ND')
4) Se crea un RasterStack a partir de los archivos de las bandas
Con el fin de realizar diferentes análisis sobre un único ráster, se unieron las 7 bandas que componen la imagen.
imagen<- paste0('C:/Users/admin/Documents/Caterina/Maestria TIG/R/Final/Imagen/LC08_L2SP_226082_20240714_20240719_02_T1_SR_B', 1:7, ".tif")
imagenlandsat<- stack(imagen)
imagenlandsat## class : RasterStack
## dimensions : 7691, 7631, 58690021, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 146085, 375015, -3630315, -3399585 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=21 +datum=WGS84 +units=m +no_defs
## names : LC08_L2SP//2_T1_SR_B1, LC08_L2SP//2_T1_SR_B2, LC08_L2SP//2_T1_SR_B3, LC08_L2SP//2_T1_SR_B4, LC08_L2SP//2_T1_SR_B5, LC08_L2SP//2_T1_SR_B6, LC08_L2SP//2_T1_SR_B75) Se cambian los nombres de las bandas y se grafica la imagen
names(imagenlandsat)<- c('B1', 'B2', 'B3', 'B4', 'B5', 'B6', 'B7')
names(imagenlandsat)## [1] "B1" "B2" "B3" "B4" "B5" "B6" "B7"plot(imagenlandsat)
6) Se realiza un histograma de la imagen para cada banda y una matriz de correlación
hist(imagenlandsat)## Warning in .hist1(raster(x, y[i]), maxpixels = maxpixels, main = main[y[i]], :
## 0% of the raster cells were used. 100000 values used.
## Warning in .hist1(raster(x, y[i]), maxpixels = maxpixels, main = main[y[i]], :
## 0% of the raster cells were used. 100000 values used.
## Warning in .hist1(raster(x, y[i]), maxpixels = maxpixels, main = main[y[i]], :
## 0% of the raster cells were used. 100000 values used.
## Warning in .hist1(raster(x, y[i]), maxpixels = maxpixels, main = main[y[i]], :
## 0% of the raster cells were used. 100000 values used.
## Warning in .hist1(raster(x, y[i]), maxpixels = maxpixels, main = main[y[i]], :
## 0% of the raster cells were used. 100000 values used.
## Warning in .hist1(raster(x, y[i]), maxpixels = maxpixels, main = main[y[i]], :
## 0% of the raster cells were used. 100000 values used.
## Warning in .hist1(raster(x, y[i]), maxpixels = maxpixels, main = main[y[i]], :
## 0% of the raster cells were used. 100000 values used.pairs(imagenlandsat, main='Matriz de correlacion')
A través del histograma, se puede apreciar cuales son los niveles de energía y como se distribuye la misma en las diferentes bandas. Y con la matriz de correlación, cual es el grado de correspondencia espacial entre bandas.
7) Se analiza cual es el sistema de coordenadas de la imagen
crs(imagenlandsat)## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=utm +zone=21 +datum=WGS84 +units=m +no_defs
## WKT2 2019 representation:
## PROJCRS["unknown",
## BASEGEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]]],
## CONVERSION["UTM zone 21N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-57,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]],
## ID["EPSG",16021]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]8) Se procede a cargar la capa vectorial, correspondiente al área de estudio
El área de estudio agrupa a los Centros Rurales de Población de Aldea San Antonio, Paso de la Arena y Quebracho. Para poder recortar la imagen satelial, es necesario que ambas capas se encuentren en el mismo sistema de coordenadas.
crp<- st_read("C:/Users/admin/Documents/Caterina/Maestria TIG/R/Final/crp")## Reading layer `CRP' from data source
## `C:\Users\admin\Documents\Caterina\Maestria TIG\R\Final\crp'
## using driver `ESRI Shapefile'
## Simple feature collection with 3 features and 10 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 5476054 ymin: 6466839 xmax: 5499938 ymax: 6490357
## Projected CRS: POSGAR 2007 / Argentina 5st_crs(crp)## Coordinate Reference System:
## User input: POSGAR 2007 / Argentina 5
## wkt:
## PROJCRS["POSGAR 2007 / Argentina 5",
## BASEGEOGCRS["POSGAR 2007",
## DATUM["Posiciones Geodesicas Argentinas 2007",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]],
## ID["EPSG",1062]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["Degree",0.0174532925199433]]],
## CONVERSION["unnamed",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",-90,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-60,
## ANGLEUNIT["Degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",1,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",5500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1,
## ID["EPSG",9001]]]]crp2<- st_transform(crp, 32621)
st_crs(crp2)## Coordinate Reference System:
## User input: EPSG:32621
## wkt:
## PROJCRS["WGS 84 / UTM zone 21N",
## BASEGEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]],
## CONVERSION["UTM zone 21N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-57,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",0.9996,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",500000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["(E)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["(N)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Engineering survey, topographic mapping."],
## AREA["Between 60°W and 54°W, northern hemisphere between equator and 84°N, onshore and offshore. Barbados. Brazil. Canada - Newfoundland and Labrador, Quebec. French Guiana. Greenland. Guyana. St Pierre and Miquelon. Suriname."],
## BBOX[0,-60,84,-54]],
## ID["EPSG",32621]]9) Se recorta el ráster a partir de una capa vectorial
imagenlandsat2<- imagenlandsat%>%crop(crp2)%>%mask(crp2)
plot(imagenlandsat2)
10) Se realizan diferentes composiciones color
Para poder identificar diversos elementos de interés en la imagen, se ejecutaron diferentes composiciones color.
Color natural
compcolor=plotRGB(imagenlandsat2,
r=4, g=3, b=2,
stretch= 'lin',
margins= TRUE,
main= 'Color natural')
A través de esta composición, que utiliza las bandas roja, verde y azul, se puede visualizar de forma aproximada la imagen en su color natural.
Infrarrojo
compinfra=plotRGB(imagenlandsat2,
r=5, g=4, b=3,
stretch= 'lin',
margins= TRUE,
main= 'Infrarrojo')
Con esta composición, que precisa de las bandas infrarrojo cercano, roja y verde, se puede observar en tonos rojos y rosados la transición de áreas con mucha vegetación a excasa,y en tono negro la presencia de cuerpos de agua.
Uso agrícola
compagricola=plotRGB(imagenlandsat2,
r=6, g=5, b=2,
stretch= 'lin',
margins= TRUE,
main= 'Uso agrícola')
Con esta última composición, en la cual se utilizan las bandas del infrarrojo térmico, infrarrojo cercano y azul, se puede indentificar a través del color verde brillante las zonas de uso de suelo agrícola.
11) Se calculan índices
Índice NDVI
ndvi=(imagenlandsat2[[5]]-imagenlandsat2[[4]])/(imagenlandsat2[[5]]+imagenlandsat2[[4]])
ndvi ## class : RasterLayer
## dimensions : 789, 798, 629622 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 192015, 215955, -3538065, -3514395 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=21 +datum=WGS84 +units=m +no_defs
## source : memory
## names : layer
## values : -0.1070056, 0.5374427 (min, max)plot(ndvi)
Como se puede observar, los valores de píxel de la imagen obtenida se extiendenden desde -0.11 a 0.53. Para poder interpretar de una mejor forma la imagen, se generó una nueva paleta de colores, y se volvió a graficar. Y además, se solicitó el histograma.
colors=colorRampPalette(c('green1', 'gold', 'darkblue'))(255)
plot(ndvi,
xlim=c(192015, 215955),
ylim=c(-3538065, -3514395),
zlim=c(-1, 1),
col=colors,
main='Índice NDVI')
Histograma índice NDVI
hist(ndvi,
main= 'Histograma NDVI',
xlab= 'NDVI',
ylab= 'Frecuencia',
col= 'green',
xlim= c(-1,1),
breaks=10)
Si se observa en conjunto la imagen con el histograma, puede inferirse que al existir mayor cantidad de valores que se ubican entre el 0 y el 0.4, el área de estudio tiende a ser de excasa cobertura vegetal.
Índice SAVI
savi=((imagenlandsat2[[5]]-imagenlandsat2[[4]])/(imagenlandsat2[[5]]+imagenlandsat2[[4]]+0.5))*1.5
savi## class : RasterLayer
## dimensions : 789, 798, 629622 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 192015, 215955, -3538065, -3514395 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=21 +datum=WGS84 +units=m +no_defs
## source : memory
## names : layer
## values : -0.1605039, 0.806153 (min, max)plot(savi,
xlim=c(192015, 215955),
ylim=c(-3538065, -3514395),
zlim=c(-1, 1),
col=colors,
main='SAVI')
Histograma índice SAVI
hist(savi,
main= 'Histograma SAVI',
xlab= 'SAVI',
ylab= 'Frecuencia',
col= 'green',
xlim= c(-1,1),
breaks=10)
En el caso de este índice, los valores se extienden desde -0.16 a 0.80, existiendo mayor presencia de numeros positivos, en comparación al otro índice calculado. El índice SAVI permite reajustar el índice NDVI, obteniendo mejores resultados. Como ya se observó anteriormente, en la imagen puede inferirse que existe un suelo con cierta vegetación, que tiende a ser escasa.
Índice NDWI
ndwi=(imagenlandsat2[[3]]-imagenlandsat2[[5]])/(imagenlandsat2[[3]]+imagenlandsat2[[5]])
ndwi## class : RasterLayer
## dimensions : 789, 798, 629622 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 192015, 215955, -3538065, -3514395 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=21 +datum=WGS84 +units=m +no_defs
## source : memory
## names : layer
## values : -0.5019678, 0.1455864 (min, max)plot(ndwi,
xlim=c(192015, 215955),
ylim=c(-3538065, -3514395),
zlim=c(-1, 1),
col=colors,
main='NDWI')
Histograma índice NDWI
hist(ndwi,
main= 'Histograma NDWI',
xlab= 'NDWI',
ylab= 'Frecuencia',
col= 'green',
xlim= c(-1,1),
breaks=10)
En el caso de este índice, los valores se extienden en el rango de -0.50 a 0.14. Al existir mayor valor de píxeles entre el rango -0.5 a 0, puede inferirse que existe mayor presencia de un suelo seco, sin tanta presencia de agua.
12) Se realiza un proceso de umbralización
A través del proceso de umbralización, se definen valores para los distintos pixeles, con el fin de poder identificar las áreas que sean de interés. En primera medida, considerando el histograma del índice NDVI, donde se indica una mayor existencia para los píxeles de valor 0.2, se optó por eliminar estos, y analizar la nueva imagen generada.
Umbralización utilizando valor escogido de índice NDVI
clases<- reclassify(ndvi, cbind(-Inf, 0.2, NA))
plot(clases,
zlim=c(-1,1),
col=colors,
main='índice NDVI: áreas de mayor vegetación')
Umbralización utilizando valores de índice NDVI de mayor
interés
En este caso, tambien analizando los valores del histograma de NDVI, se eligieron 3 clases diferentes para poder interpretar de una mejor forma la existencia de un suelo totalmente escaso de vegetación, otro intermedio y la presencia de mayor vegetación.
clases2<- reclassify(ndvi, c(-Inf,0.1,1, 0,0.4,2, 0.4,Inf,3))
plot(clases2,
col=rev(terrain.colors(3)),
main='NDVI clases')
Conclusión final
Como se puede observar, en los diferentes análisis que se llevaron a cabo, el software RStudio posee diversas herramientas que permiten la aplicación de distintos análisis sobre un objeto Raster. En el caso de las obras viales, es imprescindible realizar una planificación en base a la caracterización del terreno. A través de la combinación de bandas y el cálculo de cierto índices, se pueden visualizar e identificar rápidamente cuáles son las áreas de poseen mayor presencia de cuerpos de agua, los niveles de vegetación, las zonas mas productivas, entre otros.