Analisis de archivos raster
VE-CC
DataIntelligenceMiercoles 12-05-2021
Abstract
El siguiente informe está enfocado al estudio y manipulación de rasterslectura de raster
librerias “raster”
Leemos el raster
lulc <- raster("1/lulc.tif")
lulc## class : RasterLayer
## dimensions : 208, 178, 37024 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 359560, 361340, 6314590, 6316670 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs
## source : C:/Users/usuario/Desktop/ds/ds_rgee/estudio/1/lulc.tif
## names : lulc
## values : 1, 3 (min, max)plot(lulc)
Cambiar coordenadas de raster
librerias “raster”
Leemos el raster
lulc <- raster("2/lulc.tif")
lulc## class : RasterLayer
## dimensions : 208, 178, 37024 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 359560, 361340, 6314590, 6316670 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs
## source : C:/Users/usuario/Desktop/ds/ds_rgee/estudio/2/lulc.tif
## names : lulc
## values : 1, 3 (min, max)Para cambiar las coordenadas utilizamos la función projectRaster
luc <- lulc %>% projectRaster(crs = "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs ")
luc## class : RasterLayer
## dimensions : 222, 192, 42624 (nrow, ncol, ncell)
## resolution : 0.000107, 9.02e-05 (x, y)
## extent : -70.50895, -70.48841, -33.29937, -33.27934 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : lulc
## values : 1, 3 (min, max)Calculos de raster
librerias “raster”
Se pueden efectuar diversas operaciones con los raster:
dem <- raster("3/dem.tif")
dem## class : RasterLayer
## dimensions : 67, 68, 4556 (nrow, ncol, ncell)
## resolution : 0.0002694946, 0.0002694946 (x, y)
## extent : -70.50799, -70.48967, -33.29835, -33.28029 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : C:/Users/usuario/Desktop/ds/ds_rgee/estudio/3/dem.tif
## names : dem
## values : 1282, 1906 (min, max)Divicion:
dem2 <- dem/2
dem2## class : RasterLayer
## dimensions : 67, 68, 4556 (nrow, ncol, ncell)
## resolution : 0.0002694946, 0.0002694946 (x, y)
## extent : -70.50799, -70.48967, -33.29835, -33.28029 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : dem
## values : 641, 953 (min, max)Aplicación de promedios:
dem %>% values %>% mean## [1] 1573.888Despliegues detallados:
dem %>% values %>% summary## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1282 1483 1578 1574 1659 1906vector a raster
librerias “raster”, “sf”
Leemos el objeto vectorial
parches <- read_sf("4/parches.shp")
plot(parches)
Creamos un raster ficticio
r <-raster(resolution = c(30,30),
crs = crs(parches), # <- se crea raster vacio
ext = extent(parches))
r## class : RasterLayer
## dimensions : 384, 741, 284544 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 352937.7, 375167.7, 8996955, 9008475 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=24 +south +datum=WGS84 +units=m +no_defsPara hacer la conversión utilizamos la función rasterize y el ráster ficticio que hicimos anteriormente
par_ras <- rasterize(parches, r, lield =parches$area)
par_ras## class : RasterLayer
## dimensions : 384, 741, 284544 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 352937.7, 375167.7, 8996955, 9008475 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=24 +south +datum=WGS84 +units=m +no_defs
## source : memory
## names : layer
## values : 1, 6 (min, max)
## attributes :
## ID area
## from: 1 3080.22561
## to : 6 45.43273plot(par_ras)
para guardar el raster:
# writeRaster(par_ras, "parches2.tif")raster a vectorial
librerias “raster”, “sf”, “sp”
Leemos el raster
lulc <- raster("5/lulc.tif")Utilizamos la función rasterToPolygons para colocar el ráster luego la función st_set_sf para convertirlo en archivo vectorial
al hacer esto cada pixel se convirtió en un vector, pero eso no es lo que se busca
shp <- rasterToPolygons(lulc) %>% st_as_sf #<- cada pixel es un vector
plot(shp)
Para hacer la conversión de manera correcta colocamos el parámetro dissolve
Esto no da 3 polígonos:
shp <- rasterToPolygons(lulc, dissolve = T) %>% st_as_sf #<- contiene 3 poligonos uno por color ## Loading required namespace: rgeosplot(shp)
Para separarlo utilizamos la función disaggregate para tener varios polígonos independientes
shp <- rasterToPolygons(lulc, dissolve = T) %>% disaggregate %>% st_as_sf #<- separa poligonos
plot(shp)
Reclasificar Raster
librerias “raster”, “sf”, “sp”
Valores Discretos
Leemos el raster
lulc <- raster("6/lulc.tif")
plot(lulc)
Creamos una matriz, en la primera columna pondremos los valores que posee el ráster y en la segunda pondremos los valores que queremos asignar
cambio <- cbind(c(1,2,3),
c(1,2,2))
cambio## [,1] [,2]
## [1,] 1 1
## [2,] 2 2
## [3,] 3 2Para efectuar el cambio se utiliza la función reclassify que junto con el parámetro rcl parece tener una función similar al merge con la asignación de valores
reclass <- reclassify(lulc, rcl = cambio)
plot(reclass)
Valores Continuos
Leemos el raster
ndvi <- raster("6/ndvi.tif")
plot(ndvi)
Para poder reclasificar el NDVI necesitamos extraer sus valores
ndvi %>% values %>% summary## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.0410 0.2850 0.3623 0.3638 0.4310 0.8282 607Se crea un vector con los valores que ocupemos
m <- c(0, 0.2850, 1, #<- valor minimo(0) hasta 1er quartil
0.2850, 0.3623, 2,#<- del primer cuartil hasta mediana
0.3623, 0.4310, 3,#<- de mediana a 3er quartil
0.4310, 1, 4)#<- 3se quartil hasta valor maximo(1)
m## [1] 0.0000 0.2850 1.0000 0.2850 0.3623 2.0000 0.3623 0.4310 3.0000 0.4310
## [11] 1.0000 4.0000Convertimos el vector a matris
mat <- matrix(m, ncol = 3, byrow = T)
mat## [,1] [,2] [,3]
## [1,] 0.0000 0.2850 1
## [2,] 0.2850 0.3623 2
## [3,] 0.3623 0.4310 3
## [4,] 0.4310 1.0000 4Utilizamos nuevamente la función reclassify para hacer la reclasificación
reclass <- reclassify(ndvi, mat)
plot(reclass)
Cambiar tamnio de pixel(Resampling)
librerias “raster”, “sf”, “sp”
Leemos el raster
lulc <- raster("7/lulc.tif")
lulc## class : RasterLayer
## dimensions : 208, 178, 37024 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 359560, 361340, 6314590, 6316670 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs
## source : C:/Users/usuario/Desktop/ds/ds_rgee/estudio/7/lulc.tif
## names : lulc
## values : 1, 3 (min, max)plot(lulc) Creamos un ráster vacío con el tamaño de pixeles que queramos
r <- raster(resolution = c(30,30),
crs = crs(lulc),
ext = extent(lulc))
r## class : RasterLayer
## dimensions : 69, 59, 4071 (nrow, ncol, ncell)
## resolution : 30, 30 (x, y)
## extent : 359560, 361330, 6314600, 6316670 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defsCon la función resample y el method ngb hacemos la asignación
ras_lulc <- resample(lulc, r, method = "ngb")
plot(ras_lulc)
Cortar raster
librerias “raster”, “sf”, “sp”
Leemos el raster
ndvi <- raster("8/ndvi.tif")
ndvi## class : RasterLayer
## dimensions : 198, 168, 33264 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 359610, 361290, 6314640, 6316620 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs
## source : C:/Users/usuario/Desktop/ds/ds_rgee/estudio/8/ndvi.tif
## names : ndviplot(ndvi) Leemos el archivo vectorial
area <- read_sf("8/area.shp")Para poder hacer el corte vamos a utilizar la función crop acompañada de la función mask, esto debido a que en algunos casos cunado el archivo vectorial es muy grande tiende a hacer el corte mal, lo que se puede evitar utilizando mask
ndvi2 <- ndvi %>% crop(area) %>% mask(area)
plot(ndvi2)
Carcular NDVI
librerias “raster”, “sf”, “sp”
Leemos el shp con area de estudio
area_est <- read_sf("9/area_est.shp")Leemos los raster y cortamos posteriormente
Banda
red <- raster("9/B04_10m.jp2") %>%
crop(area_est) %>% mask(area_est)
plot(red)
Banda 8
nir <- raster("9/B08_10m.jp2") %>%
crop(area_est) %>% mask(area_est)
plot(red)
Calculamos el NDVI
ndvi <- ((nir/10000 - red/10000)/(nir/10000 + red/10000))
plot(ndvi)
Extraer informacion de raster
librerias “raster”, “sf”, “sp”
Leemos los archivos necesarios
Es importante saber que para poder extraer información, todos los archivos deben tener el mismo formato de coordenadas.
En nuestro caso el raster dem tiene un formato diferente al resto por lo que tendremos que modificarlo, deben estar en 19 Sur
ndvi <- raster("10/ndvi.tif")
ndvi## class : RasterLayer
## dimensions : 198, 168, 33264 (nrow, ncol, ncell)
## resolution : 10, 10 (x, y)
## extent : 359610, 361290, 6314640, 6316620 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs
## source : C:/Users/usuario/Desktop/ds/ds_rgee/estudio/10/ndvi.tif
## names : ndvidem <- raster("10/dem.tif")
dem <- dem %>% projectRaster(crs = "+proj=utm +zone=19 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
dem## class : RasterLayer
## dimensions : 78, 81, 6318 (nrow, ncol, ncell)
## resolution : 25.1, 29.9 (x, y)
## extent : 359419.6, 361452.7, 6314467, 6316799 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=19 +south +datum=WGS84 +units=m +no_defs
## source : memory
## names : dem
## values : 1284.318, 1905.975 (min, max)puntos <- read_sf("10/ptos_ext.shp")
puntos## Simple feature collection with 25 features and 1 field
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 359648 ymin: 6314753 xmax: 361194.9 ymax: 6316476
## Projected CRS: WGS 84 / UTM zone 19S## Warning: `...` is not empty.
##
## We detected these problematic arguments:
## * `needs_dots`
##
## These dots only exist to allow future extensions and should be empty.
## Did you misspecify an argument?## # A tibble: 25 x 2
## id geometry
## <chr> <POINT [m]>
## 1 1a (360741.6 6315008)
## 2 1b (360094.1 6316379)
## 3 1c (359777.5 6314972)
## 4 2a (360086.9 6315835)
## 5 2b (361194.9 6316476)
## 6 3a (361166.1 6314807)
## 7 4a (360533 6315857)
## 8 4b (360040.1 6315310)
## 9 4c (359648 6316238)
## 10 4d (360586.9 6316321)
## # ... with 15 more rowsExtraemos los puntos de nuestros raster con la función extract especificando que es del paquete raster
puntos_ndvi <- ndvi %>% raster::extract(puntos)
puntos_dem <- dem %>% raster::extract(puntos)Extracción de nombres de los puntos
ptos_info <- tibble(id = puntos$id,
ndvi = puntos_ndvi,
altitud = puntos_dem)
ptos_info## Warning: `...` is not empty.
##
## We detected these problematic arguments:
## * `needs_dots`
##
## These dots only exist to allow future extensions and should be empty.
## Did you misspecify an argument?## # A tibble: 25 x 3
## id ndvi altitud
## <chr> <dbl> <dbl>
## 1 1a 0.734 1329.
## 2 1b 0.109 1592.
## 3 1c 0.375 1598.
## 4 2a 0.193 1699.
## 5 2b 0.351 1861.
## 6 3a 0.386 1452.
## 7 4a 0.275 1647.
## 8 4b 0.336 1602.
## 9 4c 0.209 1506.
## 10 4d 0.244 1750.
## # ... with 15 more rowsExtraemos el primer registro sin el id
cor(ptos_info[,-1])## ndvi altitud
## ndvi 1.0000000 -0.6818164
## altitud -0.6818164 1.0000000Del siguiente grafico podemos deducir que mayor altura, el índice NDVI es mas bajo
ptos_info %>% ggplot(aes(x = altitud, y = ndvi)) +
geom_line()+
stat_smooth(method = "lm", se = F)## `geom_smooth()` using formula 'y ~ x'