-
Notifications
You must be signed in to change notification settings - Fork 4
/
03-tareas_basicas.Rmd
158 lines (108 loc) · 4.88 KB
/
03-tareas_basicas.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
---
title: "Tareas básicas"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
## Empezando a trabajar
Ahora que ya tenemos un set de paquetes que nos puede ser muy útil (existen muchos más relacionados al trabajo con datos espaciales) tenemos que cargar las librerias.
```{r message=FALSE, warning=FALSE, eval=FALSE}
# Cargo las librerías
library(tidyverse)
library(sf) #trabaja con datos vectoriales
library(rgee)
library(mapedit) #trabaja con mapas interactivos
library(raster) #trabaja con datos raster
library(tmap)
```
Y ahora inicializamos GEE, para eso vamos a necesitar nuestro usuario habilitado. La primera vez que lo inicialicemos nos va a solicitar permiso para acceder a GEE por medio de una serie de pantallas en el nevegador de Internet predeterminado. Autorizamos y se nos brindará una API Key para ingresar en la consola de R. Con eso ya estamos autenticados y con acceso a la plataforma.
```{r eval=FALSE}
ee_Initialize('nombre de usuario')
# si queres trabajar con google drive podes inicializar así: ee_Initialize('nombre de usuario', drive = TRUE)
```
GEE tiene colecciones de productos e imágenes. ¿Cómo las consulto?. Lo primero es ir al catálogo de GEE para obtener el nombre de la colección. Yo les dejo los nombres de LANDSAT y SENTINEL, que son dos de las mas utilizadas para el tipo de problema que tenemos que resolver.
* Sentinel 2: COPERNICUS/S2
* LandSat 8: LANDSAT/LC08
Vamos a consultar SENTINEL porque tiene un pixel más pequeño que LANSAT y por ende nos puede brindar mayor cantidad de información:
## Sentinel
```{r eval=FALSE}
# Seleccionando las bandas
bandas <- c('B8A','B4','B11', 'B2', 'B3', 'B5','B6','B7','B8')
# Filtrando los metadatos: usar Abril para mostrar la diferencia
# en la selección de acuerdo al porcentaje de nubes (40, 10, 90).
imagenes_sentinel <- ee$ImageCollection('COPERNICUS/S2')$
select(bandas)$
filterDate('2020-10-01','2020-10-30')$
filterMetadata('CLOUDY_PIXEL_PERCENTAGE','less_than', 40)$
mean()
escala_viz <- list(
bands = c('B8A', 'B11', 'B4'),
min = 0,
max = 10000)
Map$setCenter(-35.662447,-63.783652)
Map$addLayer(imagenes_sentinel, visParams = escala_viz)
```
Ahora bien, no es muy útil procesar todo el mundo cuando nosotros solo necesitamos una región, ahora vamos a ver como recortar un área de estudio y como cambiar algunos parámetros del filtro, como por ejemplo el porcentaje de nubes:
## Cortando el área de estudio
```{r eval=FALSE}
# Definiendo un limite
este_de_la_pampa <- ee$FeatureCollection('users/yabellini/zona_estudio')
# Seleccionando las bandas
bandas <- c('B8A','B4','B11', 'B2', 'B3', 'B5','B6','B7','B8')
# Filtrando los metadatos: usar Abril para mostrar la diferencia
# en la selección de acuerdo al porcentaje de nubes (40, 10, 90).
imagenes_sentinel <- ee$ImageCollection('COPERNICUS/S2')$
select(bandas)$
filterDate('2017-04-01','2017-04-30')$
filterMetadata('CLOUDY_PIXEL_PERCENTAGE','less_than', 40)$
mean()$
clip(este_de_la_pampa)
# Armando una escala de visualización
escala_viz <- list(
bands = c('B8A', 'B11', 'B4'),
min = 0,
max = 10000)
Map$centerObject(este_de_la_pampa, 7)
Map$addLayer(imagenes_sentinel, visParams = escala_viz)
```
## Buscando imágenes
```{r,eval=FALSE}
disponible <- ee$ImageCollection('LANDSAT/LC08/C01/T1_TOA')$
filterDate('2020-04-01','2020-06-30')$
filterBounds(ee$Geometry$Point(-63.783652,-35.662447))
ee_get_date_ic(disponible)
viz = list(min = 0,
max = 0.7,
bands = c('B7','B5','B4'),
gamma = 1.75)
landsat <- ee$Image('LANDSAT/LC08/C01/T1_TOA/LC08_228085_20200428') #Este ID lo saqué del listado anterior
Map$centerObject(eeObject = landsat,zoom = 8)
Map$addLayer(eeObject = landsat,visParams = viz)
```
## Indices multiespectrales
Los indices multiespectrales son diferentes combinaciones de bandas que nos permiten enfocarnos en un tipo de información para analizar una cobertura o un fenómeno en particular. Si no sabemos muy bien que tipo de indices se pueden calcular y para que sirven la herramienta [LandViewer](https://eos.com/landviewer/) es un muy buen lugar para consultar.
## Como calcular NDVI Sentinel 2
```{r eval=FALSE}
# Buscar imagen
coleccion_sen2 <- ee$ImageCollection('COPERNICUS/S2')$
filterDate('2020-10-01','2020-10-30')$
filterBounds(latlon)$
filterMetadata('CLOUDY_PIXEL_PERCENTAGE','less_than',5)
ee_get_date_ic(coleccion_sen2)
# Seleccionar una del listado
id <- 'COPERNICUS/S2/20201026T141049_20201026T142106_T20HMF'
sen2 <- ee$Image(id)
latlon <- ee$Geometry$Point(-63.783652,-35.662447)
Map$centerObject(latlon,zoom = 12)
# Definir paleta de colores
viz <- list(palette = c(
"#d73027", "#f46d43",
"#fdae61", "#fee08b",
"#d9ef8b", "#a6d96a",
"#66bd63", "#1a9850")
)
# Calcular indice
sen2$normalizedDifference(c('B8A','B4')) %>%
Map$addLayer(visParams = viz)
```