-
Notifications
You must be signed in to change notification settings - Fork 2
/
04-rgee.Rmd
169 lines (113 loc) · 5.37 KB
/
04-rgee.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
159
160
161
162
163
164
165
166
167
168
169
---
title: "Datos rasters y rgee"
output:
html_document:
code_download: true
highlight: tango
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
## Datos raster
Las imágenes satelitales son en general datos que se presentan en formto raster. Google Earth Engine es una plataforma que nos permite acceder a muchos datos históricos de diversos satélites por medio de su API (application programming interface, o interfaz de programación de aplicaciones).
Para poder utilizarla programando es necesario que contemos con un usuario y tengamos instalado todo el software detallado en la página del curso.
## rgee y Google Earth Engine
**rgee** es un paquete que hace de interfaz entre R y Google Earth Engine. Como hemos hecho con los otros paquetes, para poder usarlo tenemos que cargarlo:
```{r}
library(reticulate)
library(rgee)
```
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:
* Sentinel 2: COPERNICUS/S2
* LandSat 8: LANDSAT/LC08
Vamos a consultar SENTINEL porque tiene un pixel más pequeño que LANDSAT 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)
```
![Mapa del mundo con todas las imágenes sentinel](img/sentinel_completo.png)
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)
```
![Recorte del este de la provincia de La Pampa](img/recorte.png)
## 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)
```
![Ejemplo de escena LANDSAT](img/imagen_landsat.png)
## 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
latlon <- ee$Geometry$Point(-63.783652,-35.662447)
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)
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)
```
![NDVI con Sentinel](img/NDVI_sentinel2.png)