PREPROCESAMIENTO DE IMÁGENES DE SATÉLITE: RECORTE DE IMÁGENES
G174: FOTOINTERPRETACIÓN Y TELEDETECCIÓN
Domingo F. Rasilla y Francisco Conde
2024-02-06
1 RECORTE DE IMÁGENES
Los procedimientos que se explicarán en este apartado son útiles para excluir aquellas áreas de las imágenes de satélite que no interesan para su análisis posterior. Con ello, se aligera su “peso”, lo que redunda en una mayor velocidad de procesamiento.
Aunque aparentemente es una acción inocua, conviene tener en cuenta ciertos aspectos que pueden generar problemas. Por ejemplo, en una determinada composición RGB (por ejemplo a color natural o falso color), el recorte afectará a todas las bandas que componen la imagen multibanda; esto influye en la profundidad de píxel de la nueva escena resultante. Si se emplean formatos de archivos (JPG, TIF, PNG, JP2, ECW…) con profundidad de píxel diferente a la del archivo original los valores originales del píxel también variarán, produciendo imágenes con contrastes cromáticos diferentes a los originales.
La resolución espacial de la imagen será otro de los factores a tener en cuenta y determinará un mayor o menor suavizado en los bordes de la imagen en función de la escala. Con imágenes de alta resolución espacial los límites de la imagen recortada serán más límpios que los de imágenes con píxeles de gran tamaño, porque un píxel no puede cortarse por la mitad.
Para recortar, además del fichero ráster que va a ser modificado, es
necesario disponer de otro fichero con los límites de la región de
interés. Este último puede estar en formato ráster
(RasterLayer) o en formato vectorial (objetos
sf), siendo lo último los más frecuente. Para conocer
diferentes formas de obtener ese fichero con los límites de la región de
interés (ROI, “region of interest” en inglés) puede consultarse la
página delimitacion_ROI
Por último, cabe mencionar que el recorte se puede efectuar de dos maneras:
Recorte espacial: es el más habitual, y tiene por objeto reducir el área original de la imagen raster a un zona de tamaño más reducido.
Recorte espectral: cuyo objetivo es trabajar con un número reducido de bandas espectrales, extraídas de una imagen raster multibanda.
2 RECORTE ESPACIAL.
Este tipo de recorte es una labor muy frecuente cuando deseamos trabajar con una área mucho más pequeña que la escena original (el tamaño de una imagen Landsat es de unos 170 km x 185 km, es decir, unoas 31450 \(km^2\)). Con ello se ahorra espacio de almacenamiento y tiempo de procesamiento. Los procedimientos son muy diversos, por lo que el investigador puede elegir el más conveniente para sus intereses.
Como en casos anteriores, vamos a trabajar con varios paquetes de R, que debemos activar.
library(raster)## Loading required package: splibrary(sf)## Linking to GEOS 3.11.1, GDAL 3.6.2, PROJ 9.1.1; sf_use_s2() is TRUE2.1 Recorte espacial de un objeto ráster.
Como en capítulos precedentes, se creará un objeto ráster, ahora formado por 625 píxeles, a las que se incorporarán las coordenadas geográficas.
raster <- raster(ncol= 25, # Número de columnas del ráster
nrow= 25, # Número de filas
xmx= -2,xmn= -5, # Coordenadas del eje X (considerémoslas la longitud)
ymx= 42,ymn= 39) # Coordenadas del eje y (idem latitud)
values(raster) <- as.integer(runif(ncell(raster),
min = 1,
max= 10))
raster## class : RasterLayer
## dimensions : 25, 25, 625 (nrow, ncol, ncell)
## resolution : 0.12, 0.12 (x, y)
## extent : -5, -2, 39, 42 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : layer
## values : 1, 9 (min, max)2.1.1 Ejemplo 1.
En ocasiones, el fichero que corresponde a la región de interés debe ser creado ad hoc. En el siguiente ejemplo, se creará esta región delimitándolo con sus coordenadas geográficas (longitud y latitud, siempre en este orden). En primer lugar, establecerán las coordenadas mediante sendos vectores.
lon <- c(-4.20, -4.20, -3.00,-3.00, -4.20)
lat <- c(40.25,40.64, 40.64, 40.25, 40.25)Para convertir esos vectores en un objeto espacial, deben ser convertidos primero en una lista.
lista <- list(cbind(lon, lat))Este objeto list debe convertirse en un posteriormente
en un objeto sfg con la función
sf:st_polygon.
lista <- st_polygon(lista)
class(lista)## [1] "XY" "POLYGON" "sfg"Los objetos sfg almacenan la geometría con la
información gráfica del objeto: sus coordenadas, sus dimensiones,
geometría etc… Una de las clases típicas es polygon: un
polígono cerrado con cero o más huecos interiores.
Al objeto anterior se le debe añadir el CRS, en este caso el EPSG:4326 (WGS84), que es el que posee el objeto ráster.
limite <- st_sf(st_sfc(lista),
crs = 4326)Por último, y dado que ambos comparten el mismo CRS, se pueden cartografiar de manera simultánea.
plot(raster, axis = TRUE)
plot(st_geometry(limite),
add = TRUE)
A continuación procederemos a enmascarar y recortar el objeto ráster
con el objeto vectorial. La función raster::mask requiere
como argumentos
• El nombre del objeto ráster que será enmascarado y recortado
• El nombre del objeto que servirá de ventana espacial (puede ser vectorial o ráster).
raster_enmascarado <- mask(raster, limite)
plot(raster_enmascarado)
Una vez enmascarado, el objeto ráster es recortado.
raster_recortado <- crop(raster_enmascarado,
limite)Por último, se puede comprobar los resultados representándolos gráficamente:
plot(raster_enmascarado)
plot(st_geometry(limite),
add = TRUE)
La diferencia entre las funciones raster::crop (recorte)
y raster::mask (enmascara ese área de interés) es la
siguiente:
La función
raster::maskmantiene las dimensiones originales del raster, pero los valores situados fuera de los límites del objeto vectorial son convertidos en NA.La función
raster:cropreemplaza esos valores por valores nulos y ajusta la ventana a los límites del objeto vectorial.
Puede observarse que la imagen está recortada, y que los píxeles que se mantienen son aquellos cuya superficie se incluye mayoritariamente dentro de la imagen.
rm(lat, lon, lista, limite, raster_enmascarado, raster_recortado)2.1.2 Ejemplo 2.
Sin embargo, en la mayor parte de los casos se utilizan objetos
vectoriales proporcionados por alguna institución. En este caso, se va a
trabajar con los límites administrativos de la Comunidad de Madrid. El
fichero vectorial con formato .shp se importa con la función
sf::st_rad.
comunidad <- st_read("D:/G174_OCW/Lab_3_preprocesamiento_imagenes_satelite/datos/proyeccion/comunidad.shp")## Reading layer `comunidad' from data source
## `D:\G174_OCW\Lab_3_preprocesamiento_imagenes_satelite\datos\proyeccion\comunidad.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 365635.7 ymin: 4415272 xmax: 495483.9 ymax: 4557310
## Projected CRS: WGS 84 / UTM zone 30NDado que el CRS de este objeto corresponde a un sistema de referencia de coordenadas proyectado (WGS 84 / UTM zone 30N), es más correcto reproyectarlo en vez de reproyectar el objeto ráster.
comunidad <- st_transform(comunidad,
4326) # El CRS correspondiente a coordenadas geográficasA continuación, se representan gráficamente ambos objetos.
plot(raster,
axis = TRUE)
plot(st_geometry(comunidad),
add = TRUE)
Por último, se utilizará el objeto vectorial para recortar y enmascarar el objeto ráster.
raster_enmascarado <- mask(raster,
comunidad)
plot(raster_enmascarado)
raster_recortado <- crop(raster_enmascarado,
comunidad)
plot(raster_recortado)
plot(st_geometry(comunidad),
add = TRUE)
rm(raster, raster_enmascarado, raster_recortado, comunidad)2.2 Recorte espacial de imágenes de satélite
A diferencia de los casos anteriores, en los que se ha trabajado con ficheros unibanda, las escenas de Landsat contienen un conjunto de bandas que deben ser recortadas de la misma manera para su posterior análisis. Este puede hacerse repitiendo, de manera consecutiva y para cada una de las imágenes, los procedimientos anteriores. Esto puede considerarse una pérdida de tiempo, ya que es posible recortar de manera simultánea todas las bandas.
A continuación se muestra un procedimiento para convertir los
ficheros individuales de cada banda en un objeto multibanda utilizando
la función raster::RasterBrick. Esto lo realizaremos en dos
pasos:
- PASO 1: creación de una lista con los ficheros que corresponden a cada una de las bandas. Esta sintaxis ya se conoce del tema anterior.
lista <- list.files("D:/G174_OCW/Lab_2_manejo_ficheros_raster_R/datos", # Directorio de los datos
pattern = ".*B[1234567]\\.tif$", # El signo del dólar lista todos los ficheros que acaban igual
ignore.case=TRUE, # Sensible a mayúsculas y minúsculas
full.names = TRUE)- PASO 2: se importan los ficheros raster de esa lista, inicialmente
como
RasterStack, y luego se transforma enRasterBrick.
stack <- stack(lista)
imagen <- brick(stack)A continuación se consultan las características de este
RasterBrick.
imagen## class : RasterBrick
## dimensions : 7911, 7791, 61634601, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 316485, 550215, 4345785, 4583115 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## source : r_tmp_2024-02-06_080647.781183_76468_10184.grd
## names : LC08_L1TP//9_02_T1_B1, LC08_L1TP//9_02_T1_B2, LC08_L1TP//9_02_T1_B3, LC08_L1TP//9_02_T1_B4, LC08_L1TP//9_02_T1_B5, LC08_L1TP//9_02_T1_B6, LC08_L1TP//9_02_T1_B7
## min values : 8486, 7594, 7003, 5998, 5207, 4755, 4954
## max values : 65535, 57121, 65535, 65535, 65535, 65535, 65535Como ya se ha señalado, el manejo de ficheros con bandas cuya
denominación es tan complicada dificulta cualquier análisis posterior.
Por ello,es conveniente renombrar las bandas para mejorar su
identificación. En este caso, se eliminará el texto
LC08_L1TP_201032_20210723_20210729_02_T1_ de cada imagen,
dejando sólo la letra B (de banda) y el número de cada
banda.
names(imagen) <- gsub(pattern = "LC08_L1TP_201032_20210723_20210729_02_T1_", # conjunto de caracteres que se desea sustituir
replace ="", # caracteres que sustituirán a los anteriores
x <- names(imagen)) # vector o dataframe que sustituirán a los anteriores
imagen## class : RasterBrick
## dimensions : 7911, 7791, 61634601, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 316485, 550215, 4345785, 4583115 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## source : r_tmp_2024-02-06_080647.781183_76468_10184.grd
## names : B1, B2, B3, B4, B5, B6, B7
## min values : 8486, 7594, 7003, 5998, 5207, 4755, 4954
## max values : 65535, 57121, 65535, 65535, 65535, 65535, 65535Por último, con el fin de verificar que todo el proceso se ha
realizado satisfactoriamente, la imagen se representa gráficamente. Como
la imagen es multibanda, se pueden representar de manera simultánea
varias bandas, mediante la función RStoolbox::ggRGB. Esta
función incluye como argumentos:
• imagen -> Un RasterStack o un
RasterBrick.
• r -> capa 1.
• g -> capa 2.
• b -> capa 3.
• stretch -> función para el realce de la imagen. ‘lin’ corresponde a realce lineal.
Para trabajar con la función RStoolbox::ggRGB es
necesario cargar el paquete correspondiente:
library(RStoolbox)La imagen se representa a continuación en una combinación específica de badnas denominada “falso color”.
ggRGB(imagen,
r = 5, g = 4 , b = 3,
stretch = "lin")
El objeto imagen, obtenido a partir de una escena de Lansat, posee un gran tamaño (sus dimensiones son 61634601 píxeles). Para el recorte de esa imagen existen diferentes posibilidades.
2.2.1 Recorte utilizando
la función extent.
Tanto los objetos Raster como los
vectoriales contienen información sobre su extensión
espacial, que adquiere denominaciones diferentes
extensión (
extent) en objetos ráster.cuadro delimitador (
bbox) en objetos vectoriales.
En el caso de los primeros, la extensión se define mediante cuatro valores, que corresponden a las cuatro esquinas del rectángulo de cierre:
las coordenadas mínima y máxima de la dimensión x (
xmin,xmax)las coordenadas minima y maxima de la dimensión y (ymin, ymax),
Si se trabaja simultáneamente con dos objetos ráster, por ejemplo una
escena de Landsat y un MDT o una base de datos de usos del territorio,
la función raster::extent puede ser útil para trasladar los
límites de uno de los objetos (normalmente el más pequeño) al otro. El
procedimiento para recortar una imagen utilizando la función
raster::extent implica los siguientes pasos:
- PASO 1: Las dimensiones originales de la imagen se puede consultar
con
raster::extent(), que devuelve un objetoextent(operativo tanto para objetos ráster como vectoriales).
extent(imagen)## class : Extent
## xmin : 316485
## xmax : 550215
## ymin : 4345785
## ymax : 4583115Los objetos Extent tienen una característica particular:
se pueden cambiar de tamaño fácilmente. Por ejemplo, multiplicar la
extensión por 0,5 no moverá las coordenadas por un 50% hacia la
izquierda y hacia abajo, sino que reducirá la extensión
proporcionalmente, pero manteniendo la ubicación geográfica del
centroide. Un ejemplo es el siguiente:
# borrame <- extent(imagen)
# plot(borrame)
# plot(borrame*0.5, add = TRUE)
# rm(borrame)Este recurso puede ser muy útil en ocasiones, por ejemplo, para una expansión/contracción rápida de las dimensiones de trazado sin tener que jugar con las coordenadas o con la extracción de un subconjunto.
- PASO 2: a continuación se concretarán los límites del área que se
desea recortar. Para ello, se creará un objeto denominado
roia partir de la funciónraster:extent. Este objeto se creará modificando las dimensiones del objetoimagen.
roi <- extent(xmin(imagen) + 50000,
xmax(imagen) - 50000,
ymin(imagen) + 45000,
ymax(imagen) - 83000)
# roi <- extent(400000,500000,4400000,4500000) # Una forma alternativa
roi## class : Extent
## xmin : 366485
## xmax : 500215
## ymin : 4390785
## ymax : 4500115Una vez disponible el objeto roi correspondiente al
nuevo área de trabajo, se usará para recortar la imagen original con la
función raster::crop.
imagen_recortada <- crop(imagen,
roi) Podemos verficar el resultado de este recorte mediante sus metadatos y compararlos con los de la imagen original.
imagen_recortada## class : RasterBrick
## dimensions : 3644, 4457, 16241308, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 366495, 500205, 4390785, 4500105 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## source : r_tmp_2024-02-06_080755.283986_76468_59321.grd
## names : B1, B2, B3, B4, B5, B6, B7
## min values : 9265, 8315, 7134, 6188, 5207, 4755, 4954
## max values : 52017, 51199, 43592, 51099, 57124, 65535, 65535Y también gráficamente mediante una imagen en color natural.
ggRGB(imagen_recortada,
r = 4, g =3, b = 2,
stretch = "lin")
rm(x, lista, stack, imagen_recortada, roi)2.2.2 Recorte a partir de un polígono.
En la mayoría de los casos, sobre todo cuando se trabaja sobre áreas con una delimitación irregular, por ejemplo, un límite administrativo (caso de un ayuntamiento, una comunidad autónoma o un parque nacional) es necesario recortar el objeto ráster inicial apoyándose en un fichero vectorial que incorpora esos límites irregulares.
El fichero vectorial puede añadirse desde el disco duro, en cuyo caso debe estar proyectado en proyectado en el mismo sistema CRS del objeto ráster, o se puede general interactivamente.
2.2.2.1 CASO 1
Si poseemos un fichero vectorial con los límites de la zona de
interés almacenado en nuestro disco duro, podemos importarlo con la
función st_readdel paquete sf.
limites <- st_read("./datos/recorte/originales/poligono.shp")## Reading layer `poligono' from data source
## `D:\G174_OCW\Lab_3_preprocesamiento_imagenes_satelite\datos\recorte\originales\poligono.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 407201.9 ymin: 4434128 xmax: 488607.1 ymax: 4482673
## CRS: NAclass(limites)## [1] "sf" "data.frame"Como se ha señalado, es necesario que tanto el objeto raster como el vectorial tengan la misma proyección. Se puede preguntar a R si ambos ficheros comparten el mismo sistema de proyección.
compareCRS(limites, imagen)## [1] FALSEPara conocer cuáles son los sistemas de proyección de cada formato, se recurrirá a las funciones:
sf::st_crs(para el fichero vectorial).
st_crs(limites) # Se comprueba que no el fichero vectorial no tiene CRS## Coordinate Reference System: NA-raster::crs para el formato ráster
crs(imagen)## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=utm +zone=30 +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 30N",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",-3,
## 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",16030]],
## 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]]]]Como el objeto vectorial carece de proyección, se le asigna un
sistema de proyección con la función sf::st_crs. Existen
varias sintaxis alternativas.
crs_raster <- crs(imagen) # Objeto con el CRS del fichero ráster
#limites2 <- st_set_crs(limites, crs_raster)Una vez establecido el CRS de interés en el objeto vectorial, se procede a su recorte y enmascaramiento.
imagen_recortada <- crop(imagen, limites)
imagen_enmascarada <- mask(imagen_recortada,limites)
plot(imagen_enmascarada)
Y también gráficamente mediante una imagen en la que se combinan las bandas 5, 6 y 4, muy utilizada para la detección de superfiice acuáticas.
ggRGB(imagen_enmascarada,
r = 5, g =6, b = 4,
stretch = "lin")
rm(crs_raster, imagen_enmascarada, imagen_recortada, limites)2.2.2.2 CASO 2
Otra posibilidad es que el fichero vectorial se importe con un sistema de coordenadas de referencia ya establecido. Para ello es necesario descargar el fichero denominado metro que está en formato shapefile.
Primero se importará el fichero vectorial con la función
sf::st_read.
metro <- st_read("D:/G174_OCW/Lab_3_preprocesamiento_imagenes_satelite/datos/recorte/originales/metro.shp")## Reading layer `metro' from data source
## `D:\G174_OCW\Lab_3_preprocesamiento_imagenes_satelite\datos\recorte\originales\metro.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -4.11172 ymin: 40.19909 xmax: -3.410273 ymax: 40.73222
## Geodetic CRS: GCS_unknownA continuación se compararán los sistemas de coordenadas del nuevo
fichero vectorial y del fichero raster con la función
raster::compareCRS`.
compareCRS(metro, imagen)## [1] FALSEA continuación, preguntaremos al sistema cuál es el CRS del fichero vectorial
crs(metro)## [1] "GEOGCRS[\"GCS_unknown\",\n DATUM[\"World Geodetic System 1984\",\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ID[\"EPSG\",6326]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"Degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"longitude\",east,\n ORDER[1],\n ANGLEUNIT[\"Degree\",0.0174532925199433]],\n AXIS[\"latitude\",north,\n ORDER[2],\n ANGLEUNIT[\"Degree\",0.0174532925199433]]]"Debe cambiarse dicho CRS con la función
sf::st_tranform.
metro <- st_transform(metro,
"EPSG:32630")
compareCRS(metro, imagen)## [1] TRUEUna vez transformado el CRS del fichero vectorial, se puede
representar la imagen en formato ráster, y superpuestos, los límites del
área de interés representados por el fichero vectorial. En este caso se
utiliza la función Rstoolbox::plotRGB solicitando que, del
objeto imagen, se representen sólo las bandas 6,5 y 4.
plotRGB(imagen[[6:4]],
stretch = "lin",
axes = TRUE)
plot(metro,
border = "blue",
lwd = 4,
add = TRUE)
Una vez que el raster y el shapefile están en el mismo CRS, se puede crear la máscara.
- PASO 1. Recorte de la imagen (crea un recuadro con los límites
máximo y mínimo de la zona de trabajo) con la función
raster::crop.
crop <- crop(imagen, metro)
ggRGB(crop,
r= 4, g= 3, b= 2,
stretch = "lin")
PASO 2. Recorte de la imagen usando como máscara el poligono con la
función raster::mask. A continuación, su correspondiente
representación gráfica en color natural (RGB 432)
imagen_recortada <- mask(crop,
mask = metro)
ggRGB(imagen_recortada,
r= 4, g= 3, b= 2,
stretch = "lin")
Si se desea grabar el objeto resultante, habría que añadir a la
función raster::maks dos argumentos:
filename = “./datos/recorte/imagen_recortada.tif” Ruta y nombre del fichero resultante.
overwrite = TRUE sobreescribe cualquier fichero anterior con el mism nombre.
Ejemplo (se exporta con formato .tif):
# imagen_recortada <- mask(crop,
# mask = metro,
# filename = "./mi_ruta/nombre_fichero.tif",
# overwrite = TRUE)Si desea agilizar todo este proceso, ambas funciones se pueden encadenar.
borrame <- mask(crop(imagen, metro), metro)
ggRGB(borrame, r= 5, g= 6, b= 2,
stretch = "lin")
2.2.3 Recorte de una imagen manteniendo el área situada fuera del polígono.
Otra posibilidad, no muy frecuente, es el recorte de la zona de
interés, pero vaciando su interior y manteniendo las áreas circundantes.
Para ello, debe incluirse el argumento inverse = TRUE en la
función raster::mask.
borrame <- mask(imagen,
metro,
inverse = TRUE)
ggRGB(borrame,
r= 4, g= 3, b= 2,
stretch = "lin")
2.2.4 Recorte mediante otro fichero raster.
También es posible utilizar como fuente del área de interés otro fichero raster. De nuevo, ambos ficheros ráster deberían tener la misma proyección para que se puedan superponer.
mascara_raster <- brick("./datos/recorte/imagen_recortada.tif")
mascara_raster## class : RasterBrick
## dimensions : 1973, 1984, 3914432, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 405675, 465195, 4450095, 4509285 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## source : imagen_recortada.tif
## names : imagen_recortada_1, imagen_recortada_2, imagen_recortada_3, imagen_recortada_4, imagen_recortada_5, imagen_recortada_6, imagen_recortada_7
## min values : 10076, 9268, 8278, 7456, 7015, 6331, 5932
## max values : 52017, 51199, 40913, 46914, 57124, 65535, 65535Siguiendo el mismo procedimiento realizado en líneas precedentes, primero representaremos toda la imagen y posteriormente supenpodremos la imagen correspondiente al objeto mascara_raster.
plotRGB(imagen[[4:2]],
stretch = "lin",
axes = TRUE,
alpha = 100) # Parámetro que define el nivel de transparencia
plotRGB(mascara_raster[[4:2]],
stretch = "lin",
add = TRUE)
Una vez confirmado el grado de superposición entre ambas imágenes, se puede proceder a crear el subcojunto.
borrame <- crop(imagen,
mascara_raster)
borrame## class : RasterBrick
## dimensions : 1973, 1984, 3914432, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 405675, 465195, 4450095, 4509285 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## source : memory
## names : B1, B2, B3, B4, B5, B6, B7
## min values : 9476, 8547, 7361, 6541, 6129, 5547, 5388
## max values : 52017, 51199, 43592, 51099, 57124, 65535, 65535ggRGB(borrame, r= 5, g= 6, b= 2,
stretch = "lin")
2.2.5 Recorte de un número de píxeles.
Si se quiere crear un subconjunto definiendo el rango de unas celdas,
podemos utilizar la función rasterFromCells. Esta función
sólo puede trabajar con una única capa raster, a diferencia de lo que
ocurre en las opciones anteriores.
borrame <- rasterFromCells(imagen[[2]], # La 2ª banda de la imagen, que sería la green
cells = c(10000:900000)) # Rango de píxeles (número)
borrame # Metadatos de la nueva imagen## class : RasterLayer
## dimensions : 115, 7791, 895965 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 316485, 550215, 4579635, 4583085 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## source : memory
## names : layer
## values : 7792, 903756 (min, max)Este nueva capa se representaría así
plot(borrame)
2.3 Recorte espectral
La selección de alguna banda específica de la imagen de trabajo se
realiza mediante la función subset. Esta(s) nueva(s)
imágen(es) pueden salvarse posteriormente en el disco duro especificando
su nombre en el argumento filename. En el siguiente ejemplo
se seleccionan las bandas correspondientes al verde y al rojo, situadas
en las posiciones 2 y 3.
recorte_espectral <- subset(imagen_recortada,
2:3)
recorte_espectral## class : RasterBrick
## dimensions : 1973, 1984, 3914432, 2 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 405675, 465195, 4450095, 4509285 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## source : memory
## names : B2, B3
## min values : 9268, 8278
## max values : 51199, 40913plot(recorte_espectral)
Se eliminan todos los objetos creados hasta el momento
rm(list=ls())3 MOSAICO DE IMÁGENES
Un mosaico es la combinación de dos o más imágenes raster en una única imagen. A veces, una sola escena no cubre completamente el área de estudio deseada, por lo que puede requerirse combinar dos o más escenas. Las escenas Landsat permiten su superposición espacial, ya que comparten franjas comunes, pero el tamaño de píxel debe ser el mismo para la creación de mosaicos.
Para este ejercicio es necesario descargar los siguientes ficheros
3.1 Importación de los ficheros raster y conversión en objetos ráster
Se crea el objeto ráster mosaico1 a partir de un fichero
en formato .tif
mosaico1 <- brick("./datos/mosaico/originales/mosaico1.tif")A continuación se revisa el contenido de este objeto.
mosaico1 ## class : RasterBrick
## dimensions : 1488, 1142, 1699296, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 414855, 449115, 4438995, 4483635 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## source : mosaico1.tif
## names : mosaico1_1, mosaico1_2, mosaico1_3, mosaico1_4, mosaico1_5, mosaico1_6, mosaico1_7
## min values : 8881, 7609, 6917, 5821, 4546, 4837, 5034
## max values : 65535, 65535, 65535, 65535, 65535, 65535, 65535Es muy habitual que los mosaicos se realicen a partir de dos o más escenas de Landsat contiguas. Es habitual que las escenas incluyan valores NA correspondientes al marco espacial sin datos de cada imagen. En consecuencia, es conveniente asignar el valor 0 a los píxeles con un valor NA de ambas escenas. Este valor sin datos se ignora al asignar las funciones.
NAvalue(mosaico1) <- 0 # Asignación del valor 0 a los píxeles NoDataA continuación se cartografía la primera imagen en falso (FCC), que luega será utilizada en la composición.
plotRGB(mosaico1,
r=3, g=2, b=1,
stretch ='lin',
axes =TRUE,
main = "Primera imagen") 
A continuación se cargan las restantes imágenes.
mosaico2 <- brick("./datos/mosaico/originales/mosaico2.tif")
mosaico3 <- brick("./datos/mosaico/originales/mosaico3.tif")
mosaico4 <- brick("./datos/mosaico/originales/mosaico4.tif")Y se cartografía todas ellas en un único gráfico
par(mfrow = c(2, 2))
plotRGB(mosaico1,
r=3, g=2, b=1,
stretch ='lin',
axes =TRUE,
main = "Primera imagen")
plotRGB(mosaico2,
r=3, g=2, b=1,
stretch ='lin',
axes =TRUE,
main = "Segunda imagen")
plotRGB(mosaico3,
r=3, g=2, b=1,
stretch ='lin',
axes =TRUE,
main = "Tercera imagen")
plotRGB(mosaico4,
r=3, g=2, b=1,
stretch ='lin',
axes =TRUE,
main = "Cuarta imagen")
Volvemos al estado original
par(mfrow = c(1, 1))Cuando el paquete raster crea un mosaico con varias
imágenes, aquella parte que es común (normalmente, en uno o varios
bordes) se superponen entre sí. Además, todos los objetos deben tener un
origen, una resolución y un sistema de coordenadas de referencia
iguales.
Para crear un nuevo objeto raster con unas dimensiones mayores que la
de los objetos individuales debe usarse la función
raster::mosaic. Esta función aplica a su vez una serie de
rutinas para computar el valor de los píxeles en las zonas de
superposición. En R, se pueden aplicar tres funciones a la región
superpuesta: la media, el mínimo y el máximo. Además, la función acepta
como argumento los valores ausentes (na.rm).
3.2 Mosaico aplicando la función promedio a la región superpuesta
Para aplicar la función promedio a la región superpuesta se asigna el
argumento fun = mean. Esta función calcula la media de los
valores de los píxel de la región superpuesta y asigna ese valor medio
al ráster de salida.
mosaico_promedio <- mosaic(mosaico1, mosaico2, mosaico3, mosaico4,
fun = mean)
mosaico_promedio ## class : RasterBrick
## dimensions : 1488, 1142, 1699296, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 414855, 449115, 4438995, 4483635 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## source : memory
## names : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6, layer.7
## min values : 8881, 7609, 6917, 5821, 4546, 4837, 5034
## max values : 65535, 65535, 65535, 65535, 65535, 65535, 65535ggRGB(mosaico_promedio,
r= 3, g= 2, b= 1,
stretch = "lin")
3.3 Mosaico aplicando la
función min() a la región superpuesta
Para aplicar esta función, se asigna el argumento
fun = min, que toma el valor mínimo de los píxeles de las
regiones superpuestas de las imágenes ráster que se van a unir y lo
asigna a la imagen de salida.
mosaico_minimo <- mosaic(mosaico1, mosaico2, mosaico3, mosaico4,
fun = min)
mosaico_minimo ## class : RasterBrick
## dimensions : 1488, 1142, 1699296, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 414855, 449115, 4438995, 4483635 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## source : memory
## names : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6, layer.7
## min values : 8881, 7609, 6917, 5821, 4546, 4837, 5034
## max values : 65535, 65535, 65535, 65535, 65535, 65535, 65535ggRGB(mosaico_minimo, r= 3, g= 2, b= 1,
stretch = "lin")
3.4 Mosaico aplicando la
función max() a la región superpuesta
En este caso se asigana el valor máximo de píxel de las regiones
superpuestas, pasando el argumento fun = max. Crea la mejor
imagen sin pérdida de datos, y por eso es la elegida habitualmente.
mosaico_maximo <- mosaic(mosaico1, mosaico2, mosaico3, mosaico4,
fun = max)
mosaico_maximo## class : RasterBrick
## dimensions : 1488, 1142, 1699296, 7 (nrow, ncol, ncell, nlayers)
## resolution : 30, 30 (x, y)
## extent : 414855, 449115, 4438995, 4483635 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs
## source : memory
## names : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6, layer.7
## min values : 8881, 7609, 6917, 5821, 4546, 4837, 5034
## max values : 65535, 65535, 65535, 65535, 65535, 65535, 65535ggRGB(mosaico_maximo,
r= 3, g= 2, b= 1,
stretch = "lin")
Como estos objetos ya no van a ser de utilidad, se eliminan del Global Enviroment.
rm(mosaico1, mosaico2, mosaico3,mosaico4, mosaico_promedio, mosaico_minimo, mosaico_maximo)