R y GEE

El paquete rgee

VE-MS-CC

DataIntelligence

12-04-2021

Abstract

Éste informe se basa en la página https://csaybar.github.io/rgee-examples/ que contiene un poco más de 250 ejemplos de códigos rgee. La intención es explorar la información contenida en las imágenees satelitales (que son simplemente bases de datos) y por supuesto, explorar las potencialidades gráficas.

El trabajo se divide en 5 partes. La primera entrega una introducción a las características básicas, la segunda, ejercicios sobre algunos modelos de aprendizaje automático (CART, SVM y KMeans) para pasar luego a revisar el tema de imágenes. La cuarta parte aborda el tema de las colecciones (ImageCollection) y una quinta, Geometry, Feature y FeatureCollection.

Catalogo global Earth Engine:

https://developers.google.com/earth-engine/datasets/catalog

Índice

1. La información contenida en las imágenes

2. Aprendizaje automático

3. Imágenes

4. ImageCollection

5. Geometry, Feature y FeatureCollection

1. La información contenida en las imágenes

Gee opera a nivel de colleciones de imágenes (escenas) llamadas ImageCollection obtenidas de satélites.

Analicemos por ejemplo la siguiente colección:
USGS Landsat 8 Collection 1 Tier 1 Raw Scenes,

cuyo Earth Engine Snippet (que opera a modo de ID) es:
ee.ImageCollection(“LANDSAT/LC08/C01/T1”)

LANDSAT/LC08/C01/T1 significa es el ID del satélite Landsat 8 lanzado por el servicio geológico de los EE.UU. (USGS). Completa su órbita de 705 km de altura cada 99 minutos, y revisita un mismo punto sobre la superficie de la Tierra cada 16 días con un desfase de 8 días con respecto al satélite Landsat 7, del mismo proyecto. Bajo estas condiciones el satélite adquiere cerca de 650 imágenes diariamente.

1.1 Los niveles

El postfijo T1 hace referencia a los niveles.

El servicio geológico de los EE.UU (USGS) produce datos en 3 niveles (tiers) para cada satélite:

Nivel 1 (T1): datos que cumplen con los requisitos de calidad geométricos y radiométricos.
Nivel 2 (T2): datos que no cumplen con los requisitos del Nivel 1.
Tiempo real (RT): datos que aún no se han evaluado (se necesitan hasta un mes).

LANDSAT/LC08/C01/T1 Landsat 8, Colección 1, Nivel 1 únicamente
LANDSAT/LC08/C01/T2 Landsat 8, Colección 1, Nivel 2 únicamente.
LANDSAT/LC08/C01/T1_RT Landsat 8, Colección 1, Nivel 1 + Tiempo real

Los valores números de datos enteros (DNs) de Tier 1 de Landsat 8 Collection 1, representan la radiancia en el sensor calibrada y escalada.

Las escenas Landsat con la mayor calidad de datos disponible se colocan en el Nivel 1 y se consideran adecuadas para el análisis de procesamiento de series de tiempo. El nivel 1 incluye datos procesados del terreno de precisión de nivel 1 (L1TP) que tienen una radiometría bien caracterizada y están intercalibrados en los diferentes sensores Landsat. El georegistro de escenas de Nivel 1 será consistente y dentro de las tolerancias prescritas [<= 12 m de error cuadrático medio (RMSE)]. Todos los datos de Tier 1 Landsat pueden considerarse consistentes e intercalibrados (independientemente del sensor) en toda la colección.

1.2 Las bandas

En astronomía, una banda espectral designa una parte del espectro electromagnético que deja pasar un filtro estándar. Una banda espectral está determinada por su perfil de transmisión, es decir, la fracción de intensidad luminosa transmitida por una longitud de onda dada.

El Landsat 8 transporta dos instrumentos:

el generador de imágenes de tierra operativa (OLI), provee acceso a nueve bandas espectrales que cubren el espectro desde los 0.433 μm a los 1.390 μm y
el sensor infrarrojo térmico (TIRS), que registra de 10.30μm a 12.50μm.

Nombre	tamaño del pixel	Longitud de onda	Descripción
B1	30 metros	0.43 - 0.45 µm	Aerosol costero
B2	30 metros	0.45 - 0.51 µm	Azul
B3	30 metros	0.53 - 0.59 µm	Verde
B4	30 metros	0.64 - 0.67 µm	rojo
B5	30 metros	0.85 - 0.88 µm	Infrarrojo cercano
B6	30 metros	1.57 - 1.65 µm	Infrarrojos de onda corta 1
B7	30 metros	2.11 - 2.29 µm	Infrarrojos de onda corta 2
B8	15 metros	0.52 - 0.90 µm	Banda 8 Pancromática
B9	15 metros	1.36 - 1.38 µm	Cirro
B10	30 metros	10.60 - 11.19 µm	Infrarrojo térmico 1, muestreado de 100 ma 30 m
B11	30 metros	11.50 - 12.51 µm	Infrarrojo térmico 2, muestreado de 100 ma 30 m
BQA			Máscara de bits de control de calidad de Landsat Collection 1 Landsat)

1.3 Los datos asociados a una colección de imágenes filtrada: información general

Necesitamos el conjunto de escenas landsat que cubran Chile

El Path y el Row de Landsat. Son parámetros numéricos que permiten identificar una imagen satelital Landsat de forma análoga a los valores de longitud y latitud. Landsat mapea de manera constante las mismas superficies a través de una malla de cuadrículas. Cada cuadrícula presenta de manera constante el mismo Path y Row, por lo que, buscando por estos parámetros, se consiguen todas las imágenes históricas de Landsat para el mismo lugar.

Ahora, cuales son el Path y el Row para Santiago? Encontramos un convertidor automatico en linea: https://landsat.usgs.gov/landsat_acq#convertPathRow

# Realiza una llamada a alguna de las colecciones de Landsat
object <- ee$ImageCollection("LANDSAT/LC08/C01/T1")$
  filter(ee$Filter()$eq("WRS_PATH", 233))$
  filter(ee$Filter()$eq("WRS_ROW", 83))$
  filterMetadata ('CLOUD_COVER', 'Less_Than', 20)$
  filterDate("2020-03-01", "2020-04-01")$aside(ee_print)

## ------------------------------------------------ Earth Engine ImageCollection --
## ImageCollection Metadata:
##  - Class                      : ee$ImageCollection
##  - Number of Images           : 2
##  - Number of Properties       : 38
##  - Number of Pixels*          : 1415544744
##  - Approximate size*          : 28.86 GB
## Image Metadata (img_index = 0):
##  - ID                         : LANDSAT/LC08/C01/T1/LC08_233083_20200314
##  - system:time_start          : 2020-03-14 14:33:32
##  - system:time_end            : 2020-03-14 14:33:32
##  - Number of Bands            : 12
##  - Bands names                : B1 B2 B3 B4 B5 B6 B7 B8 B9 B10 B11 BQA
##  - Number of Properties       : 119
##  - Number of Pixels*          : 707772372
##  - Approximate size*          : 14.43 GB
## Band Metadata (img_band = 'B1'):
##  - EPSG (SRID)                : WGS 84 / UTM zone 19N (EPSG:32619)
##  - proj4string                : +proj=utm +zone=19 +datum=WGS84 +units=m +no_defs
##  - Geotransform               : 30 0 227985 0 -30 -3556485
##  - Nominal scale (meters)     : 30
##  - Dimensions                 : 7643 7717
##  - Number of Pixels           : 58981031
##  - Data type                  : INT
##  - Approximate size           : 1.20 GB
##  --------------------------------------------------------------------------------
##  NOTE: (*) Properties calculated considering a constant  geotransform and data type.

1.3 Los datos asociados a una colección de imágenes filtrada: información de las bandas.

La tabla rescatada de una imagen posee una serie de 11 conjuntos de datos de 7 columnas cada una que contiene información relativa a las bandas. Son 11 bandas.

first <- object$first()
a_001 <- first$getInfo()
a_001 <- as.data.frame(a_001)
a_001 %>%  kbl() %>%
 kable_material(c("striped", "hover"), font_size = 12)%>%
  scroll_box(width = "100%", height = "200px")

type	bands.id	bands.data_type.type	bands.data_type.precision	bands.data_type.max	bands.dimensions	bands.crs	bands.crs_transform	bands.id.1	bands.data_type.type.1	bands.data_type.precision.1	bands.data_type.max.1	bands.dimensions.1	bands.crs.1	bands.crs_transform.1	bands.id.2	bands.data_type.type.2	bands.data_type.precision.2	bands.data_type.max.2	bands.dimensions.2	bands.crs.2	bands.crs_transform.2	bands.id.3	bands.data_type.type.3	bands.data_type.precision.3	bands.data_type.max.3	bands.dimensions.3	bands.crs.3	bands.crs_transform.3	bands.id.4	bands.data_type.type.4	bands.data_type.precision.4	bands.data_type.max.4	bands.dimensions.4	bands.crs.4	bands.crs_transform.4	bands.id.5	bands.data_type.type.5	bands.data_type.precision.5	bands.data_type.max.5	bands.dimensions.5	bands.crs.5	bands.crs_transform.5	bands.id.6	bands.data_type.type.6	bands.data_type.precision.6	bands.data_type.max.6	bands.dimensions.6	bands.crs.6	bands.crs_transform.6	bands.id.7	bands.data_type.type.7	bands.data_type.precision.7	bands.data_type.max.7	bands.dimensions.7	bands.crs.7	bands.crs_transform.15L	bands.crs_transform.227992.5	bands.crs_transform..15L	bands.crs_transform..3556492.5	bands.id.8	bands.data_type.type.8	bands.data_type.precision.8	bands.data_type.max.8	bands.dimensions.8	bands.crs.8	bands.crs_transform.7	bands.id.9	bands.data_type.type.9	bands.data_type.precision.9	bands.data_type.max.9	bands.dimensions.9	bands.crs.9	bands.crs_transform.8	bands.id.10	bands.data_type.type.10	bands.data_type.precision.10	bands.data_type.max.10	bands.dimensions.10	bands.crs.10	bands.crs_transform.9	bands.id.11	bands.data_type.type.11	bands.data_type.precision.11	bands.data_type.max.11	bands.dimensions.11	bands.crs.11	bands.crs_transform.10	id	version	properties.RADIANCE_MULT_BAND_5	properties.RADIANCE_MULT_BAND_6	properties.RADIANCE_MULT_BAND_3	properties.RADIANCE_MULT_BAND_4	properties.RADIANCE_MULT_BAND_1	properties.RADIANCE_MULT_BAND_2	properties.K2_CONSTANT_BAND_11	properties.K2_CONSTANT_BAND_10	properties.system.footprint.type	properties.system.footprint.coordinates.c..71.6959135149003…33.0059378964372.	properties.system.footprint.coordinates.c..71.7318863847827…33.1312388061595.	properties.system.footprint.coordinates.c..71.7908117374045…33.3359154182762.	properties.system.footprint.coordinates.c..71.8683582044524…33.6041280077762.	properties.system.footprint.coordinates.c..71.9124777716106…33.7560471367911.	properties.system.footprint.coordinates.c..71.9339349062747…33.829775367711.	properties.system.footprint.coordinates.c..71.9356987447135…33.8361791821452.	properties.system.footprint.coordinates.c..71.5455053426038…33.9148443298207.	properties.system.footprint.coordinates.c..70.6854706429476…34.0834168578594.	properties.system.footprint.coordinates.c..69.9052530178163…34.2306187877162.	properties.system.footprint.coordinates.c..69.8959384104283…34.2019950853064.	properties.system.footprint.coordinates.c..69.6856899441946…33.4123402109338.	properties.system.footprint.coordinates.c..69.4793373747696…32.6221900876053.	properties.system.footprint.coordinates.c..69.4612194171066…32.5500004218985.	properties.system.footprint.coordinates.c..69.4513225578524…32.5102001573257.	properties.system.footprint.coordinates.c..71.4088930930099…32.1293142151537.	properties.system.footprint.coordinates.c..71.4443396714738…32.1221360647958.	properties.system.footprint.coordinates.c..71.4841366382819…32.2621766095897.	properties.system.footprint.coordinates.c..71.5141106873013…32.367989956654.	properties.system.footprint.coordinates.c..71.5790943258888…32.596968452028.	properties.system.footprint.coordinates.c..71.6499735469383…32.8454094439389.	properties.system.footprint.coordinates.c..71.6959135149003…33.0059378964372..1	properties.REFLECTIVE_SAMPLES	properties.SUN_AZIMUTH	properties.CPF_NAME	properties.DATE_ACQUIRED	properties.ELLIPSOID	properties.STATION_ID	properties.RESAMPLING_OPTION	properties.ORIENTATION	properties.WRS_ROW	properties.RADIANCE_MULT_BAND_9	properties.TARGET_WRS_ROW	properties.RADIANCE_MULT_BAND_7	properties.RADIANCE_MULT_BAND_8	properties.IMAGE_QUALITY_TIRS	properties.TRUNCATION_OLI	properties.CLOUD_COVER	properties.GEOMETRIC_RMSE_VERIFY	properties.COLLECTION_CATEGORY	properties.GRID_CELL_SIZE_REFLECTIVE	properties.CLOUD_COVER_LAND	properties.GEOMETRIC_RMSE_MODEL	properties.COLLECTION_NUMBER	properties.IMAGE_QUALITY_OLI	properties.LANDSAT_SCENE_ID	properties.WRS_PATH	properties.PANCHROMATIC_SAMPLES	properties.PANCHROMATIC_LINES	properties.GEOMETRIC_RMSE_MODEL_Y	properties.REFLECTIVE_LINES	properties.TIRS_STRAY_LIGHT_CORRECTION_SOURCE	properties.GEOMETRIC_RMSE_MODEL_X	properties.system.asset_size	properties.system.index	properties.REFLECTANCE_ADD_BAND_1	properties.REFLECTANCE_ADD_BAND_2	properties.DATUM	properties.REFLECTANCE_ADD_BAND_3	properties.REFLECTANCE_ADD_BAND_4	properties.RLUT_FILE_NAME	properties.REFLECTANCE_ADD_BAND_5	properties.REFLECTANCE_ADD_BAND_6	properties.REFLECTANCE_ADD_BAND_7	properties.REFLECTANCE_ADD_BAND_8	properties.BPF_NAME_TIRS	properties.GROUND_CONTROL_POINTS_VERSION	properties.DATA_TYPE	properties.UTM_ZONE	properties.system.time_end	properties.LANDSAT_PRODUCT_ID	properties.REFLECTANCE_ADD_BAND_9	properties.GRID_CELL_SIZE_PANCHROMATIC	properties.RADIANCE_ADD_BAND_4	properties.REFLECTANCE_MULT_BAND_7	properties.system.time_start	properties.RADIANCE_ADD_BAND_5	properties.REFLECTANCE_MULT_BAND_6	properties.RADIANCE_ADD_BAND_6	properties.REFLECTANCE_MULT_BAND_9	properties.PROCESSING_SOFTWARE_VERSION	properties.RADIANCE_ADD_BAND_7	properties.REFLECTANCE_MULT_BAND_8	properties.RADIANCE_ADD_BAND_1	properties.RADIANCE_ADD_BAND_2	properties.RADIANCE_ADD_BAND_3	properties.REFLECTANCE_MULT_BAND_1	properties.RADIANCE_ADD_BAND_8	properties.REFLECTANCE_MULT_BAND_3	properties.RADIANCE_ADD_BAND_9	properties.REFLECTANCE_MULT_BAND_2	properties.REFLECTANCE_MULT_BAND_5	properties.REFLECTANCE_MULT_BAND_4	properties.THERMAL_LINES	properties.TIRS_SSM_POSITION_STATUS	properties.GRID_CELL_SIZE_THERMAL	properties.NADIR_OFFNADIR	properties.RADIANCE_ADD_BAND_11	properties.REQUEST_ID	properties.EARTH_SUN_DISTANCE	properties.TIRS_SSM_MODEL	properties.FILE_DATE	properties.SCENE_CENTER_TIME	properties.SUN_ELEVATION	properties.BPF_NAME_OLI	properties.RADIANCE_ADD_BAND_10	properties.ROLL_ANGLE	properties.K1_CONSTANT_BAND_10	properties.SATURATION_BAND_1	properties.SATURATION_BAND_2	properties.SATURATION_BAND_3	properties.SATURATION_BAND_4	properties.SATURATION_BAND_5	properties.MAP_PROJECTION	properties.SATURATION_BAND_6	properties.SENSOR_ID	properties.SATURATION_BAND_7	properties.K1_CONSTANT_BAND_11	properties.SATURATION_BAND_8	properties.SATURATION_BAND_9	properties.TARGET_WRS_PATH	properties.RADIANCE_MULT_BAND_11	properties.RADIANCE_MULT_BAND_10	properties.GROUND_CONTROL_POINTS_MODEL	properties.SPACECRAFT_ID	properties.ELEVATION_SOURCE	properties.THERMAL_SAMPLES	properties.GROUND_CONTROL_POINTS_VERIFY
Image	B1	PixelType	int	65535	7671	EPSG:32619	30	B2	PixelType	int	65535	7671	EPSG:32619	30	B3	PixelType	int	65535	7671	EPSG:32619	30	B4	PixelType	int	65535	7671	EPSG:32619	30	B5	PixelType	int	65535	7671	EPSG:32619	30	B6	PixelType	int	65535	7671	EPSG:32619	30	B7	PixelType	int	65535	7671	EPSG:32619	30	B8	PixelType	int	65535	15341	EPSG:32619	15	227992.5	-15	-3556493	B9	PixelType	int	65535	7671	EPSG:32619	30	B10	PixelType	int	65535	7671	EPSG:32619	30	B11	PixelType	int	65535	7671	EPSG:32619	30	BQA	PixelType	int	65535	7671	EPSG:32619	30	LANDSAT/LC08/C01/T1/LC08_233083_20200314	-1	0.0061836	0.0015378	0.011983	0.010105	0.012699	0.013004	1201.144	1321.079	LinearRing	-71.69591	-71.73189	-71.79081	-71.86836	-71.91248	-71.93393	-71.93570	-71.54551	-70.68547	-69.90525	-69.89594	-69.68569	-69.47934	-69.46122	-69.45132	-71.40889	-71.44434	-71.48414	-71.51411	-71.57909	-71.64997	-71.69591	7671	53.57328	LC08CPF_20200101_20200331_01.04	2020-03-14	WGS84	LGN	CUBIC_CONVOLUTION	NORTH_UP	83	0.0024167	83	0.0005183	0.011436	9	UPPER	2.75	3.333	T1	30	0.99	7.016	1	9	LC82330832020074LGN00	233	15341	15461	5.422	7731	TIRS	4.453	1291195701	LC08_233083_20200314	-0.1	-0.1	WGS84	-0.1	-0.1	LC08RLUT_20150303_20431231_01_12.h5	-0.1	-0.1	-0.1	-0.1	LT8BPF20200310060739_20200324104153.01	4	L1TP	19	-1	LC08_L1TP_233083_20200314_20200325_01_T1	-0.1	15	-50.5235	2e-05	-1	-30.91786	2e-05	-7.68899	2e-05	LPGS_13.1.0	-2.5916	2e-05	-63.49469	-65.01934	-59.91476	2e-05	-57.1787	2e-05	-12.0834	2e-05	2e-05	2e-05	7731	ESTIMATED	30	NADIR	0.1	0702003256647_00027	0.9943464	FINAL	-1	14:33:32.3286510Z	45.21481	LO8BPF20200314140125_20200314144715.02	0.1	-0.001	774.8853	N	N	N	N	N	UTM	Y	OLI_TIRS	Y	480.8883	N	N	233	0.0003342	0.0003342	545	LANDSAT_8	GLS2000	7671	144
Image	B1	PixelType	int	65535	7731	EPSG:32619	0	B2	PixelType	int	65535	7731	EPSG:32619	0	B3	PixelType	int	65535	7731	EPSG:32619	0	B4	PixelType	int	65535	7731	EPSG:32619	0	B5	PixelType	int	65535	7731	EPSG:32619	0	B6	PixelType	int	65535	7731	EPSG:32619	0	B7	PixelType	int	65535	7731	EPSG:32619	0	B8	PixelType	int	65535	15461	EPSG:32619	15	227992.5	-15	-3556493	B9	PixelType	int	65535	7731	EPSG:32619	0	B10	PixelType	int	65535	7731	EPSG:32619	0	B11	PixelType	int	65535	7731	EPSG:32619	0	BQA	PixelType	int	65535	7731	EPSG:32619	0	LANDSAT/LC08/C01/T1/LC08_233083_20200314	-1	0.0061836	0.0015378	0.011983	0.010105	0.012699	0.013004	1201.144	1321.079	LinearRing	-33.00594	-33.13124	-33.33592	-33.60413	-33.75605	-33.82978	-33.83618	-33.91484	-34.08342	-34.23062	-34.20200	-33.41234	-32.62219	-32.55000	-32.51020	-32.12931	-32.12214	-32.26218	-32.36799	-32.59697	-32.84541	-33.00594	7671	53.57328	LC08CPF_20200101_20200331_01.04	2020-03-14	WGS84	LGN	CUBIC_CONVOLUTION	NORTH_UP	83	0.0024167	83	0.0005183	0.011436	9	UPPER	2.75	3.333	T1	30	0.99	7.016	1	9	LC82330832020074LGN00	233	15341	15461	5.422	7731	TIRS	4.453	1291195701	LC08_233083_20200314	-0.1	-0.1	WGS84	-0.1	-0.1	LC08RLUT_20150303_20431231_01_12.h5	-0.1	-0.1	-0.1	-0.1	LT8BPF20200310060739_20200324104153.01	4	L1TP	19	-1	LC08_L1TP_233083_20200314_20200325_01_T1	-0.1	15	-50.5235	2e-05	-1	-30.91786	2e-05	-7.68899	2e-05	LPGS_13.1.0	-2.5916	2e-05	-63.49469	-65.01934	-59.91476	2e-05	-57.1787	2e-05	-12.0834	2e-05	2e-05	2e-05	7731	ESTIMATED	30	NADIR	0.1	0702003256647_00027	0.9943464	FINAL	-1	14:33:32.3286510Z	45.21481	LO8BPF20200314140125_20200314144715.02	0.1	-0.001	774.8853	N	N	N	N	N	UTM	Y	OLI_TIRS	Y	480.8883	N	N	233	0.0003342	0.0003342	545	LANDSAT_8	GLS2000	7671	144
Image	B1	PixelType	int	65535	7671	EPSG:32619	227985	B2	PixelType	int	65535	7671	EPSG:32619	227985	B3	PixelType	int	65535	7671	EPSG:32619	227985	B4	PixelType	int	65535	7671	EPSG:32619	227985	B5	PixelType	int	65535	7671	EPSG:32619	227985	B6	PixelType	int	65535	7671	EPSG:32619	227985	B7	PixelType	int	65535	7671	EPSG:32619	227985	B8	PixelType	int	65535	15341	EPSG:32619	15	227992.5	-15	-3556493	B9	PixelType	int	65535	7671	EPSG:32619	227985	B10	PixelType	int	65535	7671	EPSG:32619	227985	B11	PixelType	int	65535	7671	EPSG:32619	227985	BQA	PixelType	int	65535	7671	EPSG:32619	227985	LANDSAT/LC08/C01/T1/LC08_233083_20200314	-1	0.0061836	0.0015378	0.011983	0.010105	0.012699	0.013004	1201.144	1321.079	LinearRing	-71.69591	-71.73189	-71.79081	-71.86836	-71.91248	-71.93393	-71.93570	-71.54551	-70.68547	-69.90525	-69.89594	-69.68569	-69.47934	-69.46122	-69.45132	-71.40889	-71.44434	-71.48414	-71.51411	-71.57909	-71.64997	-71.69591	7671	53.57328	LC08CPF_20200101_20200331_01.04	2020-03-14	WGS84	LGN	CUBIC_CONVOLUTION	NORTH_UP	83	0.0024167	83	0.0005183	0.011436	9	UPPER	2.75	3.333	T1	30	0.99	7.016	1	9	LC82330832020074LGN00	233	15341	15461	5.422	7731	TIRS	4.453	1291195701	LC08_233083_20200314	-0.1	-0.1	WGS84	-0.1	-0.1	LC08RLUT_20150303_20431231_01_12.h5	-0.1	-0.1	-0.1	-0.1	LT8BPF20200310060739_20200324104153.01	4	L1TP	19	-1	LC08_L1TP_233083_20200314_20200325_01_T1	-0.1	15	-50.5235	2e-05	-1	-30.91786	2e-05	-7.68899	2e-05	LPGS_13.1.0	-2.5916	2e-05	-63.49469	-65.01934	-59.91476	2e-05	-57.1787	2e-05	-12.0834	2e-05	2e-05	2e-05	7731	ESTIMATED	30	NADIR	0.1	0702003256647_00027	0.9943464	FINAL	-1	14:33:32.3286510Z	45.21481	LO8BPF20200314140125_20200314144715.02	0.1	-0.001	774.8853	N	N	N	N	N	UTM	Y	OLI_TIRS	Y	480.8883	N	N	233	0.0003342	0.0003342	545	LANDSAT_8	GLS2000	7671	144
Image	B1	PixelType	int	65535	7731	EPSG:32619	0	B2	PixelType	int	65535	7731	EPSG:32619	0	B3	PixelType	int	65535	7731	EPSG:32619	0	B4	PixelType	int	65535	7731	EPSG:32619	0	B5	PixelType	int	65535	7731	EPSG:32619	0	B6	PixelType	int	65535	7731	EPSG:32619	0	B7	PixelType	int	65535	7731	EPSG:32619	0	B8	PixelType	int	65535	15461	EPSG:32619	15	227992.5	-15	-3556493	B9	PixelType	int	65535	7731	EPSG:32619	0	B10	PixelType	int	65535	7731	EPSG:32619	0	B11	PixelType	int	65535	7731	EPSG:32619	0	BQA	PixelType	int	65535	7731	EPSG:32619	0	LANDSAT/LC08/C01/T1/LC08_233083_20200314	-1	0.0061836	0.0015378	0.011983	0.010105	0.012699	0.013004	1201.144	1321.079	LinearRing	-33.00594	-33.13124	-33.33592	-33.60413	-33.75605	-33.82978	-33.83618	-33.91484	-34.08342	-34.23062	-34.20200	-33.41234	-32.62219	-32.55000	-32.51020	-32.12931	-32.12214	-32.26218	-32.36799	-32.59697	-32.84541	-33.00594	7671	53.57328	LC08CPF_20200101_20200331_01.04	2020-03-14	WGS84	LGN	CUBIC_CONVOLUTION	NORTH_UP	83	0.0024167	83	0.0005183	0.011436	9	UPPER	2.75	3.333	T1	30	0.99	7.016	1	9	LC82330832020074LGN00	233	15341	15461	5.422	7731	TIRS	4.453	1291195701	LC08_233083_20200314	-0.1	-0.1	WGS84	-0.1	-0.1	LC08RLUT_20150303_20431231_01_12.h5	-0.1	-0.1	-0.1	-0.1	LT8BPF20200310060739_20200324104153.01	4	L1TP	19	-1	LC08_L1TP_233083_20200314_20200325_01_T1	-0.1	15	-50.5235	2e-05	-1	-30.91786	2e-05	-7.68899	2e-05	LPGS_13.1.0	-2.5916	2e-05	-63.49469	-65.01934	-59.91476	2e-05	-57.1787	2e-05	-12.0834	2e-05	2e-05	2e-05	7731	ESTIMATED	30	NADIR	0.1	0702003256647_00027	0.9943464	FINAL	-1	14:33:32.3286510Z	45.21481	LO8BPF20200314140125_20200314144715.02	0.1	-0.001	774.8853	N	N	N	N	N	UTM	Y	OLI_TIRS	Y	480.8883	N	N	233	0.0003342	0.0003342	545	LANDSAT_8	GLS2000	7671	144
Image	B1	PixelType	int	65535	7671	EPSG:32619	-30	B2	PixelType	int	65535	7671	EPSG:32619	-30	B3	PixelType	int	65535	7671	EPSG:32619	-30	B4	PixelType	int	65535	7671	EPSG:32619	-30	B5	PixelType	int	65535	7671	EPSG:32619	-30	B6	PixelType	int	65535	7671	EPSG:32619	-30	B7	PixelType	int	65535	7671	EPSG:32619	-30	B8	PixelType	int	65535	15341	EPSG:32619	15	227992.5	-15	-3556493	B9	PixelType	int	65535	7671	EPSG:32619	-30	B10	PixelType	int	65535	7671	EPSG:32619	-30	B11	PixelType	int	65535	7671	EPSG:32619	-30	BQA	PixelType	int	65535	7671	EPSG:32619	-30	LANDSAT/LC08/C01/T1/LC08_233083_20200314	-1	0.0061836	0.0015378	0.011983	0.010105	0.012699	0.013004	1201.144	1321.079	LinearRing	-71.69591	-71.73189	-71.79081	-71.86836	-71.91248	-71.93393	-71.93570	-71.54551	-70.68547	-69.90525	-69.89594	-69.68569	-69.47934	-69.46122	-69.45132	-71.40889	-71.44434	-71.48414	-71.51411	-71.57909	-71.64997	-71.69591	7671	53.57328	LC08CPF_20200101_20200331_01.04	2020-03-14	WGS84	LGN	CUBIC_CONVOLUTION	NORTH_UP	83	0.0024167	83	0.0005183	0.011436	9	UPPER	2.75	3.333	T1	30	0.99	7.016	1	9	LC82330832020074LGN00	233	15341	15461	5.422	7731	TIRS	4.453	1291195701	LC08_233083_20200314	-0.1	-0.1	WGS84	-0.1	-0.1	LC08RLUT_20150303_20431231_01_12.h5	-0.1	-0.1	-0.1	-0.1	LT8BPF20200310060739_20200324104153.01	4	L1TP	19	-1	LC08_L1TP_233083_20200314_20200325_01_T1	-0.1	15	-50.5235	2e-05	-1	-30.91786	2e-05	-7.68899	2e-05	LPGS_13.1.0	-2.5916	2e-05	-63.49469	-65.01934	-59.91476	2e-05	-57.1787	2e-05	-12.0834	2e-05	2e-05	2e-05	7731	ESTIMATED	30	NADIR	0.1	0702003256647_00027	0.9943464	FINAL	-1	14:33:32.3286510Z	45.21481	LO8BPF20200314140125_20200314144715.02	0.1	-0.001	774.8853	N	N	N	N	N	UTM	Y	OLI_TIRS	Y	480.8883	N	N	233	0.0003342	0.0003342	545	LANDSAT_8	GLS2000	7671	144
Image	B1	PixelType	int	65535	7731	EPSG:32619	-3556485	B2	PixelType	int	65535	7731	EPSG:32619	-3556485	B3	PixelType	int	65535	7731	EPSG:32619	-3556485	B4	PixelType	int	65535	7731	EPSG:32619	-3556485	B5	PixelType	int	65535	7731	EPSG:32619	-3556485	B6	PixelType	int	65535	7731	EPSG:32619	-3556485	B7	PixelType	int	65535	7731	EPSG:32619	-3556485	B8	PixelType	int	65535	15461	EPSG:32619	15	227992.5	-15	-3556493	B9	PixelType	int	65535	7731	EPSG:32619	-3556485	B10	PixelType	int	65535	7731	EPSG:32619	-3556485	B11	PixelType	int	65535	7731	EPSG:32619	-3556485	BQA	PixelType	int	65535	7731	EPSG:32619	-3556485	LANDSAT/LC08/C01/T1/LC08_233083_20200314	-1	0.0061836	0.0015378	0.011983	0.010105	0.012699	0.013004	1201.144	1321.079	LinearRing	-33.00594	-33.13124	-33.33592	-33.60413	-33.75605	-33.82978	-33.83618	-33.91484	-34.08342	-34.23062	-34.20200	-33.41234	-32.62219	-32.55000	-32.51020	-32.12931	-32.12214	-32.26218	-32.36799	-32.59697	-32.84541	-33.00594	7671	53.57328	LC08CPF_20200101_20200331_01.04	2020-03-14	WGS84	LGN	CUBIC_CONVOLUTION	NORTH_UP	83	0.0024167	83	0.0005183	0.011436	9	UPPER	2.75	3.333	T1	30	0.99	7.016	1	9	LC82330832020074LGN00	233	15341	15461	5.422	7731	TIRS	4.453	1291195701	LC08_233083_20200314	-0.1	-0.1	WGS84	-0.1	-0.1	LC08RLUT_20150303_20431231_01_12.h5	-0.1	-0.1	-0.1	-0.1	LT8BPF20200310060739_20200324104153.01	4	L1TP	19	-1	LC08_L1TP_233083_20200314_20200325_01_T1	-0.1	15	-50.5235	2e-05	-1	-30.91786	2e-05	-7.68899	2e-05	LPGS_13.1.0	-2.5916	2e-05	-63.49469	-65.01934	-59.91476	2e-05	-57.1787	2e-05	-12.0834	2e-05	2e-05	2e-05	7731	ESTIMATED	30	NADIR	0.1	0702003256647_00027	0.9943464	FINAL	-1	14:33:32.3286510Z	45.21481	LO8BPF20200314140125_20200314144715.02	0.1	-0.001	774.8853	N	N	N	N	N	UTM	Y	OLI_TIRS	Y	480.8883	N	N	233	0.0003342	0.0003342	545	LANDSAT_8	GLS2000	7671	144

library(rgee)
library(mapedit) # (OPTIONAL) Interactive editing of vector data
library(raster)  # Manipulate raster data
library(scales)  # Scale functions for visualization
library(cptcity)  # cptcity color gradients!
#library(tmap)    # Thematic Map Visualization <3

ee_Initialize()

## -- rgee 1.0.9 --------------------------------------- earthengine-api 0.1.259 -- 
##  v email: not_defined
##  v Initializing Google Earth Engine:
 v Initializing Google Earth Engine:  DONE!
## 
 v Earth Engine user: users/tarredwall 
## --------------------------------------------------------------------------------

object <- ee$ImageCollection("LANDSAT/LC08/C01/T1")$
  filter(ee$Filter()$eq("WRS_PATH", 233))$
  filter(ee$Filter()$eq("WRS_ROW", 83))$
  filterMetadata ('CLOUD_COVER', 'Less_Than', 20)$
  filterDate("2020-03-01", "2020-04-01")$aside(ee_print)

## ------------------------------------------------ Earth Engine ImageCollection --
## ImageCollection Metadata:
##  - Class                      : ee$ImageCollection
##  - Number of Images           : 2
##  - Number of Properties       : 38
##  - Number of Pixels*          : 1415544744
##  - Approximate size*          : 28.86 GB
## Image Metadata (img_index = 0):
##  - ID                         : LANDSAT/LC08/C01/T1/LC08_233083_20200314
##  - system:time_start          : 2020-03-14 14:33:32
##  - system:time_end            : 2020-03-14 14:33:32
##  - Number of Bands            : 12
##  - Bands names                : B1 B2 B3 B4 B5 B6 B7 B8 B9 B10 B11 BQA
##  - Number of Properties       : 119
##  - Number of Pixels*          : 707772372
##  - Approximate size*          : 14.43 GB
## Band Metadata (img_band = 'B1'):
##  - EPSG (SRID)                : WGS 84 / UTM zone 19N (EPSG:32619)
##  - proj4string                : +proj=utm +zone=19 +datum=WGS84 +units=m +no_defs
##  - Geotransform               : 30 0 227985 0 -30 -3556485
##  - Nominal scale (meters)     : 30
##  - Dimensions                 : 7643 7717
##  - Number of Pixels           : 58981031
##  - Data type                  : INT
##  - Approximate size           : 1.20 GB
##  --------------------------------------------------------------------------------
##  NOTE: (*) Properties calculated considering a constant  geotransform and data type.

library(kableExtra)
addNDVI <- function(image) {
  return(image$addBands(image$normalizedDifference(c("B5", "B4"))))
}
ndviCollection <- object$map(addNDVI)
first <- ndviCollection$first()
a_001 <- first$getInfo()
a_001 <- as.data.frame(a_001)

kbl(a_001[,c(2:9)]) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")

bands.id	bands.data_type.type	bands.data_type.precision	bands.data_type.max	bands.dimensions	bands.crs	bands.crs_transform
B1	PixelType	int	65535	7671	EPSG:32619	30
B1	PixelType	int	65535	7731	EPSG:32619	0
B1	PixelType	int	65535	7671	EPSG:32619	227985
B1	PixelType	int	65535	7731	EPSG:32619	0
B1	PixelType	int	65535	7671	EPSG:32619	-30
B1	PixelType	int	65535	7731	EPSG:32619	-3556485

kbl(a_001[,c(10:17)])%>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")

bands.id.1	bands.data_type.type.1	bands.data_type.precision.1	bands.data_type.max.1	bands.dimensions.1	bands.crs.1	bands.crs_transform.1
B2	PixelType	int	65535	7671	EPSG:32619	30
B2	PixelType	int	65535	7731	EPSG:32619	0
B2	PixelType	int	65535	7671	EPSG:32619	227985
B2	PixelType	int	65535	7731	EPSG:32619	0
B2	PixelType	int	65535	7671	EPSG:32619	-30
B2	PixelType	int	65535	7731	EPSG:32619	-3556485

kbl(a_001[,c(18:25)])%>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")

bands.id.2	bands.data_type.type.2	bands.data_type.precision.2	bands.data_type.max.2	bands.dimensions.2	bands.crs.2	bands.crs_transform.2
B3	PixelType	int	65535	7671	EPSG:32619	30
B3	PixelType	int	65535	7731	EPSG:32619	0
B3	PixelType	int	65535	7671	EPSG:32619	227985
B3	PixelType	int	65535	7731	EPSG:32619	0
B3	PixelType	int	65535	7671	EPSG:32619	-30
B3	PixelType	int	65535	7731	EPSG:32619	-3556485

kbl(a_001[,c(26:33)])%>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")

bands.id.3	bands.data_type.type.3	bands.data_type.precision.3	bands.data_type.max.3	bands.dimensions.3	bands.crs.3	bands.crs_transform.3
B4	PixelType	int	65535	7671	EPSG:32619	30
B4	PixelType	int	65535	7731	EPSG:32619	0
B4	PixelType	int	65535	7671	EPSG:32619	227985
B4	PixelType	int	65535	7731	EPSG:32619	0
B4	PixelType	int	65535	7671	EPSG:32619	-30
B4	PixelType	int	65535	7731	EPSG:32619	-3556485

kbl(a_001[,c(34:41)])%>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")

bands.id.4	bands.data_type.type.4	bands.data_type.precision.4	bands.data_type.max.4	bands.dimensions.4	bands.crs.4	bands.crs_transform.4
B5	PixelType	int	65535	7671	EPSG:32619	30
B5	PixelType	int	65535	7731	EPSG:32619	0
B5	PixelType	int	65535	7671	EPSG:32619	227985
B5	PixelType	int	65535	7731	EPSG:32619	0
B5	PixelType	int	65535	7671	EPSG:32619	-30
B5	PixelType	int	65535	7731	EPSG:32619	-3556485

kbl(a_001[,c(42:49)])%>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")

bands.id.5	bands.data_type.type.5	bands.data_type.precision.5	bands.data_type.max.5	bands.dimensions.5	bands.crs.5	bands.crs_transform.5
B6	PixelType	int	65535	7671	EPSG:32619	30
B6	PixelType	int	65535	7731	EPSG:32619	0
B6	PixelType	int	65535	7671	EPSG:32619	227985
B6	PixelType	int	65535	7731	EPSG:32619	0
B6	PixelType	int	65535	7671	EPSG:32619	-30
B6	PixelType	int	65535	7731	EPSG:32619	-3556485

kbl(a_001[,c(50:57)])%>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")

bands.id.6	bands.data_type.type.6	bands.data_type.precision.6	bands.data_type.max.6	bands.dimensions.6	bands.crs.6	bands.crs_transform.6
B7	PixelType	int	65535	7671	EPSG:32619	30
B7	PixelType	int	65535	7731	EPSG:32619	0
B7	PixelType	int	65535	7671	EPSG:32619	227985
B7	PixelType	int	65535	7731	EPSG:32619	0
B7	PixelType	int	65535	7671	EPSG:32619	-30
B7	PixelType	int	65535	7731	EPSG:32619	-3556485

kbl(a_001[,c(58:70)]) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")

bands.id.7	bands.data_type.type.7	bands.data_type.precision.7	bands.data_type.max.7	bands.dimensions.7	bands.crs.7	bands.crs_transform.15L	bands.crs_transform.227992.5	bands.crs_transform..15L	bands.crs_transform..3556492.5
B8	PixelType	int	65535	15341	EPSG:32619	15	227992.5	-15	-3556493
B8	PixelType	int	65535	15461	EPSG:32619	15	227992.5	-15	-3556493
B8	PixelType	int	65535	15341	EPSG:32619	15	227992.5	-15	-3556493
B8	PixelType	int	65535	15461	EPSG:32619	15	227992.5	-15	-3556493
B8	PixelType	int	65535	15341	EPSG:32619	15	227992.5	-15	-3556493
B8	PixelType	int	65535	15461	EPSG:32619	15	227992.5	-15	-3556493

kbl(a_001[,c(71:78)]) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")

bands.id.8	bands.data_type.type.8	bands.data_type.precision.8	bands.data_type.max.8	bands.dimensions.8	bands.crs.8	bands.crs_transform.7
B9	PixelType	int	65535	7671	EPSG:32619	30
B9	PixelType	int	65535	7731	EPSG:32619	0
B9	PixelType	int	65535	7671	EPSG:32619	227985
B9	PixelType	int	65535	7731	EPSG:32619	0
B9	PixelType	int	65535	7671	EPSG:32619	-30
B9	PixelType	int	65535	7731	EPSG:32619	-3556485

kbl(a_001[,c(79:86)]) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")

bands.id.9	bands.data_type.type.9	bands.data_type.precision.9	bands.data_type.max.9	bands.dimensions.9	bands.crs.9	bands.crs_transform.8
B10	PixelType	int	65535	7671	EPSG:32619	30
B10	PixelType	int	65535	7731	EPSG:32619	0
B10	PixelType	int	65535	7671	EPSG:32619	227985
B10	PixelType	int	65535	7731	EPSG:32619	0
B10	PixelType	int	65535	7671	EPSG:32619	-30
B10	PixelType	int	65535	7731	EPSG:32619	-3556485

kbl(a_001[,c(87:110)]) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")

bands.id.10	bands.data_type.type.10	bands.data_type.precision.10	bands.data_type.max.10	bands.dimensions.10	bands.crs.10	bands.crs_transform.9	bands.id.11	bands.data_type.type.11	bands.data_type.precision.11	bands.data_type.max.11	bands.dimensions.11	bands.crs.11	bands.crs_transform.10	bands.id.12	bands.data_type.type.12	bands.data_type.precision.12	bands.data_type.min.12	bands.data_type.max.12	bands.dimensions.12	bands.crs.12	bands.crs_transform.11
B11	PixelType	int	65535	7671	EPSG:32619	30	BQA	PixelType	int	65535	7671	EPSG:32619	30	nd	PixelType	float	-1	1	7671	EPSG:32619	30
B11	PixelType	int	65535	7731	EPSG:32619	0	BQA	PixelType	int	65535	7731	EPSG:32619	0	nd	PixelType	float	-1	1	7731	EPSG:32619	0
B11	PixelType	int	65535	7671	EPSG:32619	227985	BQA	PixelType	int	65535	7671	EPSG:32619	227985	nd	PixelType	float	-1	1	7671	EPSG:32619	227985
B11	PixelType	int	65535	7731	EPSG:32619	0	BQA	PixelType	int	65535	7731	EPSG:32619	0	nd	PixelType	float	-1	1	7731	EPSG:32619	0
B11	PixelType	int	65535	7671	EPSG:32619	-30	BQA	PixelType	int	65535	7671	EPSG:32619	-30	nd	PixelType	float	-1	1	7671	EPSG:32619	-30
B11	PixelType	int	65535	7731	EPSG:32619	-3556485	BQA	PixelType	int	65535	7731	EPSG:32619	-3556485	nd	PixelType	float	-1	1	7731	EPSG:32619	-3556485

La columna 111 del dataset nos da la identificación de la scene:

kbl(a_001[,c(111)]) %>%
kable_styling(bootstrap_options = c("striped", "hover")) %>%
kable_paper() %>%
scroll_box(width = "100%", height = "300px")

x
LANDSAT/LC08/C01/T1/LC08_233083_20200314
LANDSAT/LC08/C01/T1/LC08_233083_20200314
LANDSAT/LC08/C01/T1/LC08_233083_20200314
LANDSAT/LC08/C01/T1/LC08_233083_20200314
LANDSAT/LC08/C01/T1/LC08_233083_20200314
LANDSAT/LC08/C01/T1/LC08_233083_20200314

1.4 Vamos a analizar la información extraíble de una imagen LANDSAT/LC08/C01/T1.

Propiedades de las imágenes

1.1.1 Descripción general de los datos de una imagen.

Despleguemos las primeras 30 columnas de la base de datos:

a_001 <- ee$Image("LANDSAT/LC08/C01/T1/LC08_044034_20140403") %>% ee$Image$getInfo()
write.csv2(a_001, "a_001.csv")
# saveRDS(g,"g.rds")
a_001 <- as.data.frame(a_001)
head(colnames(a_001),30)

##  [1] "type"                        "bands.id"                   
##  [3] "bands.data_type.type"        "bands.data_type.precision"  
##  [5] "bands.data_type.min"         "bands.data_type.max"        
##  [7] "bands.dimensions"            "bands.crs"                  
##  [9] "bands.crs_transform"         "bands.id.1"                 
## [11] "bands.data_type.type.1"      "bands.data_type.precision.1"
## [13] "bands.data_type.min.1"       "bands.data_type.max.1"      
## [15] "bands.dimensions.1"          "bands.crs.1"                
## [17] "bands.crs_transform.1"       "bands.id.2"                 
## [19] "bands.data_type.type.2"      "bands.data_type.precision.2"
## [21] "bands.data_type.min.2"       "bands.data_type.max.2"      
## [23] "bands.dimensions.2"          "bands.crs.2"                
## [25] "bands.crs_transform.2"       "bands.id.3"                 
## [27] "bands.data_type.type.3"      "bands.data_type.precision.3"
## [29] "bands.data_type.min.3"       "bands.data_type.max.3"

Que se corresponden con las propiedades de las imagenes contenidas en: USGS Landsat 8 Collection 1 Tier 1 Raw Scenes https://developers.google.com/earth-engine/datasets/catalog/LANDSAT_LC08_C01_T1#image-properties

1.5 Obtención de una escena

object <- ee$ImageCollection("LANDSAT/LC08/C01/T1")$
  filter(ee$Filter()$eq("WRS_PATH", 233))$
  filter(ee$Filter()$eq("WRS_ROW", 83))$
  filterMetadata ('CLOUD_COVER', 'Less_Than', 20)$
  filterDate("2020-03-01", "2020-04-01")
# Representa la imagen a traves de una composicion RGB
first <- object$first()

vizParams <- list(
  bands = c("B5", "B4", "B3"),
  min = 5000,
  max = 15000,
  gamma = 1.3
)
Map$setCenter(lon = -70.64724, lat = -33.47269, zoom = 8)
Map$addLayer(first, vizParams, "Landsat 8")

library(rgee)
ee_Initialize()

## -- rgee 1.0.9 --------------------------------------- earthengine-api 0.1.259 -- 
##  v email: not_defined
##  v Initializing Google Earth Engine:
 v Initializing Google Earth Engine:  DONE!
## 
 v Earth Engine user: users/tarredwall 
## --------------------------------------------------------------------------------

l8 <- ee$ImageCollection("LANDSAT/LC08/C01/T1")
image_2015 <- ee$Algorithms$Landsat$simpleComposite(
  collection = l8$filterDate("2015-01-01", "2015-12-31"),
  asFloat = TRUE
)
bands <- c( "B5", "B6")
getNDBI <- function(image) {
  image$normalizedDifference(c("B5", "B4"))
}
ndbi1 <- getNDBI(image_2015)
ndbiParams <- list(palette = c(
  "#000000", "#ff0000", "#ffffff","#0000ff"
))
landMask <- ee$Image("CGIAR/SRTM90_V4")$mask()
maskedDifference <- ndbi1$updateMask(landMask)
Map$setCenter(lon = -70.64724, lat = -33.47269, zoom = 10)
Map$addLayer(maskedDifference, ndbiParams, "NDBI")

Descarga de precipitación desde GEE https://rpubs.com/Bonilla-Escoto/720385

library(rgee)
library(sf)

## Warning: package 'sf' was built under R version 4.0.5

## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1

require(tidyverse)

## Loading required package: tidyverse

## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --

## v ggplot2 3.3.2     v purrr   0.3.4
## v tibble  3.0.1     v dplyr   1.0.2
## v tidyr   1.1.0     v stringr 1.4.0
## v readr   1.3.1     v forcats 0.5.0

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x readr::col_factor() masks scales::col_factor()
## x purrr::discard()    masks scales::discard()
## x tidyr::extract()    masks raster::extract()
## x dplyr::filter()     masks stats::filter()
## x dplyr::group_rows() masks kableExtra::group_rows()
## x dplyr::lag()        masks stats::lag()
## x dplyr::select()     masks raster::select()

punto <- st_point(c(617175,1428076), dim="XY")
puntoUTM <- st_sfc(punto,crs=32616)
puntoGeo <- st_transform(puntoUTM,4326)
punto_ee <- sf_as_ee(puntoGeo)

gee1 <- ee$ImageCollection("ECMWF/ERA5/MONTHLY")$select("total_precipitation")$
  filterDate("2005-01-01", "2020-01-01")

datos_gee1 <- ee_extract(x=gee1,y=punto_ee,scale = 250,sf=T)
st_geometry(datos_gee1) <- NULL

datos_t1 <- as.data.frame(t(datos_gee1[,-1]))
head(datos_t1)

##                                     1
## X200502_total_precipitation 0.0029145
## X200503_total_precipitation 0.0158447
## X200504_total_precipitation 0.0500956
## X200505_total_precipitation 0.1628986
## X200506_total_precipitation 0.2862051
## X200507_total_precipitation 0.1587650

sal <- rownames_to_column(datos_t1, var="fecha_apilada")

colnames(sal) <- c("Fecha_apilada","Valores")
Fecha <- seq.Date(from = as.Date("01-01-2005", format="%d-%m-%Y", sep="-"),
                  to=as.Date("01-12-2019", format="%d-%m-%Y",sep="-"),
                  by="month")

#sal$Fecha <- Fecha

sal$Valores_mm <- sal$Valores*1000

mensual <- as.data.frame(matrix(sal$Valores_mm, ncol = 12, byrow = T))

## Warning in matrix(sal$Valores_mm, ncol = 12, byrow = T): la longitud de los
## datos [179] no es un submúltiplo o múltiplo del número de filas [15] en la
## matriz

colnames(mensual) <- c("Ene","Feb","Mar","Abr","May","Jun","Jul","Ago","Sep","Oct","Nov","Dic")
mensual$Anio <- seq.int(from=2005, to=2019,by=1) 
head(mensual)

##       Ene     Feb      Mar      Abr      May      Jun      Jul      Ago
## 1  2.9145 15.8447  50.0956 162.8986 286.2051 158.7650 156.5607 130.7581
## 2 14.5964  7.5814  17.5079  92.0442 180.7454  95.7995  77.7299  99.9864
## 3  5.1926 17.1376  55.7835 193.5994 126.5605  97.5545 176.3386 203.9109
## 4 13.0933  9.3654  16.8935 255.3760 156.1875 162.6278 152.4938 219.4597
## 5  9.9886  7.9723  12.8181 257.8217 223.0067  73.8023  82.5352 109.9282
## 6 24.8376 17.5173 131.0137 431.3926 170.1693 179.8437 423.4943 269.9235
##        Sep      Oct     Nov     Dic Anio
## 1 195.7654 101.2689 23.9303 22.9142 2005
## 2 155.8129  43.0326 46.5991 16.9111 2006
## 3 321.9736  49.9877 18.0156 18.9402 2007
## 4 287.7165  33.1901 15.6414 10.3437 2008
## 5 108.7499  42.6289 16.7708  9.9737 2009
## 6 117.7062  61.2735  5.1110 14.4748 2010

p<-ggplot(data=mensual, aes(x=Anio, y=Dic)) +
  geom_bar(stat="identity", fill="steelblue")+
  geom_text(aes(label=Dic), vjust=1.6, color="black", size=3.5)+
  theme_minimal()
p

Construimos un histograma del Sentinel

Quiero muestrear píxeles aleatoriamente dentro de una imagen que tiene tres bandas, en un conjunto de regiones diferentes. Con esos datos muestreados, quiero crear un gráfico que contenga histogramas para cada región de muestra, para una banda determinada. Entonces tendré tres parcelas (una para cada banda).

En mi caso, una región representa un tipo específico de hábitat dentro de la imagen, por lo que tengo un polígono para cada tipo de hábitat.

library(ggplot2)
###############################################################
# Create a geometry for each habitat type that I want to sample
###############################################################

unimproved_grazing = ee$Geometry$Polygon(
  list(
    c(-1.345864517292552, 59.97338078318311),
    c(-1.345864517292552, 59.97237144577356),
    c(-1.3438904114575911, 59.97237144577356),
    c(-1.3438904114575911, 59.97338078318311)
  )
)

improved_cut = ee$Geometry$Polygon(
  list(
    c(-1.329604032461742, 59.92455155105425),
    c(-1.329604032461742, 59.923922477208166),
    c(-1.3280268935609851, 59.923922477208166),
    c(-1.3280268935609851, 59.92455155105425)
  )
)

heath = ee$Geometry$Polygon(
  list(
    c(-1.351189448199126, 59.94125887709115),
    c(-1.351189448199126, 59.940151826756974),
    c(-1.3490865963314502, 59.940151826756974),
    c(-1.3490865963314502, 59.94125887709115)
  )
)

arable = ee$Geometry$Polygon(
  list(
    c(-1.3263107736431534, 59.969782903189916),
    c(-1.3263107736431534, 59.96911710113049),
    c(-1.325023313316005, 59.96911710113049),
    c(-1.325023313316005, 59.969782903189916)
  )
)

##################################################
# Create feature collection of all habitat regions
##################################################

habitats = ee$FeatureCollection(
  list(
    ee$Feature(unimproved_grazing, list(name = "unimproved")),
    ee$Feature(improved_cut,       list(name = "improved")),
    ee$Feature(heath,              list(name = "heath")),
    ee$Feature(arable,             list(name = "arable"))  
  )
)


sentinel1 = ee$ImageCollection('COPERNICUS/S1_GRD')$
  # Filter date range - one year of images
  filterDate("2019-03-01", "2020-03-01")

# Filter by metadata properties
vvvh = sentinel1$
  filter(ee$Filter$listContains('transmitterReceiverPolarisation', 'VV'))$
  filter(ee$Filter$listContains('transmitterReceiverPolarisation', 'VH'))$
  filter(ee$Filter$eq('instrumentMode', 'IW'))$
  filter(ee$Filter$eq('resolution_meters', 10))

# Filter to get images from different look angles.
vhAscending = vvvh$filter(ee$Filter$eq('orbitProperties_pass', 'ASCENDING'))
vhDescending = vvvh$filter(ee$Filter$eq('orbitProperties_pass', 'DESCENDING'))

# Create a composite from means at different polarizations and look angles.
# Calculate mean of VH ascending polarisation
VH_Ascending_mean <- vhAscending$select("VH")$mean()
# Calculate mean of VH descending polarisation
VH_Descending_mean <- vhDescending$select("VH")$mean()
# Merge VV and VH and calculate mean
VV_Ascending_Descending_mean <- vhAscending$select("VV") %>%
  ee$ImageCollection$merge(vhDescending$select("VV")) %>%
  ee$ImageCollection$mean()

# Create single image containing bands of interest
collection <- ee$ImageCollection$fromImages(list(
  VH_Ascending_mean,
  VV_Ascending_Descending_mean,
  VH_Descending_mean
))$toBands()

######################################
# Sample pixels in each habitat region
######################################
mySamples = collection$sampleRegions(
                    collection = habitats,
                    scale= 10
                    )


# ee_help(collection$sampleRegions)

for (index in 0:3) {
  habitat_f <- habitats %>% ee_get(index)
  mySamples <- collection$sampleRegions(
    collection = ee$Feature(habitat_f$first()),
    scale= 10,
    geometries=TRUE
  )
  cat(sprintf("Processing geom: %02d \n", index + 1))
  mySamples_sf <- ee_as_sf(mySamples)
  if (index == 0) {
    mySamples_sf_batch <- mySamples_sf   
  } else {
    mySamples_sf_batch <- rbind(mySamples_sf_batch, mySamples_sf)
  }
}

## Processing geom: 01 
## Processing geom: 02 
## Processing geom: 03 
## Processing geom: 04

ggplot(mySamples_sf_batch, aes(X1_VV)) +
  geom_histogram(color="darkblue", fill="lightblue") +
  facet_grid(~name)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Obtener los valores de la base de datos:

La misión Sentinel-1 proporciona datos de un instrumento de radar de apertura sintética (SAR) de banda C de doble polarización. Esta colección incluye las escenas S1 Ground Range Detected (GRD), procesadas con la caja de herramientas Sentinel-1 para generar un producto calibrado y orto-corregido. La colección se actualiza diariamente. Los nuevos activos se ingieren dentro de los dos días posteriores a su disponibilidad.

Esta colección contiene todas las escenas GRD. Cada escena tiene una de 3 resoluciones (10, 25 o 40 metros), 4 combinaciones de bandas (correspondientes a la polarización de la escena) y 3 modos de instrumento. El uso de la colección en un contexto de mosaico probablemente requerirá filtrar hasta un conjunto homogéneo de bandas y parámetros. Consulte este artículo para obtener detalles sobre el uso y preprocesamiento de la colección. Cada escena contiene 1 o 2 de las 4 posibles bandas de polarización, según la configuración de polarización del instrumento. Las combinaciones posibles son VV o HH de banda única y VV + VH y HH + HV de doble banda:

VV: copolarización única, transmisión vertical / recepción vertical
HH: copolarización única, transmisión horizontal / recepción horizontal
VV + VH: polarización cruzada de doble banda, transmisión vertical / recepción horizontal
HH + HV: polarización cruzada de doble banda, transmisión horizontal / recepción vertical

Cada escena también incluye una banda de ‘ángulo’ adicional que contiene el ángulo de incidencia de visión aproximado en grados en cada punto. Esta banda se genera interpolando la propiedad ‘IncideAngle’ del campo cuadriculado ‘geolocationGridPoint’ proporcionado con cada activo.

Cada escena se preprocesó con Sentinel-1 Toolbox siguiendo los siguientes pasos:

Eliminación de ruido térmico
Calibración radiométrica
Corrección de terreno usando SRTM 30 o ASTER DEM para áreas de más de 60 grados de latitud, donde SRTM no está disponible. Los valores finales corregidos por el terreno se convierten a decibelios mediante una escala logarítmica (10 * log10 (x)).

Para obtener más información sobre estos pasos de procesamiento previo, consulte el artículo de procesamiento previo de Sentinel-1. Para obtener más consejos sobre cómo trabajar con imágenes de Sentinel-1, consulte el tutorial de Guido Lemoine sobre conceptos básicos de SAR y el tutorial de Mort Canty sobre detección de cambios de SAR.

Esta colección se calcula sobre la marcha. Si desea utilizar la colección subyacente con valores de potencia sin procesar (que se actualiza más rápido), consulte COPERNICUS / S1_GRD_FLOAT.

  # Construímos una geometría que describe un polígono:
  habitats = ee$FeatureCollection(
    list(
      ee$Feature(ee$Geometry$Polygon(
        list(
          c(-1.345864517292552, 59.97338078318311),
          c(-1.345864517292552, 59.97237144577356),
          c(-1.3438904114575911, 59.97237144577356),
          c(-1.3438904114575911, 59.97338078318311)
        )
      ), list(name = "unimproved"))
    )
  )

  sentinel1 = ee$ImageCollection('COPERNICUS/S1_GRD')$filterDate("2019-03-01", "2020-03-01")

  # Filter by metadata properties
  vvvh = sentinel1$
    filter(ee$Filter$listContains('transmitterReceiverPolarisation', 'VV'))$
    filter(ee$Filter$listContains('transmitterReceiverPolarisation', 'VH'))$
    filter(ee$Filter$eq('instrumentMode', 'IW'))$
    filter(ee$Filter$eq('resolution_meters', 10))

  # Filter to get images from different look angles.
  vhAscending = vvvh$filter(ee$Filter$eq('orbitProperties_pass', 'ASCENDING'))
  vhDescending = vvvh$filter(ee$Filter$eq('orbitProperties_pass', 'DESCENDING'))

  # Calculate mean of VH ascending polarisation
  VH_Ascending_mean <- vhAscending$select("VH")$mean()
  
  # Calculate mean of VH descending polarisation
  VH_Descending_mean <- vhDescending$select("VH")$mean()
  
  # Merge VV and VH and calculate mean
  VV_Ascending_Descending_mean <- vhAscending$select("VV") %>%
    ee$ImageCollection$merge(vhDescending$select("VV")) %>%
    ee$ImageCollection$mean()

  # Create single image containing bands of interest
  collection <- ee$ImageCollection$fromImages(list(VH_Ascending_mean,VV_Ascending_Descending_mean,VH_Descending_mean))$toBands()

  # importante: para obtener ayuda de R:
  ee_help(ee$Geometry$Polygon)

  mySamples <- collection$sampleRegions(
    collection = ee$Feature(habitats$first()),
    scale= 10,
    geometries=TRUE
  )

  mySamples_sf_batch <-  ee_as_sf(mySamples)  

  ggplot(mySamples_sf_batch, aes(X1_VV)) +
    geom_histogram(color="darkblue", fill="lightblue") +
    facet_grid(~name)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

a <- as.data.frame(mySamples_sf)
a  %>%  kbl() %>%
 kable_material(c("striped", "hover"), font_size = 12)%>%
  scroll_box(width = "100%", height = "200px")

X0_VH	X1_VV	X2_VH	name	geometry
-20.02172	-12.89675	-20.32706	arable	POINT (-1.326228 59.96915)
-20.12092	-13.00383	-20.41705	arable	POINT (-1.326138 59.96915)
-20.21582	-13.20390	-20.36710	arable	POINT (-1.326048 59.96915)
-20.45528	-13.35599	-20.69379	arable	POINT (-1.325958 59.96915)
-20.52344	-13.37268	-20.48843	arable	POINT (-1.325868 59.96915)
-20.56013	-13.30978	-20.57392	arable	POINT (-1.325779 59.96915)
-20.35753	-13.20408	-20.52587	arable	POINT (-1.325689 59.96915)
-20.16272	-13.09865	-20.35889	arable	POINT (-1.325599 59.96915)
-20.14003	-12.92395	-20.36939	arable	POINT (-1.325509 59.96915)
-19.91871	-12.75071	-20.24559	arable	POINT (-1.325419 59.96915)
-19.82414	-12.71546	-20.02239	arable	POINT (-1.32533 59.96915)
-19.71656	-12.62275	-19.61642	arable	POINT (-1.32524 59.96915)
-19.60558	-12.64772	-19.75683	arable	POINT (-1.32515 59.96915)
-19.60592	-12.68841	-19.62177	arable	POINT (-1.32506 59.96915)
-20.26387	-12.98642	-20.18087	arable	POINT (-1.326228 59.96924)
-20.23379	-13.04805	-20.24608	arable	POINT (-1.326138 59.96924)
-20.24295	-13.11691	-20.27763	arable	POINT (-1.326048 59.96924)
-20.25515	-13.33413	-20.74608	arable	POINT (-1.325958 59.96924)
-20.38536	-13.32222	-20.54083	arable	POINT (-1.325868 59.96924)
-20.41444	-13.34241	-20.36454	arable	POINT (-1.325779 59.96924)
-20.35862	-13.30075	-20.10031	arable	POINT (-1.325689 59.96924)
-20.32328	-13.18406	-20.10581	arable	POINT (-1.325599 59.96924)
-20.32923	-13.00434	-20.34064	arable	POINT (-1.325509 59.96924)
-20.21219	-12.82865	-20.30770	arable	POINT (-1.325419 59.96924)
-20.20127	-12.76990	-20.18783	arable	POINT (-1.32533 59.96924)
-20.12188	-12.62799	-19.97489	arable	POINT (-1.32524 59.96924)
-19.89861	-12.70538	-19.88010	arable	POINT (-1.32515 59.96924)
-19.76500	-12.63888	-19.66829	arable	POINT (-1.32506 59.96924)
-20.46296	-13.24092	-20.25372	arable	POINT (-1.326228 59.96933)
-20.47604	-13.33872	-20.30657	arable	POINT (-1.326138 59.96933)
-20.52329	-13.28309	-20.23974	arable	POINT (-1.326048 59.96933)
-20.47729	-13.37945	-20.45973	arable	POINT (-1.325958 59.96933)
-20.31831	-13.32646	-20.42256	arable	POINT (-1.325868 59.96933)
-20.30166	-13.40445	-20.19303	arable	POINT (-1.325779 59.96933)
-20.22769	-13.32022	-20.05328	arable	POINT (-1.325689 59.96933)
-20.23566	-13.23012	-20.20655	arable	POINT (-1.325599 59.96933)
-20.13489	-13.12533	-20.54698	arable	POINT (-1.325509 59.96933)
-20.07454	-13.01608	-20.72703	arable	POINT (-1.325419 59.96933)
-20.13259	-12.95435	-20.53161	arable	POINT (-1.32533 59.96933)
-20.15818	-12.80701	-20.29710	arable	POINT (-1.32524 59.96933)
-20.03196	-12.91492	-20.14620	arable	POINT (-1.32515 59.96933)
-19.76623	-12.77865	-19.82271	arable	POINT (-1.32506 59.96933)
-20.91416	-13.36608	-20.43523	arable	POINT (-1.326228 59.96942)
-20.85911	-13.52073	-20.46222	arable	POINT (-1.326138 59.96942)
-20.77483	-13.38483	-20.37609	arable	POINT (-1.326048 59.96942)
-20.79339	-13.34576	-20.23170	arable	POINT (-1.325958 59.96942)
-20.48414	-13.23775	-20.26956	arable	POINT (-1.325868 59.96942)
-20.40943	-13.26491	-20.10286	arable	POINT (-1.325779 59.96942)
-20.29028	-13.32862	-20.21315	arable	POINT (-1.325689 59.96942)
-20.24129	-13.31490	-20.49951	arable	POINT (-1.325599 59.96942)
-20.11229	-13.32264	-20.77931	arable	POINT (-1.325509 59.96942)
-19.93085	-13.29105	-20.89190	arable	POINT (-1.325419 59.96942)
-19.94877	-13.21406	-20.72653	arable	POINT (-1.32533 59.96942)
-19.83694	-12.96664	-20.35517	arable	POINT (-1.32524 59.96942)
-20.04012	-13.01437	-20.20085	arable	POINT (-1.32515 59.96942)
-19.82112	-12.90779	-19.91867	arable	POINT (-1.32506 59.96942)
-21.18834	-13.36220	-20.41664	arable	POINT (-1.326228 59.96951)
-21.01876	-13.43764	-20.53008	arable	POINT (-1.326138 59.96951)
-20.69177	-13.43672	-20.50918	arable	POINT (-1.326048 59.96951)
-20.59111	-13.37023	-20.49506	arable	POINT (-1.325958 59.96951)
-20.43378	-13.29198	-20.56033	arable	POINT (-1.325868 59.96951)
-20.41118	-13.28322	-20.10121	arable	POINT (-1.325779 59.96951)
-20.36944	-13.34338	-20.23467	arable	POINT (-1.325689 59.96951)
-20.41864	-13.32802	-20.39319	arable	POINT (-1.325599 59.96951)
-20.31956	-13.35080	-20.43049	arable	POINT (-1.325509 59.96951)
-20.22145	-13.27553	-20.44900	arable	POINT (-1.325419 59.96951)
-20.24731	-13.20028	-20.29174	arable	POINT (-1.32533 59.96951)
-20.00794	-12.97447	-19.97122	arable	POINT (-1.32524 59.96951)
-20.06283	-12.89607	-19.79867	arable	POINT (-1.32515 59.96951)
-19.88746	-12.84259	-19.64511	arable	POINT (-1.32506 59.96951)
-20.85770	-13.32837	-20.40653	arable	POINT (-1.326228 59.9696)
-20.68879	-13.37659	-20.53830	arable	POINT (-1.326138 59.9696)
-20.39563	-13.33136	-20.66514	arable	POINT (-1.326048 59.9696)
-20.23980	-13.36294	-20.56236	arable	POINT (-1.325958 59.9696)
-20.24941	-13.28733	-20.61870	arable	POINT (-1.325868 59.9696)
-20.22343	-13.30267	-20.57296	arable	POINT (-1.325779 59.9696)
-20.22694	-13.30041	-20.45197	arable	POINT (-1.325689 59.9696)
-20.32863	-13.17939	-20.53450	arable	POINT (-1.325599 59.9696)
-20.14040	-13.19141	-20.57992	arable	POINT (-1.325509 59.9696)
-20.26803	-13.08915	-20.37702	arable	POINT (-1.325419 59.9696)
-20.11907	-12.96191	-20.25704	arable	POINT (-1.32533 59.9696)
-20.07743	-12.81463	-20.24693	arable	POINT (-1.32524 59.9696)
-19.99936	-12.67353	-19.96078	arable	POINT (-1.32515 59.9696)
-19.77626	-12.72000	-19.80505	arable	POINT (-1.32506 59.9696)
-20.23510	-13.34584	-20.27369	arable	POINT (-1.326228 59.96969)
-20.19379	-13.30125	-20.41511	arable	POINT (-1.326138 59.96969)
-20.16432	-13.21008	-20.75431	arable	POINT (-1.326048 59.96969)
-20.09491	-13.15004	-20.44683	arable	POINT (-1.325958 59.96969)
-19.97745	-13.12019	-20.62176	arable	POINT (-1.325868 59.96969)
-19.93991	-13.17738	-20.62721	arable	POINT (-1.325779 59.96969)
-20.01657	-13.15896	-20.55720	arable	POINT (-1.325689 59.96969)
-20.10550	-12.94943	-20.66532	arable	POINT (-1.325599 59.96969)
-20.04866	-12.98604	-20.53429	arable	POINT (-1.325509 59.96969)
-20.19087	-12.95102	-20.53580	arable	POINT (-1.325419 59.96969)
-20.07817	-12.83333	-20.42425	arable	POINT (-1.32533 59.96969)
-20.07566	-12.77687	-20.42853	arable	POINT (-1.32524 59.96969)
-19.92796	-12.66635	-20.17609	arable	POINT (-1.32515 59.96969)
-19.81896	-12.73411	-20.10009	arable	POINT (-1.32506 59.96969)
-20.03543	-13.38211	-20.28672	arable	POINT (-1.326228 59.96978)
-19.97436	-13.31415	-20.30629	arable	POINT (-1.326138 59.96978)
-20.07100	-13.33250	-20.44164	arable	POINT (-1.326048 59.96978)
-20.05391	-13.24441	-20.47105	arable	POINT (-1.325958 59.96978)
-20.06691	-13.14589	-20.65197	arable	POINT (-1.325868 59.96978)
-20.07415	-13.19373	-20.73987	arable	POINT (-1.325779 59.96978)
-20.19356	-13.08856	-20.66456	arable	POINT (-1.325689 59.96978)
-20.23397	-13.07550	-20.54150	arable	POINT (-1.325599 59.96978)
-20.27905	-13.00442	-20.41921	arable	POINT (-1.325509 59.96978)
-20.41404	-12.97102	-20.35579	arable	POINT (-1.325419 59.96978)
-20.30415	-12.94876	-20.26441	arable	POINT (-1.32533 59.96978)
-20.09733	-12.94789	-20.37518	arable	POINT (-1.32524 59.96978)
-19.93495	-12.84096	-20.21013	arable	POINT (-1.32515 59.96978)
-19.85271	-12.87528	-20.35620	arable	POINT (-1.32506 59.96978)

a <- as.data.frame(mySamples_sf_batch)
a  %>%  kbl() %>%
 kable_material(c("striped", "hover"), font_size = 12)%>%
  scroll_box(width = "100%", height = "200px")

X0_VH	X1_VV	X2_VH	name	geometry
-22.34926	-15.91189	-22.87572	unimproved	POINT (-1.345811 59.97238)
-22.21591	-15.88941	-22.73153	unimproved	POINT (-1.345721 59.97238)
-21.97966	-16.01451	-22.78574	unimproved	POINT (-1.345631 59.97238)
-21.89113	-16.02274	-22.80749	unimproved	POINT (-1.345541 59.97238)
-21.76509	-15.97965	-22.74236	unimproved	POINT (-1.345452 59.97238)
-21.80516	-16.13436	-22.81662	unimproved	POINT (-1.345362 59.97238)
-21.89018	-16.11657	-22.58584	unimproved	POINT (-1.345272 59.97238)
-21.96362	-16.09216	-22.49100	unimproved	POINT (-1.345182 59.97238)
-22.31687	-15.96285	-22.25296	unimproved	POINT (-1.345092 59.97238)
-22.26144	-15.83632	-21.99900	unimproved	POINT (-1.345003 59.97238)
-22.69117	-15.79825	-21.94824	unimproved	POINT (-1.344913 59.97238)
-22.86474	-15.62087	-21.95547	unimproved	POINT (-1.344823 59.97238)
-22.79661	-15.47186	-22.14489	unimproved	POINT (-1.344733 59.97238)
-23.09050	-15.27697	-22.12411	unimproved	POINT (-1.344643 59.97238)
-22.82591	-15.13440	-22.31815	unimproved	POINT (-1.344553 59.97238)
-22.71445	-14.87813	-22.21076	unimproved	POINT (-1.344464 59.97238)
-22.49065	-14.74864	-22.00986	unimproved	POINT (-1.344374 59.97238)
-22.07535	-14.73045	-22.03094	unimproved	POINT (-1.344284 59.97238)
-21.82844	-14.81384	-21.82220	unimproved	POINT (-1.344194 59.97238)
-21.71726	-14.91948	-21.73916	unimproved	POINT (-1.344104 59.97238)
-22.03812	-15.09589	-21.70903	unimproved	POINT (-1.344014 59.97238)
-22.04823	-15.15580	-21.67277	unimproved	POINT (-1.343925 59.97238)
-22.29961	-15.91687	-22.88622	unimproved	POINT (-1.345811 59.97247)
-22.16317	-15.89559	-22.56592	unimproved	POINT (-1.345721 59.97247)
-22.22898	-15.89353	-22.59414	unimproved	POINT (-1.345631 59.97247)
-22.04797	-15.91060	-22.71691	unimproved	POINT (-1.345541 59.97247)
-21.97182	-16.07256	-22.53795	unimproved	POINT (-1.345452 59.97247)
-21.97494	-16.11200	-22.66757	unimproved	POINT (-1.345362 59.97247)
-22.07728	-16.05436	-22.54982	unimproved	POINT (-1.345272 59.97247)
-22.11634	-15.97489	-22.93400	unimproved	POINT (-1.345182 59.97247)
-22.32816	-15.84062	-22.64100	unimproved	POINT (-1.345092 59.97247)
-22.27004	-15.76314	-22.29578	unimproved	POINT (-1.345003 59.97247)
-22.36537	-15.72136	-22.17396	unimproved	POINT (-1.344913 59.97247)
-22.29385	-15.50280	-21.85018	unimproved	POINT (-1.344823 59.97247)
-22.28096	-15.34687	-22.10721	unimproved	POINT (-1.344733 59.97247)
-22.51557	-15.17942	-22.08745	unimproved	POINT (-1.344643 59.97247)
-22.75653	-15.12734	-22.25224	unimproved	POINT (-1.344553 59.97247)
-22.87004	-15.07264	-22.23837	unimproved	POINT (-1.344464 59.97247)
-22.86395	-14.97103	-22.18748	unimproved	POINT (-1.344374 59.97247)
-22.69602	-14.94965	-22.17328	unimproved	POINT (-1.344284 59.97247)
-22.16186	-14.96355	-22.12732	unimproved	POINT (-1.344194 59.97247)
-21.97112	-15.00835	-22.09475	unimproved	POINT (-1.344104 59.97247)
-21.99987	-15.22618	-22.14131	unimproved	POINT (-1.344014 59.97247)
-21.92010	-15.26917	-22.13194	unimproved	POINT (-1.343925 59.97247)
-21.92852	-15.79740	-22.95455	unimproved	POINT (-1.345811 59.97256)
-21.90721	-15.80009	-22.37697	unimproved	POINT (-1.345721 59.97256)
-22.18528	-15.75753	-22.25433	unimproved	POINT (-1.345631 59.97256)
-22.20729	-15.77410	-22.26069	unimproved	POINT (-1.345541 59.97256)
-22.51404	-15.88464	-22.58650	unimproved	POINT (-1.345452 59.97256)
-22.53036	-15.86774	-22.94224	unimproved	POINT (-1.345362 59.97256)
-22.53977	-15.83087	-22.68218	unimproved	POINT (-1.345272 59.97256)
-22.52488	-15.79566	-22.98419	unimproved	POINT (-1.345182 59.97256)
-22.55057	-15.57543	-22.75748	unimproved	POINT (-1.345092 59.97256)
-22.30694	-15.50906	-22.25421	unimproved	POINT (-1.345003 59.97256)
-22.17115	-15.47788	-22.15552	unimproved	POINT (-1.344913 59.97256)
-21.95561	-15.27768	-21.92158	unimproved	POINT (-1.344823 59.97256)
-21.79101	-15.28003	-22.26744	unimproved	POINT (-1.344733 59.97256)
-22.06672	-15.15738	-22.35577	unimproved	POINT (-1.344643 59.97256)
-22.44712	-15.20579	-22.39809	unimproved	POINT (-1.344553 59.97256)
-22.61144	-15.23101	-22.71306	unimproved	POINT (-1.344464 59.97256)
-22.63982	-15.23158	-22.72155	unimproved	POINT (-1.344374 59.97256)
-22.54517	-15.29802	-22.72000	unimproved	POINT (-1.344284 59.97256)
-22.24314	-15.27920	-22.88174	unimproved	POINT (-1.344194 59.97256)
-21.91469	-15.30846	-22.75703	unimproved	POINT (-1.344104 59.97256)
-21.76489	-15.38044	-22.50904	unimproved	POINT (-1.344014 59.97256)
-21.78670	-15.40461	-22.20937	unimproved	POINT (-1.343925 59.97256)
-21.97779	-15.48170	-22.48384	unimproved	POINT (-1.345811 59.97265)
-21.95936	-15.54435	-22.17528	unimproved	POINT (-1.345721 59.97265)
-22.11548	-15.42745	-22.10846	unimproved	POINT (-1.345631 59.97265)
-22.25962	-15.42753	-22.15615	unimproved	POINT (-1.345541 59.97265)
-22.52362	-15.36528	-22.32281	unimproved	POINT (-1.345452 59.97265)
-22.59466	-15.32997	-22.58449	unimproved	POINT (-1.345362 59.97265)
-22.83997	-15.43620	-22.62401	unimproved	POINT (-1.345272 59.97265)
-22.76557	-15.48969	-22.79838	unimproved	POINT (-1.345182 59.97265)
-22.73486	-15.39360	-22.65099	unimproved	POINT (-1.345092 59.97265)
-22.55798	-15.35668	-22.37452	unimproved	POINT (-1.345003 59.97265)
-22.18621	-15.32172	-22.40649	unimproved	POINT (-1.344913 59.97265)
-21.99974	-15.31439	-22.29078	unimproved	POINT (-1.344823 59.97265)
-21.69493	-15.32975	-22.64503	unimproved	POINT (-1.344733 59.97265)
-21.75851	-15.21288	-22.80644	unimproved	POINT (-1.344643 59.97265)
-21.91536	-15.18926	-22.96097	unimproved	POINT (-1.344553 59.97265)
-22.07930	-15.27092	-23.08847	unimproved	POINT (-1.344464 59.97265)
-22.09493	-15.28085	-22.92289	unimproved	POINT (-1.344374 59.97265)
-22.11447	-15.28952	-22.99562	unimproved	POINT (-1.344284 59.97265)
-22.03848	-15.30274	-22.91792	unimproved	POINT (-1.344194 59.97265)
-21.88705	-15.22280	-22.79154	unimproved	POINT (-1.344104 59.97265)
-21.80391	-15.30598	-22.42220	unimproved	POINT (-1.344014 59.97265)
-21.83127	-15.37444	-22.14390	unimproved	POINT (-1.343925 59.97265)
-21.86570	-15.20233	-22.66075	unimproved	POINT (-1.345811 59.97274)
-21.78771	-15.18262	-22.71818	unimproved	POINT (-1.345721 59.97274)
-21.76648	-15.12879	-22.45075	unimproved	POINT (-1.345631 59.97274)
-21.89388	-15.09136	-22.29077	unimproved	POINT (-1.345541 59.97274)
-22.04722	-15.05048	-22.04886	unimproved	POINT (-1.345452 59.97274)
-22.15751	-15.10574	-22.21241	unimproved	POINT (-1.345362 59.97274)
-22.31850	-15.24818	-22.40979	unimproved	POINT (-1.345272 59.97274)
-22.42901	-15.34082	-22.66044	unimproved	POINT (-1.345182 59.97274)
-22.46221	-15.33020	-22.59297	unimproved	POINT (-1.345092 59.97274)
-22.29596	-15.34494	-22.58145	unimproved	POINT (-1.345003 59.97274)
-22.21360	-15.32532	-22.80264	unimproved	POINT (-1.344913 59.97274)
-22.05430	-15.32682	-22.72462	unimproved	POINT (-1.344823 59.97274)
-21.85202	-15.35268	-23.04824	unimproved	POINT (-1.344733 59.97274)
-21.57052	-15.21132	-23.00643	unimproved	POINT (-1.344643 59.97274)
-21.37291	-15.12446	-22.96794	unimproved	POINT (-1.344553 59.97274)
-21.55231	-15.17566	-22.82254	unimproved	POINT (-1.344464 59.97274)
-21.76842	-15.11345	-22.48793	unimproved	POINT (-1.344374 59.97274)
-21.97525	-15.02065	-22.34484	unimproved	POINT (-1.344284 59.97274)
-22.10091	-15.02157	-21.95459	unimproved	POINT (-1.344194 59.97274)
-22.07181	-14.95599	-22.00109	unimproved	POINT (-1.344104 59.97274)
-22.07840	-14.92686	-22.06321	unimproved	POINT (-1.344014 59.97274)
-22.19837	-15.02488	-22.00927	unimproved	POINT (-1.343925 59.97274)
-21.76553	-15.30813	-22.80079	unimproved	POINT (-1.345811 59.97283)
-21.62665	-15.27709	-22.78941	unimproved	POINT (-1.345721 59.97283)
-21.56074	-15.31464	-22.64466	unimproved	POINT (-1.345631 59.97283)
-21.60453	-15.24232	-22.46348	unimproved	POINT (-1.345541 59.97283)
-21.49362	-15.14964	-22.22831	unimproved	POINT (-1.345452 59.97283)
-21.67327	-15.07540	-22.33174	unimproved	POINT (-1.345362 59.97283)
-21.85949	-15.11551	-22.24978	unimproved	POINT (-1.345272 59.97283)
-22.04067	-15.16645	-22.46252	unimproved	POINT (-1.345182 59.97283)
-22.40039	-15.22902	-22.72918	unimproved	POINT (-1.345092 59.97283)
-22.34913	-15.34380	-22.73451	unimproved	POINT (-1.345003 59.97283)
-22.25212	-15.35941	-23.08429	unimproved	POINT (-1.344913 59.97283)
-22.09622	-15.39595	-23.11525	unimproved	POINT (-1.344823 59.97283)
-21.81326	-15.36927	-22.94700	unimproved	POINT (-1.344733 59.97283)
-21.42869	-15.25391	-22.80808	unimproved	POINT (-1.344643 59.97283)
-21.32290	-15.18844	-22.63324	unimproved	POINT (-1.344553 59.97283)
-21.69392	-15.11559	-22.45803	unimproved	POINT (-1.344464 59.97283)
-22.13406	-15.02552	-22.11996	unimproved	POINT (-1.344374 59.97283)
-22.30202	-14.76567	-21.86035	unimproved	POINT (-1.344284 59.97283)
-22.53209	-14.63251	-21.46585	unimproved	POINT (-1.344194 59.97283)
-22.40516	-14.57303	-21.56973	unimproved	POINT (-1.344104 59.97283)
-22.38931	-14.48128	-21.77187	unimproved	POINT (-1.344014 59.97283)
-22.39138	-14.71058	-21.91156	unimproved	POINT (-1.343925 59.97283)
-22.26588	-15.65654	-22.81160	unimproved	POINT (-1.345811 59.97292)
-22.19764	-15.71238	-23.09331	unimproved	POINT (-1.345721 59.97292)
-22.00024	-15.70143	-23.17772	unimproved	POINT (-1.345631 59.97292)
-21.94204	-15.66371	-23.06460	unimproved	POINT (-1.345541 59.97292)
-21.83799	-15.51219	-23.00098	unimproved	POINT (-1.345452 59.97292)
-21.85356	-15.29634	-22.80835	unimproved	POINT (-1.345362 59.97292)
-21.87769	-15.25407	-22.36239	unimproved	POINT (-1.345272 59.97292)
-21.95982	-15.18768	-22.43488	unimproved	POINT (-1.345182 59.97292)
-22.04289	-15.31008	-22.67168	unimproved	POINT (-1.345092 59.97292)
-22.15053	-15.40118	-22.82156	unimproved	POINT (-1.345003 59.97292)
-22.25697	-15.44189	-22.89163	unimproved	POINT (-1.344913 59.97292)
-22.12453	-15.46480	-23.03866	unimproved	POINT (-1.344823 59.97292)
-21.94612	-15.40550	-22.72031	unimproved	POINT (-1.344733 59.97292)
-21.89403	-15.26791	-22.73055	unimproved	POINT (-1.344643 59.97292)
-21.82125	-15.16462	-22.48966	unimproved	POINT (-1.344553 59.97292)
-22.10257	-15.01990	-22.08998	unimproved	POINT (-1.344464 59.97292)
-22.37544	-14.85929	-21.96502	unimproved	POINT (-1.344374 59.97292)
-22.38999	-14.70191	-21.68476	unimproved	POINT (-1.344284 59.97292)
-22.50571	-14.60607	-21.42378	unimproved	POINT (-1.344194 59.97292)
-22.51082	-14.53075	-21.59657	unimproved	POINT (-1.344104 59.97292)
-22.59015	-14.43234	-21.84483	unimproved	POINT (-1.344014 59.97292)
-22.59395	-14.58702	-21.75558	unimproved	POINT (-1.343925 59.97292)
-22.59976	-15.82260	-22.92425	unimproved	POINT (-1.345811 59.97301)
-22.54212	-15.78813	-23.02629	unimproved	POINT (-1.345721 59.97301)
-22.47538	-15.89120	-22.99235	unimproved	POINT (-1.345631 59.97301)
-22.44055	-15.85456	-23.01211	unimproved	POINT (-1.345541 59.97301)
-22.37847	-15.73307	-22.90147	unimproved	POINT (-1.345452 59.97301)
-22.25150	-15.59374	-22.43795	unimproved	POINT (-1.345362 59.97301)
-21.83943	-15.52617	-22.24057	unimproved	POINT (-1.345272 59.97301)
-21.79085	-15.46399	-22.22201	unimproved	POINT (-1.345182 59.97301)
-21.78509	-15.47956	-22.28311	unimproved	POINT (-1.345092 59.97301)
-21.99629	-15.50004	-22.58030	unimproved	POINT (-1.345003 59.97301)
-22.07003	-15.44787	-22.76644	unimproved	POINT (-1.344913 59.97301)
-22.24255	-15.35868	-23.06027	unimproved	POINT (-1.344823 59.97301)
-22.32707	-15.36473	-22.70823	unimproved	POINT (-1.344733 59.97301)
-22.40869	-15.19686	-22.41702	unimproved	POINT (-1.344643 59.97301)
-22.38243	-15.03614	-22.30198	unimproved	POINT (-1.344553 59.97301)
-22.17149	-14.95536	-22.15409	unimproved	POINT (-1.344464 59.97301)
-22.29183	-14.83017	-22.06424	unimproved	POINT (-1.344374 59.97301)
-22.15220	-14.91144	-22.03322	unimproved	POINT (-1.344284 59.97301)
-22.31998	-14.90031	-22.02175	unimproved	POINT (-1.344194 59.97301)
-22.40842	-14.75952	-21.81305	unimproved	POINT (-1.344104 59.97301)
-22.39792	-14.69019	-21.83295	unimproved	POINT (-1.344014 59.97301)
-22.44975	-14.62622	-21.76186	unimproved	POINT (-1.343925 59.97301)
-22.42614	-15.78419	-23.11120	unimproved	POINT (-1.345811 59.9731)
-22.38357	-15.71853	-23.06984	unimproved	POINT (-1.345721 59.9731)
-22.34121	-15.72388	-22.79326	unimproved	POINT (-1.345631 59.9731)
-22.41764	-15.70565	-22.51386	unimproved	POINT (-1.345541 59.9731)
-22.35609	-15.72229	-22.44312	unimproved	POINT (-1.345452 59.9731)
-22.28237	-15.65711	-22.35886	unimproved	POINT (-1.345362 59.9731)
-21.99147	-15.56860	-22.27400	unimproved	POINT (-1.345272 59.9731)
-22.07032	-15.45841	-22.17037	unimproved	POINT (-1.345182 59.9731)
-21.89464	-15.39335	-22.04005	unimproved	POINT (-1.345092 59.9731)
-21.99809	-15.37769	-22.15215	unimproved	POINT (-1.345003 59.9731)
-22.15370	-15.45542	-22.32369	unimproved	POINT (-1.344913 59.9731)
-22.28458	-15.37987	-22.60357	unimproved	POINT (-1.344823 59.9731)
-22.34969	-15.28493	-22.55929	unimproved	POINT (-1.344733 59.9731)
-22.33060	-15.21955	-22.36315	unimproved	POINT (-1.344643 59.9731)
-22.18611	-14.99822	-22.35174	unimproved	POINT (-1.344553 59.9731)
-21.97796	-14.85605	-22.29448	unimproved	POINT (-1.344464 59.9731)
-21.97427	-14.81289	-22.27880	unimproved	POINT (-1.344374 59.9731)
-21.98473	-14.83447	-22.40398	unimproved	POINT (-1.344284 59.9731)
-22.13085	-14.88455	-22.22783	unimproved	POINT (-1.344194 59.9731)
-22.38240	-14.85527	-22.01012	unimproved	POINT (-1.344104 59.9731)
-22.53206	-14.79811	-22.04463	unimproved	POINT (-1.344014 59.9731)
-22.53541	-14.76501	-21.76268	unimproved	POINT (-1.343925 59.9731)
-22.23757	-15.60706	-23.14522	unimproved	POINT (-1.345811 59.97319)
-22.07327	-15.55352	-23.07737	unimproved	POINT (-1.345721 59.97319)
-22.16035	-15.51023	-22.89436	unimproved	POINT (-1.345631 59.97319)
-22.30915	-15.59514	-22.44612	unimproved	POINT (-1.345541 59.97319)
-22.33782	-15.56478	-22.18060	unimproved	POINT (-1.345452 59.97319)
-22.30500	-15.61097	-22.17947	unimproved	POINT (-1.345362 59.97319)
-22.09543	-15.56551	-22.12238	unimproved	POINT (-1.345272 59.97319)
-22.02195	-15.51946	-22.14096	unimproved	POINT (-1.345182 59.97319)
-21.96520	-15.45304	-22.11039	unimproved	POINT (-1.345092 59.97319)
-21.98657	-15.37232	-22.02296	unimproved	POINT (-1.345003 59.97319)
-22.16159	-15.29753	-22.23248	unimproved	POINT (-1.344913 59.97319)
-22.21102	-15.20583	-22.38738	unimproved	POINT (-1.344823 59.97319)
-22.11585	-15.05681	-22.38600	unimproved	POINT (-1.344733 59.97319)
-22.14490	-14.93161	-22.32081	unimproved	POINT (-1.344643 59.97319)
-21.95658	-14.92937	-22.25154	unimproved	POINT (-1.344553 59.97319)
-21.85696	-14.82848	-22.24205	unimproved	POINT (-1.344464 59.97319)
-21.82417	-14.81088	-22.42299	unimproved	POINT (-1.344374 59.97319)
-22.02679	-14.78798	-22.23210	unimproved	POINT (-1.344284 59.97319)
-22.20547	-14.78178	-22.26198	unimproved	POINT (-1.344194 59.97319)
-22.42942	-14.79622	-21.87109	unimproved	POINT (-1.344104 59.97319)
-22.77939	-14.85581	-21.74714	unimproved	POINT (-1.344014 59.97319)
-22.49362	-14.91334	-21.68610	unimproved	POINT (-1.343925 59.97319)
-22.21044	-15.52502	-22.83231	unimproved	POINT (-1.345811 59.97328)
-21.91846	-15.46858	-22.63349	unimproved	POINT (-1.345721 59.97328)
-22.02031	-15.47524	-22.57350	unimproved	POINT (-1.345631 59.97328)
-22.05058	-15.53830	-22.36777	unimproved	POINT (-1.345541 59.97328)
-22.03050	-15.53598	-22.33141	unimproved	POINT (-1.345452 59.97328)
-22.16771	-15.60631	-22.27188	unimproved	POINT (-1.345362 59.97328)
-22.15777	-15.42857	-22.28691	unimproved	POINT (-1.345272 59.97328)
-22.06865	-15.35641	-22.49895	unimproved	POINT (-1.345182 59.97328)
-22.05546	-15.36470	-22.57702	unimproved	POINT (-1.345092 59.97328)
-21.96615	-15.21841	-22.71653	unimproved	POINT (-1.345003 59.97328)
-21.76507	-15.11049	-22.88488	unimproved	POINT (-1.344913 59.97328)
-21.81775	-15.00730	-22.79723	unimproved	POINT (-1.344823 59.97328)
-21.91159	-15.01587	-22.74424	unimproved	POINT (-1.344733 59.97328)
-21.89335	-15.00588	-22.64776	unimproved	POINT (-1.344643 59.97328)
-21.90457	-15.06631	-22.50162	unimproved	POINT (-1.344553 59.97328)
-22.02905	-15.12495	-22.53719	unimproved	POINT (-1.344464 59.97328)
-22.15016	-15.07007	-22.65560	unimproved	POINT (-1.344374 59.97328)
-22.31467	-15.03613	-22.39747	unimproved	POINT (-1.344284 59.97328)
-22.39164	-14.96686	-22.34867	unimproved	POINT (-1.344194 59.97328)
-22.28125	-14.95092	-21.89616	unimproved	POINT (-1.344104 59.97328)
-22.23808	-14.94860	-21.67926	unimproved	POINT (-1.344014 59.97328)
-22.02255	-14.90802	-21.45511	unimproved	POINT (-1.343925 59.97328)
-22.28496	-15.44309	-22.46307	unimproved	POINT (-1.345811 59.97337)
-22.18744	-15.49726	-22.41188	unimproved	POINT (-1.345721 59.97337)
-21.99264	-15.48308	-22.46228	unimproved	POINT (-1.345631 59.97337)
-21.88923	-15.54486	-22.44620	unimproved	POINT (-1.345541 59.97337)
-21.92425	-15.52790	-22.35361	unimproved	POINT (-1.345452 59.97337)
-21.83892	-15.54099	-22.24788	unimproved	POINT (-1.345362 59.97337)
-22.14189	-15.51085	-22.29165	unimproved	POINT (-1.345272 59.97337)
-22.12796	-15.49753	-22.41546	unimproved	POINT (-1.345182 59.97337)
-22.05475	-15.46135	-22.73773	unimproved	POINT (-1.345092 59.97337)
-22.01917	-15.40319	-22.88660	unimproved	POINT (-1.345003 59.97337)
-21.87995	-15.25486	-22.92831	unimproved	POINT (-1.344913 59.97337)
-21.79203	-15.18836	-23.04165	unimproved	POINT (-1.344823 59.97337)
-22.05940	-15.30185	-22.97576	unimproved	POINT (-1.344733 59.97337)
-22.10435	-15.30093	-22.78901	unimproved	POINT (-1.344643 59.97337)
-22.43526	-15.26185	-22.87282	unimproved	POINT (-1.344553 59.97337)
-22.39315	-15.28386	-22.74946	unimproved	POINT (-1.344464 59.97337)
-22.31948	-15.23491	-22.64637	unimproved	POINT (-1.344374 59.97337)
-22.36651	-15.18727	-22.48144	unimproved	POINT (-1.344284 59.97337)
-22.03854	-15.21850	-22.09712	unimproved	POINT (-1.344194 59.97337)
-22.14178	-15.18509	-21.73988	unimproved	POINT (-1.344104 59.97337)
-22.22406	-15.13094	-21.54485	unimproved	POINT (-1.344014 59.97337)
-22.31737	-14.97945	-21.23326	unimproved	POINT (-1.343925 59.97337)

mySamples_sf_batch$X1_VV

##   [1] -15.91189 -15.88941 -16.01451 -16.02274 -15.97965 -16.13436 -16.11657
##   [8] -16.09216 -15.96285 -15.83632 -15.79825 -15.62087 -15.47186 -15.27697
##  [15] -15.13440 -14.87813 -14.74864 -14.73045 -14.81384 -14.91948 -15.09589
##  [22] -15.15580 -15.91687 -15.89559 -15.89353 -15.91060 -16.07256 -16.11200
##  [29] -16.05436 -15.97489 -15.84062 -15.76314 -15.72136 -15.50280 -15.34687
##  [36] -15.17942 -15.12734 -15.07264 -14.97103 -14.94965 -14.96355 -15.00835
##  [43] -15.22618 -15.26917 -15.79740 -15.80009 -15.75753 -15.77410 -15.88464
##  [50] -15.86774 -15.83087 -15.79566 -15.57543 -15.50906 -15.47788 -15.27768
##  [57] -15.28003 -15.15738 -15.20579 -15.23101 -15.23158 -15.29802 -15.27920
##  [64] -15.30846 -15.38044 -15.40461 -15.48170 -15.54435 -15.42745 -15.42753
##  [71] -15.36528 -15.32997 -15.43620 -15.48969 -15.39360 -15.35668 -15.32172
##  [78] -15.31439 -15.32975 -15.21288 -15.18926 -15.27092 -15.28085 -15.28952
##  [85] -15.30274 -15.22280 -15.30598 -15.37444 -15.20233 -15.18262 -15.12879
##  [92] -15.09136 -15.05048 -15.10574 -15.24818 -15.34082 -15.33020 -15.34494
##  [99] -15.32532 -15.32682 -15.35268 -15.21132 -15.12446 -15.17566 -15.11345
## [106] -15.02065 -15.02157 -14.95599 -14.92686 -15.02488 -15.30813 -15.27709
## [113] -15.31464 -15.24232 -15.14964 -15.07540 -15.11551 -15.16645 -15.22902
## [120] -15.34380 -15.35941 -15.39595 -15.36927 -15.25391 -15.18844 -15.11559
## [127] -15.02552 -14.76567 -14.63251 -14.57303 -14.48128 -14.71058 -15.65654
## [134] -15.71238 -15.70143 -15.66371 -15.51219 -15.29634 -15.25407 -15.18768
## [141] -15.31008 -15.40118 -15.44189 -15.46480 -15.40550 -15.26791 -15.16462
## [148] -15.01990 -14.85929 -14.70191 -14.60607 -14.53075 -14.43234 -14.58702
## [155] -15.82260 -15.78813 -15.89120 -15.85456 -15.73307 -15.59374 -15.52617
## [162] -15.46399 -15.47956 -15.50004 -15.44787 -15.35868 -15.36473 -15.19686
## [169] -15.03614 -14.95536 -14.83017 -14.91144 -14.90031 -14.75952 -14.69019
## [176] -14.62622 -15.78419 -15.71853 -15.72388 -15.70565 -15.72229 -15.65711
## [183] -15.56860 -15.45841 -15.39335 -15.37769 -15.45542 -15.37987 -15.28493
## [190] -15.21955 -14.99822 -14.85605 -14.81289 -14.83447 -14.88455 -14.85527
## [197] -14.79811 -14.76501 -15.60706 -15.55352 -15.51023 -15.59514 -15.56478
## [204] -15.61097 -15.56551 -15.51946 -15.45304 -15.37232 -15.29753 -15.20583
## [211] -15.05681 -14.93161 -14.92937 -14.82848 -14.81088 -14.78798 -14.78178
## [218] -14.79622 -14.85581 -14.91334 -15.52502 -15.46858 -15.47524 -15.53830
## [225] -15.53598 -15.60631 -15.42857 -15.35641 -15.36470 -15.21841 -15.11049
## [232] -15.00730 -15.01587 -15.00588 -15.06631 -15.12495 -15.07007 -15.03613
## [239] -14.96686 -14.95092 -14.94860 -14.90802 -15.44309 -15.49726 -15.48308
## [246] -15.54486 -15.52790 -15.54099 -15.51085 -15.49753 -15.46135 -15.40319
## [253] -15.25486 -15.18836 -15.30185 -15.30093 -15.26185 -15.28386 -15.23491
## [260] -15.18727 -15.21850 -15.18509 -15.13094 -14.97945

Global Surface Water

library(rgee)
# ee_reattach() # reattach ee as a reserved word

ee_Initialize()

## -- rgee 1.0.9 --------------------------------------- earthengine-api 0.1.259 -- 
##  v email: not_defined
##  v Initializing Google Earth Engine:
 v Initializing Google Earth Engine:  DONE!
## 
 v Earth Engine user: users/tarredwall 
## --------------------------------------------------------------------------------

###############################
# Asset List
###############################

gsw <- ee$Image("JRC/GSW1_1/GlobalSurfaceWater")
occurrence <- gsw$select("occurrence")

###############################
# Constants
###############################

VIS_OCCURRENCE <- list(
  min = 0,
  max = 100,
  palette = c("red", "blue")
)

VIS_WATER_MASK <- list(
  palette = c("white", "black")
)

###############################
# Calculations
###############################

# Create a water mask layer, and set the image mask so that non-water areas
# are opaque.
water_mask <- occurrence$gt(90)$selfMask()

###############################
# Initialize Map Location
###############################

# Uncomment one of the following statements to center the map.
# Map$setCenter(-90.162, 29.8597, 10)   # New Orleans, USA
# Map$setCenter(-114.9774, 31.9254, 10) # Mouth of the Colorado River, Mexico
# Map$setCenter(-111.1871, 37.0963, 11) # Lake Powell, USA
# Map$setCenter(149.412, -35.0789, 11)  # Lake George, Australia
# Map$setCenter(105.26, 11.2134, 9)     # Mekong River Basin, SouthEast Asia
# Map$setCenter(90.6743, 22.7382, 10)   # Meghna River, Bangladesh
# Map$setCenter(81.2714, 16.5079, 11)   # Godavari River Basin Irrigation Project, India
# Map$setCenter(14.7035, 52.0985, 12)   # River Oder, Germany & Poland
# Map$setCenter(-59.1696, -33.8111, 9)  # Buenos Aires, Argentina
Map$setCenter(-74.4557, -8.4289, 11)  # Ucayali River, Peru

###############################
# Map Layers
###############################
Map$addLayer(occurrence$updateMask(occurrence$divide(100)), VIS_OCCURRENCE, "Water Occurrence (1984-2018)") +
Map$addLayer(water_mask, VIS_WATER_MASK, "90% occurrence water mask", FALSE)

2_water_occurrence_change_intensity.R

library(rgee)
# ee_reattach() # reattach ee as a reserved word

ee_Initialize()

## -- rgee 1.0.9 --------------------------------------- earthengine-api 0.1.259 -- 
##  v email: not_defined
##  v Initializing Google Earth Engine:
 v Initializing Google Earth Engine:  DONE!
## 
 v Earth Engine user: users/tarredwall 
## --------------------------------------------------------------------------------

###############################
# Asset List
###############################

gsw = ee$Image('JRC/GSW1_1/GlobalSurfaceWater')
occurrence = gsw$select('occurrence')
change = gsw$select("change_abs")
roi = ee$Geometry$Polygon(
  list(
    c(-74.17213, -8.65569),
    c(-74.17419, -8.39222),
    c(-74.38362, -8.36980),
    c(-74.43031, -8.61293)
  )
)

###############################
# Constants
###############################

VIS_OCCURRENCE = list(
  min = 0,
  max = 100,
  palette = c('red', 'blue')
)

VIS_CHANGE = list(
  min = -50,
  max = 50,
  palette = c('red', 'black', 'limegreen')
)

VIS_WATER_MASK = list(
  palette = c('white', 'black')
)

###############################
# Calculations
###############################

# Create a water mask layer, and set the image mask so that non-water areas are transparent.
water_mask = occurrence$gt(90)$mask(1)

###############################
# Initialize Map Location
###############################

# Uncomment one of the following statements to center the map on
# a particular location.
# Map$setCenter(-90.162, 29.8597, 10)   # New Orleans, USA
# Map$setCenter(-114.9774, 31.9254, 10) # Mouth of the Colorado River, Mexico
# Map$setCenter(-111.1871, 37.0963, 11) # Lake Powell, USA
# Map$setCenter(149.412, -35.0789, 11)  # Lake George, Australia
# Map$setCenter(105.26, 11.2134, 9)     # Mekong River Basin, SouthEast Asia
# Map$setCenter(90.6743, 22.7382, 10)   # Meghna River, Bangladesh
# Map$setCenter(81.2714, 16.5079, 11)   # Godavari River Basin Irrigation Project, India
# Map$setCenter(14.7035, 52.0985, 12)   # River Oder, Germany & Poland
# Map$setCenter(-59.1696, -33.8111, 9)  # Buenos Aires, Argentina\
Map$setCenter(-74.4557, -8.4289, 11)  # Ucayali River, Peru

###############################
# Map Layers
###############################

Map$addLayer(water_mask, VIS_WATER_MASK, '90% occurrence water mask', FALSE) +
Map$addLayer(occurrence$updateMask(occurrence$divide(100)),  VIS_OCCURRENCE, "Water Occurrence (1984-2015)") +
Map$addLayer(change, VIS_CHANGE,'occurrence change intensity')

Keiko/fire_australia.R

library(rgee)
# ee_reattach() # reattach ee as a reserved word

ee_Initialize()

## -- rgee 1.0.9 --------------------------------------- earthengine-api 0.1.259 -- 
##  v email: not_defined
##  v Initializing Google Earth Engine:
 v Initializing Google Earth Engine:  DONE!
## 
 v Earth Engine user: users/tarredwall 
## --------------------------------------------------------------------------------

# Credits to: Keiko Nomura, Senior Analyst, Space Intelligence Ltd
# Source: https://medium.com/google-earth/10-tips-for-becoming-an-earth-engine-expert-b11aad9e598b
# GEE JS: https://code.earthengine.google.com/?scriptPath=users%2Fnkeikon%2Fmedium%3Afire_australia

geometry <- ee$Geometry$Polygon(
  list(
    c(153.02512376008724, -28.052192238512877),
    c(153.02512376008724, -28.702237664294238),
    c(153.65683762727474, -28.702237664294238),
    c(153.65683762727474, -28.052192238512877)
  )
)

Map$centerObject(ee$FeatureCollection(geometry), zoom = 10)

## NOTE: Center obtained from the first element.

# Use clear images from May and Dec 2019
imageMay <- ee$Image("COPERNICUS/S2_SR/20190506T235259_20190506T235253_T56JNP")
imageDec <- ee$Image("COPERNICUS/S2_SR/20191202T235239_20191202T235239_T56JNP")

Map$addLayer(
  eeObject = imageMay,
  visParams = list(
    bands = c("B4", "B3", "B2"),
    min = 0,
    max = 1800
  ),
  name = "May 2019 (True colours)"
) +
  Map$addLayer(
    eeObject = imageDec,
    visParams = list(
      bands = c("B4", "B3", "B2"),
      min = 0,
      max = 1800
    ),
    name = "Dec 2019 (True colours)"
  )

# Compute NDVI and use grey colour for areas with NDVI < 0.8 in May 2019
NDVI <- imageMay$normalizedDifference(c("B8", "B4"))$rename("NDVI")
grey <- imageMay$mask(NDVI$select("NDVI")$lt(0.8))

Map$addLayer(
  eeObject = grey,
  visParams = list(
    bands = c("B3", "B3", "B3"),
    min = 0,
    max = 1800,
    gamma = 1.5
  ),
  name = "grey (base)"
)

# Export as mosaic. Alternatively you can also use blend().
mosaicDec <- ee$ImageCollection(
  c(
    imageDec$visualize(
      bands = c("B4", "B3", "B2"),
      min = 0,
      max = 1800
    ),
    grey$visualize(
      bands = c("B3", "B3", "B3"),
      min = 0,
      max = 1800
    )
  )
)$mosaic()

mosaicMay <- ee$ImageCollection(c(
  imageMay$visualize(
    bands = c("B4", "B3", "B2"),
    min = 0,
    max = 1800
  ),
  grey$visualize(
    bands = c("B3", "B3", "B3"),
    min = 0,
    max = 1800
  )
))$mosaic()

# ee_image_to_drive(
#   image = mosaicMay,
#   description = 'May',
#   region = geometry,
#   crs = 'EPSG:3857',
#   scale = 10
# )

# ee_image_to_drive(
#   image = mosaicDec,
#   description = 'Dec',
#   region = geometry,
#   crs = 'EPSG:3857',
#   scale = 10
# )

# ============ #
#  Topography  #
# ============ #

# Add topography by computing a hillshade using the terrain algorithms
elev <- ee$Image("USGS/SRTMGL1_003")
shadeAll <- ee$Terrain$hillshade(elev)
shade <- shadeAll$mask(elev$gt(0)) # mask the sea

mayTR <- ee$ImageCollection(c(
  imageMay$visualize(
    bands = c("B4", "B3", "B2"),
    min = 0,
    max = 1800
  ),
  shade$visualize(
    bands = c("hillshade", "hillshade", "hillshade"),
    opacity = 0.2
  )
))$mosaic()

highVeg <- NDVI$gte(0.8)$visualize(
  min = 0,
  max = 1
)

Map$addLayer(mayTR$mask(highVeg), list(gamma = 0.8), "May (with topography)")

# Convert the visualized elevation to HSV, first converting to [0, 1] data.
hsv <- mayTR$divide(255)$rgbToHsv()
# Select only the hue and saturation bands.
hs <- hsv$select(0, 1)
# Convert the hillshade to [0, 1] data, as expected by the HSV algorithm.
v <- shade$divide(255)
# Create a visualization image by converting back to RGB from HSV.
# Note the cast to byte in order to export the image correctly.
rgb <- hs$addBands(v)$hsvToRgb()$multiply(255)$byte()

Map$addLayer(rgb$mask(highVeg), list(gamma = 0.5), "May (topography visualised)")

# Export the image
mayTRMosaic <- ee$ImageCollection(c(
  rgb$mask(highVeg)$visualize(gamma = 0.5),
  grey$visualize(
    bands = c("B3", "B3", "B3"),
    min = 0,
    max = 1800
  )
))$mosaic()

# ee_image_to_drive(
#   image = mayTRMosaic,
#   description = 'MayTerrain',
#   region = geometry,
#   crs = 'EPSG:3857',
#   scale = 10
# )

decTR <- ee$ImageCollection(c(
  imageDec$visualize(
    bands = c("B4", "B3", "B2"),
    min = 0,
    max = 1800
  ),
  shade$visualize(
    bands = c("hillshade", "hillshade", "hillshade"),
    opacity = 0.2
  )
))$mosaic()

Map$addLayer(decTR$mask(highVeg), list(gamma = 0.8), "Dec (with topography)")

# Convert the visualized elevation to HSV, first converting to [0, 1] data.
hsv <- decTR$divide(255)$rgbToHsv()
# Select only the hue and saturation bands.
hs <- hsv$select(0, 1)
# Convert the hillshade to [0, 1] data, as expected by the HSV algorithm.
v <- shade$divide(255)
# Create a visualization image by converting back to RGB from HSV.
# Note the cast to byte in order to export the image correctly.
rgb <- hs$addBands(v)$hsvToRgb()$multiply(255)$byte()

Map$addLayer(rgb$mask(highVeg), list(gamma = 0.5), "Dec (topography visualised)")

# Export the image
decTRMosaic <- ee$ImageCollection(c(
  rgb$mask(highVeg)$visualize(
    gamma = 0.5
  ),
  grey$visualize(
    bands = c("B3", "B3", "B3"),
    min = 0,
    max = 1800
  )
))$mosaic()

# ee_image_to_drive(
#   image = decTRMosaic,
#   description = 'DecTerrain',
#   region = geometry,
#   crs = 'EPSG:3857',
#   scale = 10
# )

Algorithms/sentinel-1_filtering.R

library(rgee)
# ee_reattach() # reattach ee as a reserved word

ee_Initialize()

## -- rgee 1.0.9 --------------------------------------- earthengine-api 0.1.259 -- 
##  v email: not_defined
##  v Initializing Google Earth Engine:
 v Initializing Google Earth Engine:  DONE!
## 
 v Earth Engine user: users/tarredwall 
## --------------------------------------------------------------------------------

# Load the Sentinel-1 ImageCollection.
sentinel1 <- ee$ImageCollection("COPERNICUS/S1_GRD")$
  filterBounds(ee$Geometry$Point(-122.37383, 37.6193))

# Filter by metadata properties.
vh <- sentinel1$
  filter(ee$Filter$listContains("transmitterReceiverPolarisation", "VV"))$
  filter(ee$Filter$listContains("transmitterReceiverPolarisation", "VH"))$
  filter(ee$Filter$eq("instrumentMode", "IW"))

# Filter to get images from different look angles.
vhAscending <- vh$filter(ee$Filter$eq("orbitProperties_pass", "ASCENDING"))
vhDescending <- vh$filter(ee$Filter$eq("orbitProperties_pass", "DESCENDING"))

# Create a composite from means at different polarizations and look angles.
VH_Ascending_mean <- vhAscending$select("VH")$mean()
VV_Ascending_Descending_mean <- vhAscending$select("VV") %>%
  ee$ImageCollection$merge(vhDescending$select("VV")) %>%
  ee$ImageCollection$mean()
VH_Descending_mean <- vhDescending$select("VH")$mean()

composite <- ee$Image$cat(list(
  VH_Ascending_mean,
  VV_Ascending_Descending_mean,
  VH_Descending_mean
))$focal_median()

# Display as a composite of polarization and backscattering characteristics.
Map$setCenter(-122.37383, 37.6193, 10)
Map$addLayer(
  composite,
  list(min = c(-25, -20, -25), max = c(0, 10, 0)),
  "composite"
)

GPWv411: Basic Demographic Characteristics

The Gridded Population of World Version 4 (GPWv4), Revisión 11 modela la distribución de la población humana global para los años 2000, 2005, 2010, 2015 y 2020 en celdas de cuadrícula de 30 segundos de arco (aproximadamente 1 km). La población se distribuye en celdas utilizando la asignación proporcional de población de las unidades censales y administrativas. Los datos de entrada de población se recopilan con la resolución espacial más detallada disponible a partir de los resultados de la ronda de censos de 2010, que tuvo lugar entre 2005 y 2014. Los datos de entrada se extrapolan para producir estimaciones de población para cada año modelado.

Estas cuadrículas de conteo de población contienen estimaciones de población, por edad y sexo, por celda de cuadrícula de 30 segundos de arco consistente con los censos nacionales y los registros de población. Hay una imagen para cada categoría de edad y sexo modelada basada en la ronda del censo de 2010.

dataset = ee$ImageCollection("CIESIN/GPWv411/GPW_Basic_Demographic_Characteristics")$first();


raster = dataset$select('basic_demographic_characteristics');
raster_vis <- list(
  max = 1000.0,
  palette = c(
    'ffffe7',
    '86a192',
    '509791',
    '307296',
    '2c4484',
    '000066'
  ),
  min= 0.0
);
Map$setCenter(lon = -70.64724, lat = -33.47269, zoom = 8);
Map$addLayer(raster, raster_vis, 'basic_demographic_characteristics');

Sentinel-5P NRTI CO: Near Real-Time Carbon Monoxide

https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S5P_NRTI_L3_CO

ee.ImageCollection(“COPERNICUS/S5P/NRTI/L3_CO”)

library(rgee)
ee_Initialize()

## -- rgee 1.0.9 --------------------------------------- earthengine-api 0.1.259 -- 
##  v email: not_defined
##  v Initializing Google Earth Engine:
 v Initializing Google Earth Engine:  DONE!
## 
 v Earth Engine user: users/tarredwall 
## --------------------------------------------------------------------------------

collection = ee$ImageCollection('COPERNICUS/S5P/NRTI/L3_CO')$select('H2O_column_number_density')$filterDate('2019-06-01', '2019-06-11');

band_viz = list(
  min= 0,
  max= 0.005,
  palette= c('black', 'blue', 'purple', 'cyan', 'green', 'yellow', 'red')
);

Map$addLayer(collection$mean(), band_viz, 'S5P CO');

Map$setCenter(lon = -70.64724, lat = -33.47269, zoom = 8);

Zonal statistics

library(rgee)
# ee_reattach() # reattach ee as a reserved word

ee_Initialize()

## -- rgee 1.0.9 --------------------------------------- earthengine-api 0.1.259 -- 
##  v email: not_defined
##  v Initializing Google Earth Engine:
 v Initializing Google Earth Engine:  DONE!
## 
 v Earth Engine user: users/tarredwall 
## --------------------------------------------------------------------------------

# Load a region representing the United States
region <- ee$FeatureCollection("USDOS/LSIB_SIMPLE/2017")$
  filter(ee$Filter$eq("country_na", "United States"))

# Load MODIS land cover categories in 2001.
landcover <- ee$Image("MODIS/051/MCD12Q1/2001_01_01")$
  select("Land_Cover_Type_1")

# Load nightlights image inputs.
nl2001 <- ee$Image("NOAA/DMSP-OLS/NIGHTTIME_LIGHTS/F152001")$
  select("stable_lights")
nl2012 <- ee$Image("NOAA/DMSP-OLS/NIGHTTIME_LIGHTS/F182012")$
  select("stable_lights")

# Compute the nightlights decadal difference, add land cover codes.
nlDiff <- nl2012$subtract(nl2001)$addBands(landcover)

# Grouped a mean 'reducer': change of nightlights by land cover category.
means <- nlDiff$reduceRegion(
  reducer = ee$Reducer$mean()$group(
    groupField = 1,
    groupName = "code"
  ),
  geometry = region$geometry(),
  scale = 1000,
  maxPixels = 1e8
)

# Print the resultant Dictionary.
print(means$getInfo())

## $groups
## $groups[[1]]
## $groups[[1]]$code
## [1] 0
## 
## $groups[[1]]$mean
## [1] -0.1964631
## 
## 
## $groups[[2]]
## $groups[[2]]$code
## [1] 1
## 
## $groups[[2]]$mean
## [1] -0.1524617
## 
## 
## $groups[[3]]
## $groups[[3]]$code
## [1] 2
## 
## $groups[[3]]$mean
## [1] 0.5306221
## 
## 
## $groups[[4]]
## $groups[[4]]$code
## [1] 3
## 
## $groups[[4]]$mean
## [1] -0.5528383
## 
## 
## $groups[[5]]
## $groups[[5]]$code
## [1] 4
## 
## $groups[[5]]$mean
## [1] -0.2330683
## 
## 
## $groups[[6]]
## $groups[[6]]$code
## [1] 5
## 
## $groups[[6]]$mean
## [1] -0.002707834
## 
## 
## $groups[[7]]
## $groups[[7]]$code
## [1] 6
## 
## $groups[[7]]$mean
## [1] 0.5096342
## 
## 
## $groups[[8]]
## $groups[[8]]$code
## [1] 7
## 
## $groups[[8]]$mean
## [1] 0.3023636
## 
## 
## $groups[[9]]
## $groups[[9]]$code
## [1] 8
## 
## $groups[[9]]$mean
## [1] 0.7629956
## 
## 
## $groups[[10]]
## $groups[[10]]$code
## [1] 9
## 
## $groups[[10]]$mean
## [1] 0.3289744
## 
## 
## $groups[[11]]
## $groups[[11]]$code
## [1] 10
## 
## $groups[[11]]$mean
## [1] 0.1542069
## 
## 
## $groups[[12]]
## $groups[[12]]$code
## [1] 11
## 
## $groups[[12]]$mean
## [1] 0.2643969
## 
## 
## $groups[[13]]
## $groups[[13]]$code
## [1] 12
## 
## $groups[[13]]$mean
## [1] -0.217288
## 
## 
## $groups[[14]]
## $groups[[14]]$code
## [1] 13
## 
## $groups[[14]]$mean
## [1] 0.2775946
## 
## 
## $groups[[15]]
## $groups[[15]]$code
## [1] 14
## 
## $groups[[15]]$mean
## [1] -0.06651379
## 
## 
## $groups[[16]]
## $groups[[16]]$code
## [1] 15
## 
## $groups[[16]]$mean
## [1] 0.08212125
## 
## 
## $groups[[17]]
## $groups[[17]]$code
## [1] 16
## 
## $groups[[17]]$mean
## [1] 0.1719363

1.4 La matemática de las bandas

Haremos una sustracción entre los Índices de Vegetación de Diferencia Normalizada (NDVI) de dos imágenes.

getNDVI <- function(image) {
  return(image$normalizedDifference(c("B4", "B3")))
}
image1 <- ee$Image("LANDSAT/LT05/C01/T1_TOA/LT05_044034_19900604")
image2 <- ee$Image("LANDSAT/LT05/C01/T1_TOA/LT05_044034_20100611")
ndvi1 <- getNDVI(image1)
ndvi2 <- getNDVI(image2)
ndviDifference <- ndvi2$subtract(ndvi1)
ndviParams <- list(palette = c(
  "#d73027", "#f46d43", "#fdae61", "#fee08b",
  "#d9ef8b", "#a6d96a", "#66bd63", "#1a9850"
))
ndwiParams <- list(min = -0.5, max = 0.5, palette = c("FF0000", "FFFFFF", "0000FF"))
Map$centerObject(ndvi1)
Map$addLayer(ndvi1, ndviParams, "NDVI 1") +
Map$addLayer(ndvi2, ndviParams, "NDVI 2") +
Map$addLayer(ndviDifference, ndwiParams, "NDVI difference")

1.5 La función de mapeo

Podemos recorrer todos los datos contenidos en una ImageCollection sin necesidad de usar un ciclo for:

library(rgee)
library(kableExtra)
ee_Initialize()

## -- rgee 1.0.9 --------------------------------------- earthengine-api 0.1.259 -- 
##  v email: not_defined
##  v Initializing Google Earth Engine:
 v Initializing Google Earth Engine:  DONE!
## 
 v Earth Engine user: users/tarredwall 
## --------------------------------------------------------------------------------

addNDVI <- function(image) {
  return(image$addBands(image$normalizedDifference(c("B5", "B4"))))
}

collection <- ee$ImageCollection("LANDSAT/LC08/C01/T1")$
  filterBounds(ee$Geometry$Point(-122.262, 37.8719))$
  filterDate("2014-06-01", "2014-10-01")

ndviCollection <- collection$map(addNDVI)

first <- ndviCollection$first()
a_004 <- first$getInfo()
f_004 <- tail(a_004)
g_004 <- head(f_004,2)

write.csv2(g_004,"g_004.csv")
s <- read.csv2("g_004.csv")
s[,c(1:20)]%>% 
 kbl(caption = "Datos contenidos en LANDSAT/LC08/C01/T1") %>%
  kable_classic(full_width = F, html_font = "Cambria")%>%
  scroll_box(width = "100%", height = "200px")

Datos contenidos en LANDSAT/LC08/C01/T1
X	type	bands.id	bands.data_type.type	bands.data_type.precision	bands.data_type.max	bands.dimensions	bands.crs	bands.crs_transform	bands.id.1	bands.data_type.type.1	bands.data_type.precision.1	bands.data_type.max.1	bands.dimensions.1	bands.crs.1	bands.crs_transform.1	bands.id.2	bands.data_type.type.2
1	Image	B1	PixelType	int	65535	7671	EPSG:32610	30	B2	PixelType	int	65535	7671	EPSG:32610	30	B3	PixelType
2	Image	B1	PixelType	int	65535	7801	EPSG:32610	0	B2	PixelType	int	65535	7801	EPSG:32610	0	B3	PixelType
3	Image	B1	PixelType	int	65535	7671	EPSG:32610	463785	B2	PixelType	int	65535	7671	EPSG:32610	463785	B3	PixelType
4	Image	B1	PixelType	int	65535	7801	EPSG:32610	0	B2	PixelType	int	65535	7801	EPSG:32610	0	B3	PixelType
5	Image	B1	PixelType	int	65535	7671	EPSG:32610	-30	B2	PixelType	int	65535	7671	EPSG:32610	-30	B3	PixelType
6	Image	B1	PixelType	int	65535	7801	EPSG:32610	4264515	B2	PixelType	int	65535	7801	EPSG:32610	4264515	B3	PixelType

1.6 La función Reduce.

Podemos calcular la mediana de cada banda de las última 5 escenas de nubes:

collection <- ee$ImageCollection("LANDSAT/LC08/C01/T1")$
  filterBounds(ee$Geometry$Point(-122.262, 37.8719))$
  filterDate("2014-01-01", "2014-12-31")$
  sort("CLOUD_COVER")

median <- collection$limit(5)$reduce(ee$Reducer$median())

vizParams <- list(
  bands = c("B5_median", "B4_median", "B3_median"),
  min = 5000, max = 15000, gamma = 1.3
)

Map$addLayer(
  eeObject = median,
  visParams = vizParams,
  name = "Median image"
)

1.7 Estadísticas de las imágenes

Cargamos y desplegamos una imagen de la parte superior de la atmósfera calibrada (TOA) del Landsat 8:

image <- ee$Image("LANDSAT/LC08/C01/T1_TOA/LC08_044034_20140318")
Map$addLayer(
  eeObject = image,
  visParams = list(bands = c("B4", "B3", "B2"), max = 30000),
  name = "Landsat 8"
)

Creamos un rectángulo arbitrario como una región y lo desplegamos:

region <- ee$Geometry$Rectangle(-122.2806, 37.1209, -122.0554, 37.2413)
Map$addLayer(
  eeObject = region,
  name = "Region"
)

Obtenemos un diccionario de las medias en la región:

mean <- image$reduceRegion(
  reducer = ee$Reducer$mean(),
  geometry = region,
  scale = 30
)

print(mean$getInfo())

## $B1
## [1] 0.1015666
## 
## $B10
## [1] 288.3895
## 
## $B11
## [1] 287.4494
## 
## $B2
## [1] 0.07771911
## 
## $B3
## [1] 0.0553003
## 
## $B4
## [1] 0.03670828
## 
## $B5
## [1] 0.2213975
## 
## $B6
## [1] 0.07888247
## 
## $B7
## [1] 0.03508945
## 
## $B8
## [1] 0.04710009
## 
## $B9
## [1] 0.00156097
## 
## $BQA
## [1] 2720.119

1.8 Enmascaramiento

Aunque la composición con la mediana es una mejora, es posible enmascarar partes de la imagen. Enmascarar píxeles de una imagen hace que esos píxeles sean transparentes y los excluye del análisis. Cada píxel de cada banda de una imagen tiene una máscara. Aquellos con un valor de máscara de 0 o menos serán transparentes. Aquellos con una máscara de cualquier valor por encima de 0 serán renderizados.

En el siguiente código creamos una función que obtiene una diferencia normalizada entre dos bandas, cargamos dos imágenes del Landsat 5, aplicamos la función, calculamos su diferencia, cargamos la máscara como archivos de demostración de videojuego (DEM) desde la Misión de topografía de radar Shuttle (SRTM), actualizamos la máscara de diferencia NDVI con la landMask y visualizamos el resultado enmascarado.

getNDVI <- function(image) {
  return(image$normalizedDifference(c("B4", "B3")))
}

image1 <- ee$Image("LANDSAT/LT05/C01/T1_TOA/LT05_044034_19900604")
image2 <- ee$Image("LANDSAT/LT05/C01/T1_TOA/LT05_044034_20100611")

ndvi1 <- getNDVI(image1)
ndvi2 <- getNDVI(image2)

ndviDifference <- ndvi2$subtract(ndvi1)
landMask <- ee$Image("CGIAR/SRTM90_V4")$mask()

maskedDifference <- ndviDifference$updateMask(landMask)

vizParams <- list(
  min = -0.5,
  max = 0.5,
  palette = c("FF0000", "FFFFFF", "0000FF")
)

Map$addLayer(
  eeObject = maskedDifference,
  visParams = vizParams,
  name = "NDVI difference"
)

2. Aprendizaje automático

(volver al índice)

2.1 Clasificación supervisada

2.1.1 Árboles de clasificación y regresión (CART)

Tuvimos que hacer un cambio en el código contenido en nuestra página madre. cart da el siguiente error:

cart: Error in py_call_impl(callable, dots$args, dots$keywords) : EEException: Classifier.cart: This classifier has been removed. For more information see: http://goo.gle/deprecated-classifiers.

Debió ser reemplazada por: smileCart

Referencia: https://developers.google.com/earth-engine/transitions#new-style-classifiers

Las imágenes de entrada son un compuesto Landsat 8 sin nubes, cuyas bandas utilizamos para la predicción. Cargamos los puntos de entrenamiento. La propiedad numérica ‘clase’ almacena etiquetas conocidas. Ésta propiedad de la tabla almacena las etiquetas de cobertura terrestre. Superponemos los puntos en las imágenes para entrenar. Entrenamos un clasificador CART con parámetros predeterminados. Clasificamos la imagen con las mismas bandas utilizadas para el entrenamiento. Establecemos los parámetros de visualización y visualizamos las entradas y los resultados.

# Input imagery is a cloud-free Landsat 8 composite.
l8 <- ee$ImageCollection("LANDSAT/LC08/C01/T1")

image <- ee$Algorithms$Landsat$simpleComposite(
  collection = l8$filterDate("2018-01-01", "2018-12-31"),
  asFloat = TRUE
)

# Use these bands for prediction.
bands <- c("B2", "B3", "B4", "B5", "B6", "B7", "B10", "B11")

# Load training points. The numeric property 'class' stores known labels.
points <- ee$FeatureCollection("GOOGLE/EE/DEMOS/demo_landcover_labels")

# This property of the table stores the land cover labels.
label <- "landcover"

# Overlay the points on the imagery to get training.
training <- image$select(bands)$sampleRegions(
  collection = points,
  properties = list(label),
  scale = 30
)

# Train a CART classifier with default parameters.
trained <- ee$Classifier$smileCart()$train(training, label, bands)

# Classify the image with the same bands used for training.
classified <- image$select(bands)$classify(trained)

# Viz parameters.
viz_img <- list(bands = c("B4", "B3", "B2"), max = 0.4)
viz_class <- list(palette = c("red", "green", "blue"), min = 0, max = 2)

# Display the inputs and the results.
Map$centerObject(points)

## NOTE: Center obtained from the first element.

Map$addLayer(image, viz_img, name = "image") +
Map$addLayer(classified, viz_class, name = "classification")

2.1.2 Máquinas de soporte vectorial (SVM)

svm: Error in py_call_impl(callable, dots$args, dots$keywords) : EEException: Classifier.svm: This classifier has been replaced. For more information see: http://goo.gle/deprecated-classifiers.

tuvo que ser reemplazado por libsvm

Referencia: https://developers.google.com/earth-engine/guides/classification

# Input imagery is a cloud-free Landsat 8 composite.
l8 = ee$ImageCollection("LANDSAT/LC08/C01/T1")

image = ee$Algorithms$Landsat$simpleComposite(
  collection = l8$filterDate("2018-01-01", "2018-12-31"),
  asFloat = TRUE
)

# Use these bands for prediction.
bands = c("B2", "B3", "B4", "B5", "B6", "B7", "B10", "B11")

# Manually created polygons.
forest1 = ee$Geometry$Rectangle(-63.0187, -9.3958, -62.9793, -9.3443)
forest2 = ee$Geometry$Rectangle(-62.8145, -9.206, -62.7688, -9.1735)
nonForest1 = ee$Geometry$Rectangle(-62.8161, -9.5001, -62.7921, -9.4486)
nonForest2 = ee$Geometry$Rectangle(-62.6788, -9.044, -62.6459, -8.9986)

# Make a FeatureCollection from the hand-made geometries.
polygons = ee$FeatureCollection(c(
  ee$Feature(nonForest1, list(class = 0)),
  ee$Feature(nonForest2, list(class = 0)),
  ee$Feature(forest1, list(class = 1)),
  ee$Feature(forest2, list(class = 1))
))

# Get the values for all pixels in each polygon in the training.
training = image$sampleRegions(
  # Get the sample from the polygons FeatureCollection.
  collection = polygons,
  # Keep this list of properties from the polygons.
  properties = list("class"),
  # Set the scale to get Landsat pixels in the polygons.
  scale = 30
)

# Create an SVM classifier with custom parameters.
classifier = ee$Classifier$libsvm(
  kernelType = "RBF",
  gamma = 0.5,
  cost = 10
)

# Train the classifier.
trained = classifier$train(training, "class", bands)

# Classify the image.
classified = image$classify(trained)

# Display the classification result and the input image.
geoviz_image = list(bands = c("B4", "B3", "B2"), max = 0.5, gamma = 2)
geoviz_class = list(min = 0, max = 1, palette = c("red", "green"))

Map$setCenter(-62.836, -9.2399, 9)
Map$addLayer(
  eeObject = image,
  visParams = geoviz_image,
  name = "image"
) +
Map$addLayer(
  eeObject = classified,
  visParams = geoviz_class,
  name = "deforestation"
) +
Map$addLayer(
  eeObject = polygons,
  name = "training polygons"
)

2.2 Clasificación no supervisada

2.2.1 Clustering por Kmeans

# Load a pre-computed Landsat composite for input.
input <- ee$Image("LANDSAT/LE7_TOA_1YEAR/2001")

# Define a region in which to generate a sample of the input.
region <- ee$Geometry$Rectangle(29.7, 30, 32.5, 31.7)

# Display the sample region.
Map$setCenter(31.5, 31.0, 7)
Map$addLayer(
  eeObject = ee$Image()$paint(region, 0, 2),
  name = "region"
)

# Make the training dataset.
training <- input$sample(
  region = region,
  scale = 30,
  numPixels = 5000
)

# Instantiate the clusterer and train it.
clusterer <- ee$Clusterer$wekaKMeans(15)$train(training)

# Cluster the input using the trained clusterer.
result <- input$cluster(clusterer)

# Display the clusters with random colors.
Map$centerObject(region)
Map$addLayer(
  eeObject = result$randomVisualizer(),
  name = "clusters"
)

3. Imágenes

(volver al índice)

3.1 Introducción

Crea una imagen constante con un valor de píxel de 1 Crea una imagen constante con un valor de píxel de 1 Concatenar imágenes Cambie la opción de impresión por: “simplemente”, “json”, “ee_print” Cambiar el nombre de las imágenes usando ee $ Image $ select

# Create a constant Image with a pixel value of 1
image1 <- ee$Image(1)
print(image1, type = "json")

## ee.Image({
##   "functionInvocationValue": {
##     "functionName": "Image.constant",
##     "arguments": {
##       "value": {
##         "constantValue": 1.0
##       }
##     }
##   }
## })

print(image1, type = "simply")

## EarthEngine Object: Image

print(image1, type = "ee_print")

## ---------------------------------------------------------- Earth Engine Image --
## Image Metadata:
##  - Class                      : ee$Image
##  - ID                         : no_id
##  - Number of Bands            : 1
##  - Bands names                : constant
##  - Number of Properties       : 0
##  - Number of Pixels*          : 64800
##  - Approximate size*          : 101.25 KB
## Band Metadata (img_band = constant):
##  - EPSG (SRID)                : WGS 84 (EPSG:4326)
##  - proj4string                : +proj=longlat +datum=WGS84 +no_defs
##  - Geotransform               : 1 0 0 0 1 0
##  - Nominal scale (meters)     : 111319.5
##  - Dimensions                 : 360 180
##  - Number of Pixels           : 64800
##  - Data type                  : INT
##  - Approximate size           : 101.25 KB
##  --------------------------------------------------------------------------------
##  NOTE: (*) Properties calculated considering a constant geotransform and data type.

# Create a constant Image with a pixel value of 1
image2 <- ee$Image(2)

# Concatenate Images
image3 <- ee$Image$cat(c(image1, image2))
ee_print(image3, clean = TRUE)

## ---------------------------------------------------------- Earth Engine Image --
## Image Metadata:
##  - Class                      : ee$Image
##  - ID                         : no_id
##  - Number of Bands            : 2
##  - Bands names                : constant constant_1
##  - Number of Properties       : 0
##  - Number of Pixels*          : 129600
##  - Approximate size*          : 202.50 KB
## Band Metadata (img_band = constant):
##  - EPSG (SRID)                : WGS 84 (EPSG:4326)
##  - proj4string                : +proj=longlat +datum=WGS84 +no_defs
##  - Geotransform               : 1 0 0 0 1 0
##  - Nominal scale (meters)     : 111319.5
##  - Dimensions                 : 360 180
##  - Number of Pixels           : 64800
##  - Data type                  : INT
##  - Approximate size           : 101.25 KB
##  --------------------------------------------------------------------------------
##  NOTE: (*) Properties calculated considering a constant geotransform and data type.

# Change print option by: "simply","json","ee_print"
options(rgee.print.option = "simply")
multiband <- ee$Image(c(1, 2, 3))
print(multiband)

## EarthEngine Object: Image

# Rename images using ee$Image$select
renamed <- multiband$select(
  opt_selectors = c("constant", "constant_1", "constant_2"),
  opt_names = c("band1", "band2", "band3")
)
ee_print(renamed)

## ---------------------------------------------------------- Earth Engine Image --
## Image Metadata:
##  - Class                      : ee$Image
##  - ID                         : no_id
##  - Number of Bands            : 3
##  - Bands names                : band1 band2 band3
##  - Number of Properties       : 0
##  - Number of Pixels*          : 194400
##  - Approximate size*          : 303.75 KB
## Band Metadata (img_band = band1):
##  - EPSG (SRID)                : WGS 84 (EPSG:4326)
##  - proj4string                : +proj=longlat +datum=WGS84 +no_defs
##  - Geotransform               : 1 0 0 0 1 0
##  - Nominal scale (meters)     : 111319.5
##  - Dimensions                 : 360 180
##  - Number of Pixels           : 64800
##  - Data type                  : INT
##  - Approximate size           : 101.25 KB
##  --------------------------------------------------------------------------------
##  NOTE: (*) Properties calculated considering a constant geotransform and data type.

3.2 Visualización de imágenes

# Simple RGB Vizualization -Landsat
image <- ee$Image("LANDSAT/LC8_L1T_TOA/LC80440342014077LGN00")
vizParams <- list(
  min = 0,
  max = 0.5,
  bands = c("B5", "B4", "B3"),
  gamma = c(0.95, 1.1, 1)
)
Map$setCenter(-122.1899, 37.5010, 5)
Map$addLayer(image, vizParams)

# Simple Single-Band Vizualization
ndwi <- image$normalizedDifference(c("B3", "B5"))
ndwiViz <- list(
  min = 0.5,
  max = 1,
  palette = c("00FFFF", "0000FF")
)
## Masking
ndwiMasked <- ndwi$updateMask(ndwi$gte(0.4))

## Produce an RGB layers.
imageRGB <- image$visualize(
  list(
    bands = c("B5", "B4", "B3"),
    max = 0.5
  )
)

ndwiRGB <- ndwiMasked$visualize(
  list(
    min = 0.5,
    max = 1,
    palette = c("00FFFF", "0000FF")
  )
)

# Mosaic the visualization layers and display( or export).
roi <- ee$Geometry$Point(c(-122.4481, 37.7599))$buffer(20000)
Map$centerObject(image$clip(roi))

Map$addLayer(
  eeObject = image$clip(roi),
  visParams = vizParams,
  name = "Landsat 8"
) +
  Map$addLayer(
    eeObject = ndwiMasked$clip(roi),
    visParams = ndwiViz,
    name = "NDWI"
  )

3.3 Información sobre las imágenes y su metadata

# # Load an Landsat8 Image - San Francisco, California.
# image <- ee$Image("LANDSAT/LC8_L1T/LC80440342014077LGN00")
# 
# # Get Band names of an ee$Image
# bandNames <- image$bandNames()
# #ee_help(bandNames)
# cat("Band names: ", paste(bandNames$getInfo(), collapse = " "))
# 
# # Get Band names of an ee$Image
# b1proj <- image$select("B1")$projection()$getInfo()
# cat("Band 1 projection: ")
# cat("type: ", b1proj$type)
# cat("crs: ", b1proj$crs)
# cat("geotransform: ", paste0(b1proj$transform, " "))
# 
# b1scale <- image$select("B1")$projection()$nominalScale()
# cat("Band 1 scale: ", b1scale$getInfo())
# 
# b8scale <- image$select("B8")$projection()$nominalScale()
# cat("Band 8 scale: ", b8scale$getInfo())
# 
# properties <- image$propertyNames()$getInfo()
# cat("Metadata properties: \n-", paste0(properties, collapse = "\n- "))
# 
# cloudiness <- image$get("CLOUD_COVER")$getInfo()
# cat("CLOUD_COVER: ", cloudiness)
# 
# iso_date <- eedate_to_rdate(image$get("system:time_start"))
# iso_timestamp <- eedate_to_rdate(
#   ee_date = image$get("system:time_start"),
#   js = TRUE
# )
# cat("ISO Date: ", as.character(iso_date))
# cat("Timestamp : ", format(iso_timestamp, scientific = FALSE))

3.4 Operaciones matemáticas

# Load two 5-year Landsat 7 composites.
landsat1999 <- ee$Image("LANDSAT/LE7_TOA_5YEAR/1999_2003")
landsat2008 <- ee$Image("LANDSAT/LE7_TOA_5YEAR/2008_2012")

# Compute NDVI the hard way.
ndvi1999 <- landsat1999$select("B4")$subtract(landsat1999$select("B3"))$divide(landsat1999$select("B4")$add(landsat1999$select("B3")))

# Compute NDVI the easy way.
ndvi2008 <- landsat2008$normalizedDifference(c("B4", "B3"))

# Compute the multi-band difference image.
diff <- landsat2008$subtract(landsat1999)

# Compute the squared difference in each band.
squaredDifference <- diff$pow(2)

Map$setZoom(zoom = 3)
Map$addLayer(
  eeObject = diff,
  visParams = list(bands = c("B4", "B3", "B2"), min = -32, max = 32),
  name = "difference"
) +
  Map$addLayer(
    eeObject = squaredDifference,
    visParams = list(bands = c("B4", "B3", "B2"), max = 1000),
    name = "squared diff."
  )

3.5 Relational, conditional and Boolean operations

# Load a Landsat 8 image.
image <- ee$Image("LANDSAT/LC08/C01/T1_TOA/LC08_044034_20140318")

# Create NDVI and NDWI spectral indices.
ndvi <- image$normalizedDifference(c("B5", "B4"))
ndwi <- image$normalizedDifference(c("B3", "B5"))

# Create a binary layer using logical operations.
bare <- ndvi$lt(0.2)$And(ndwi$lt(0))

# Mask and display the binary layer.
Map$setCenter(lon = -122.3578, lat = 37.7726)
Map$setZoom(zoom = 10)

Map$addLayer(
  eeObject = bare$updateMask(bare),
  visParams = list(min = 0, max = 20),
  name = "bare"
)

3.6 Convoluciones

# Load and display an image.
image <- ee$Image("LANDSAT/LC08/C01/T1_TOA/LC08_044034_20140318")

# Define a boxcar or low-pass kernel.
# boxcar = ee$Kernel$square(list(
#   radius = 7, units = 'pixels', 'normalize' = T
# ))

boxcar <- ee$Kernel$square(7, "pixels", T)

# Smooth the image by convolving with the boxcar kernel.
smooth <- image$convolve(boxcar)

# Define a Laplacian, or edge-detection kernel.
laplacian <- ee$Kernel$laplacian8(1, F)

# Apply the edge-detection kernel.
edgy <- image$convolve(laplacian)

Map$setCenter(lon = -121.9785, lat = 37.8694)
Map$setZoom(zoom = 11)

Map$addLayer(
  eeObject = image,
  visParams = list(bands = c("B5", "B4", "B3"), max = 0.5),
  name = "input image"
) +
  Map$addLayer(
    eeObject = smooth,
    visParams = list(bands = c("B5", "B4", "B3"), max = 0.5),
    name = "smoothed"
  ) +
  Map$addLayer(
    eeObject = edgy,
    visParams = list(bands = c("B5", "B4", "B3"), max = 0.5),
    name = "edges"
  )

3.7 Operaciones morfológicas

# Load a Landsat 8 image, select the NIR band, threshold, display.
image <- ee$Image("LANDSAT/LC08/C01/T1_TOA/LC08_044034_20140318")$select(4)$gt(0.2)

# Define a kernel.
kernel <- ee$Kernel$circle(radius = 1)

# Perform an erosion followed by a dilation, display.
opened <- image$focal_min(kernel = kernel, iterations = 2)$focal_max(kernel = kernel, iterations = 2)

Map$setCenter(lon = -122.1899, lat = 37.5010)
Map$setZoom(zoom = 13)

Map$addLayer(
  eeObject = image,
  visParams = list(min = 0, max = 1),
  name = "NIR threshold"
) +
  Map$addLayer(
    eeObject = opened,
    visParams = list(min = 0, max = 1),
    name = "opened"
  )

3.8 Gradientes

# Load a Landsat 8 image and select the panchromatic band.
image <- ee$Image("LANDSAT/LC08/C01/T1/LC08_044034_20140318")$select("B8")

# Compute the image gradient in the X and Y directions.
xyGrad <- image$gradient()

# Compute the magnitude of the gradient.
gradient <- xyGrad$select("x")$pow(2)$add(xyGrad$select("y")$pow(2))$sqrt()

# Compute the direction of the gradient.
direction <- xyGrad$select("y")$atan2(xyGrad$select("x"))

# Display the results.
Map$setCenter(lon = -122.054, lat = 37.7295)
Map$setZoom(zoom = 10)

Map$addLayer(
  eeObject = direction,
  visParams = list(min = -2, max = 2),
  name = "direction"
) +
  Map$addLayer(
    eeObject = gradient,
    visParams = list(min = -7, max = 7),
    name = "opened"
  )

3.9 Detección de edges

# Load a Landsat 8 image, select the panchromatic band$
image <- ee$Image("LANDSAT/LC08/C01/T1/LC08_044034_20140318")$select("B8")

# Perform Canny edge detection and display the result$
canny <- ee$Algorithms$CannyEdgeDetector(image, threshold = 10, sigma = 1)

# Perform Hough transform of the Canny result and display$
hough <- ee$Algorithms$HoughTransform(canny, 256, 600, 100)

# Load a Landsat 8 image, select the panchromatic band$
image <- ee$Image("LANDSAT/LC08/C01/T1/LC08_044034_20140318")$select("B8")

# Define a "fat" Gaussian kernel$
fat <- ee$Kernel$gaussian(
  radius = 3,
  sigma = 3,
  units = "pixels",
  normalize = T,
  magnitude = -1
)

# Define a "skinny" Gaussian kernel.
skinny <- ee$Kernel$gaussian(
  radius = 3,
  sigma = 1,
  units = "pixels",
  normalize = T
)

# Compute a difference-of-Gaussians (DOG) kernel.
dog <- fat$add(skinny)

# Compute the zero crossings of the second derivative, display.
zeroXings <- image$convolve(dog)$zeroCrossing()

Map$setCenter(lon = -122.054, lat = 37.7295)
Map$setZoom(zoom = 10)

Map$addLayer(
  eeObject = canny,
  visParams = list(min = 0, max = 1),
  name = "canny"
) +
  Map$addLayer(
    eeObject = hough,
    visParams = list(min = 0, max = 1),
    name = "hough"
  ) +
  Map$addLayer(
    eeObject = image,
    visParams = list(min = 0, max = 12000),
    name = "L_B8"
  ) +
  Map$addLayer(
    eeObject = zeroXings$updateMask(zeroXings),
    visParams = list(palette = "FF0000"),
    name = "zero crossings"
  )

continuará…

4. ImageCollection

(volver al índice)

4.2 ImageCollection Information and Metadata

# # Load a Landsat 8 ImageCollection for a single path-row.
# 
# collection <- ee$ImageCollection("LANDSAT/LC08/C01/T1_TOA")$
#   filter(ee$Filter$eq("WRS_PATH", 44))$
#   filter(ee$Filter$eq("WRS_ROW", 34))$
#   filterDate("2014-03-01", "2014-09-01")
# ee_print(collection)
# 
# # Get the number of images.
# 
# count <- collection$size()
# cat("Count: ", count$getInfo())
# 
# # Get the date range of images in the collection.
# 
# range <- collection$reduceColumns(
#   ee$Reducer$minMax(),
#   list("system:time_start")
# )
# 
# col_min <- eedate_to_rdate(range$get("min"))
# col_max <- eedate_to_rdate(range$get("max"))
# cat("Date range: ", as.character(col_min), as.character(col_max))
# 
# # Get statistics for a property of the images in the collection.
# 
# sunStats <- collection$aggregate_stats("SUN_ELEVATION")
# cat("Sun elevation statistics: ")
# sunStats$getInfo()
# 
# # Sort by a cloud cover property, get the least cloudy image.
# 
# image <- ee$Image(collection$sort("CLOUD_COVER")$first())
# cat("Least cloudy image: ")
# image$getInfo()
# 
# # Limit the collection to the 10 most recent images.
# 
# recent <- collection$sort("system:time_start", FALSE)$limit(10)
# cat("Recent images: ")
# recent$getInfo()

4.3 Filtering an ImageCollection

library(rgee)
# ee_reattach() # reattach ee as a reserved word

ee_Initialize()

## -- rgee 1.0.9 --------------------------------------- earthengine-api 0.1.259 -- 
##  v email: not_defined
##  v Initializing Google Earth Engine:
 v Initializing Google Earth Engine:  DONE!
## 
 v Earth Engine user: users/tarredwall 
## --------------------------------------------------------------------------------

# Load Landsat 5 data, filter by date and bounds.
collection <- ee$ImageCollection("LANDSAT/LT05/C01/T2")$
  filterDate("1987-01-01", "1990-05-01")$
  filterBounds(ee$Geometry$Point(25.8544, -18.08874))

# Also filter the collection by the IMAGE_QUALITY property.
filtered <- collection$filterMetadata("IMAGE_QUALITY", "equals", 9)

# Create two composites to check the effect of filtering by IMAGE_QUALITY.
badComposite <- ee$Algorithms$Landsat$simpleComposite(collection, 75, 3)
goodComposite <- ee$Algorithms$Landsat$simpleComposite(filtered, 75, 3)

# Display the composites.
Map$setCenter(lon = 25.8544, lat = -18.08874)
Map$setZoom(zoom = 13)

Map$addLayer(
  eeObject = badComposite,
  visParams = list(bands = c("B3", "B2", "B1"), gain = 3.5),
  name = "bad composite"
) +
Map$addLayer(
  eeObject = goodComposite,
  visParams = list(bands = c("B3", "B2", "B1"), gain = 3.5),
  name = "good composite"
)

4.4 Mapping over an ImageCollection

# 
# library(rgee)
# ee_Initialize()
# 
# # This function adds a band representing the image timestamp.
# addTime <- function(image) {
#   return(image$addBands(image$metadata("system:time_start")))
# }
# 
# conditional <- function(image) {
#   return(ee$Algorithms$If(
#     ee$Number(image$get("SUN_ELEVATION"))$gt(40),
#     image,
#     ee$Image(0)
#   ))
# }
# 
# # Load a Landsat 8 collection for a single path-row.
# collection <- ee$ImageCollection("LANDSAT/LC08/C01/T1_TOA")$
#   filter(ee$Filter$eq("WRS_PATH", 44))$
#   filter(ee$Filter$eq("WRS_ROW", 34))
# 
# # Map the function over the collection and display the result.
# print(collection$map(addTime)$getInfo())
# 
# 
# # Load a Landsat 8 collection for a single path-row.
# collection <- ee$ImageCollection("LANDSAT/LC8_L1T_TOA")$
#   filter(ee$Filter$eq("WRS_PATH", 44))$
#   filter(ee$Filter$eq("WRS_ROW", 34))
# 
# # This function uses a conditional statement to return the image if
# # the solar elevation > 40 degrees.  Otherwise it returns a zero image.
# # conditional = function(image) {
# #   return ee.Algorithms.If(ee.Number(image.get('SUN_ELEVATION')).gt(40),
# #                           image,
# #                           ee.Image(0))
# # }
# 
# # Map the function over the collection, convert to a List and print the result.
# print(collection$map(conditional)$getInfo())

4.5 Reducing an ImageCollection

library(rgee)
# ee_reattach() # reattach ee as a reserved word

ee_Initialize()

## -- rgee 1.0.9 --------------------------------------- earthengine-api 0.1.259 -- 
##  v email: not_defined
##  v Initializing Google Earth Engine:
 v Initializing Google Earth Engine:  DONE!
## 
 v Earth Engine user: users/tarredwall 
## --------------------------------------------------------------------------------

addTime <- function(image) {
  image$addBands(image$metadata("system:time_start")$divide(1000 * 60 * 60 * 24 * 365))
}

# Load a Landsat 8 collection for a single path-row.
collection <- ee$ImageCollection("LANDSAT/LC08/C01/T1_TOA")$
  filter(ee$Filter$eq("WRS_PATH", 44))$
  filter(ee$Filter$eq("WRS_ROW", 34))$
  filterDate("2014-01-01", "2015-01-01")

# Compute a median image and display.
median <- collection$median()
Map$setCenter(lon = -122.3578, lat = 37.7726)
Map$setZoom(zoom = 12)

Map$addLayer(
  eeObject = median,
  visParams = list(bands = c("B4", "B3", "B2"), max = 0.3),
  name = "median"
)

# Reduce the collection with a median reducer.
median <- collection$reduce(ee$Reducer$median())

# Display the median image.
Map$addLayer(
  eeObject = median,
  visParams = list(bands = c("B4_median", "B3_median", "B2_median"), max = 0.3),
  name = "also median"
)

# # This function adds a band representing the image timestamp.
# addTime = function(image) {
#   return image.addBands(image.metadata('system:time_start')
#     # Convert milliseconds from epoch to years to aid in
#     # interpretation of the following trend calculation. \
#     .divide(1000 * 60 * 60 * 24 * 365))
# }

# Load a MODIS collection, filter to several years of 16 day mosaics,
# and map the time band function over it.
collection <- ee$ImageCollection("MODIS/006/MYD13A1")$
  filterDate("2004-01-01", "2010-10-31")$
  map(addTime)

# Select the bands to model with the independent variable first.
trend <- collection$select(list("system:time_start", "EVI"))$
  reduce(ee$Reducer$linearFit())

# Display the trend with increasing slopes in green, decreasing in red.
Map$setCenter(lon = -96.943, lat = 39.436)
Map$setZoom(zoom = 5)

Map$addLayer(
  eeObject = trend,
  visParams = list(bands = c("scale", "scale", "offset"), min = 0, max = c(-100, 100, 10000)),
  name = "EVI trend"
)

4.6 Compositing and Mosaicking

# Load three NAIP quarter quads in the same location, different times.
naip2004_2012 <- ee$ImageCollection("USDA/NAIP/DOQQ")$
  filterBounds(ee$Geometry$Point(-71.08841, 42.39823))$
  filterDate("2004-07-01", "2012-12-31")$
  select(c("R", "G", "B"))

# Temporally composite the images with a maximum value function.
composite <- naip2004_2012$max()
Map$setCenter(lon = -71.12532, lat = 42.3712)
Map$setZoom(zoom = 12)

Map$addLayer(
  eeObject = composite,
  visParams = list(),
  name = "max value composite"
)

# Load four 2012 NAIP quarter quads, different locations.
naip2012 <- ee$ImageCollection("USDA/NAIP/DOQQ")$
  filterBounds(ee$Geometry$Rectangle(-71.17965, 42.35125, -71.08824, 42.40584))$
  filterDate("2012-01-01", "2012-12-31")

# Spatially mosaic the images in the collection and display.
mosaic <- naip2012$mosaic()
Map$addLayer(
  eeObject = mosaic,
  visParams = list(),
  name = "spatial mosaic"
)

# Load a NAIP quarter quad, display.
naip <- ee$Image("USDA/NAIP/DOQQ/m_4207148_nw_19_1_20120710")
Map$setCenter(lon = -71.0915, lat = 42.3443)
Map$setZoom(zoom = 14)

Map$addLayer(
  eeObject = naip,
  visParams = list(),
  name = "NAIP DOQQ"
)

# Create the NDVI and NDWI spectral indices.
ndvi <- naip$normalizedDifference(c("N", "R"))
ndwi <- naip$normalizedDifference(c("G", "N"))

# Create some binary images from thresholds on the indices.
# This threshold is designed to detect bare land.
bare1 <- ndvi$lt(0.2)$And(ndwi$lt(0.3))
# This detects bare land with lower sensitivity. It also detects shadows.
bare2 <- ndvi$lt(0.2)$And(ndwi$lt(0.8))

# Define visualization parameters for the spectral indices.
ndviViz <- list(min = -1, max = 1, palette = c("FF0000", "00FF00"))
ndwiViz <- list(min = 0.5, max = 1, palette = c("00FFFF", "0000FF"))

# Mask and mosaic visualization images.  The last layer is on top.
mosaic <- ee$ImageCollection(list(
  # NDWI > 0.5 is water.  Visualize it with a blue palette.
  ndwi$updateMask(ndwi$gte(0.5))$visualize(ndwiViz),
  # NDVI > 0.2 is vegetation.  Visualize it with a green palette.
  ndvi$updateMask(ndvi$gte(0.2))$visualize(ndviViz),
  # Visualize bare areas with shadow (bare2 but not bare1) as gray.
  bare2$updateMask(bare2$And(bare1$Not()))$visualize(list(palette = c("AAAAAA"))),
  # Visualize the other bare areas as white.
  bare1$updateMask(bare1)$visualize(list(palette = c("FFFFFF")))
))$mosaic()

Map$addLayer(
  eeObject = ndwi,
  visParams = list(),
  name = "Visualization mosaic"
)

Continuará…

5. Geometry, Feature y FeatureCollection

(volver al índice)

5.2 Geodesic vs. Planar Geometries

# Create a geodesic polygon.
polygon <- ee$Geometry$Polygon(
  list(
    c(-5, 40), c(65, 40), c(65, 60), c(-5, 60), c(-5, 60)
  )
)

# Create a planar polygon.
planarPolygon <- ee$Geometry(polygon, {}, FALSE)
polygon <- ee$FeatureCollection(polygon)
planarPolygon <- ee$FeatureCollection(planarPolygon)

# Display the polygons by adding them to the map.
Map$centerObject(polygon, zoom = 2)

## NOTE: Center obtained from the first element.

Map$addLayer(polygon, list(color = "FF0000"), "geodesic polygon") +
Map$addLayer(planarPolygon, list(color = "000000"), "planar polygon")

5.3 Geometry Visualization and Information

# Create a geodesic polygon.
polygon <- ee$Geometry$Polygon(
  list(
    c(-5, 40), c(65, 40), c(65, 60), c(-5, 60), c(-5, 60)
  )
)

# Create a planar polygon.
planarPolygon <- ee$Geometry(polygon, {}, FALSE)
polygon <- ee$FeatureCollection(polygon)
planarPolygon <- ee$FeatureCollection(planarPolygon)

# Display the polygons by adding them to the map.
Map$centerObject(polygon, zoom = 2)

## NOTE: Center obtained from the first element.

Map$addLayer(polygon, list(color = "FF0000"), "geodesic polygon") +
Map$addLayer(planarPolygon, list(color = "000000"), "planar polygon")

# Create two circular geometries.
poly1 <- ee$Geometry$Point(c(-50, 30))$buffer(1e6)
poly2 <- ee$Geometry$Point(c(-40, 30))$buffer(1e6)

# Display polygon 1 in red and polygon 2 in blue.
Map$setCenter(-45, 30)
Map$addLayer(poly1, list(color = "FF0000"), "poly1") +
Map$addLayer(poly2, list(color = "0000FF"), "poly2")

# Compute the intersection, display it in blue.
intersection <- poly1$intersection(poly2, ee$ErrorMargin(1))
Map$addLayer(intersection, list(color = "00FF00"), "intersection")

# Compute the union, display it in magenta.
union <- poly1$union(poly2, ee$ErrorMargin(1))
Map$addLayer(union, list(color = "FF00FF"), "union")

# Compute the difference, display in yellow.
diff1 <- poly1$difference(poly2, ee$ErrorMargin(1))
Map$addLayer(diff1, list(color = "FFFF00"), "diff1")

# Compute symmetric difference, display in black.
symDiff <- poly1$symmetricDifference(poly2, ee$ErrorMargin(1))
Map$addLayer(symDiff, list(color = "000000"), "symmetric difference")

# Create a geodesic polygon.
polygon <- ee$Geometry$Polygon(
  list(
    c(-5, 40), c(65, 40), c(65, 60), c(-5, 60), c(-5, 60)
  )
)

# Create a planar polygon.
planarPolygon <- ee$Geometry(polygon, {}, FALSE)
polygon <- ee$FeatureCollection(polygon)
planarPolygon <- ee$FeatureCollection(planarPolygon)

# Display the polygons by adding them to the map.
Map$centerObject(polygon, zoom = 2)

## NOTE: Center obtained from the first element.

Map$addLayer(polygon, list(color = "FF0000"), "geodesic polygon") +
Map$addLayer(planarPolygon, list(color = "000000"), "planar polygon")

library(rgee)
# ee_reattach() # reattach ee as a reserved word

ee_Initialize()

## -- rgee 1.0.9 --------------------------------------- earthengine-api 0.1.259 -- 
##  v email: not_defined
##  v Initializing Google Earth Engine:
 v Initializing Google Earth Engine:  DONE!
## 
 v Earth Engine user: users/tarredwall 
## --------------------------------------------------------------------------------

# This function gets NDVI from Landsat 8 imagery.
addNDVI <- function(image) {
  return(image$addBands(image$normalizedDifference(c("B5", "B4"))))
}

# This function masks cloudy pixels.
cloudMask <- function(image) {
  clouds <- ee$Algorithms$Landsat$simpleCloudScore(image)$select("cloud")
  return(image$updateMask(clouds$lt(10)))
}

# Load a Landsat collection, map the NDVI and cloud masking functions over it.
collection <- ee$ImageCollection("LANDSAT/LC08/C01/T1_TOA")$
  filterBounds(ee$Geometry$Point(c(-122.262, 37.8719)))$
  filterDate("2014-03-01", "2014-05-31")$
  map(addNDVI)$
  map(cloudMask)

# Reduce the collection to the mean of each pixel and display.
meanImage <- collection$reduce(ee$Reducer$mean())
vizParams <- list(
  bands = c("B5_mean", "B4_mean", "B3_mean"),
  min = 0,
  max = 0.5
)

Map$addLayer(
  eeObject = meanImage,
  visParams = vizParams,
  name = "mean"
)

# Load a region in which to compute the mean and display it.
counties <- ee$FeatureCollection("TIGER/2016/Counties")
santaClara <- ee$Feature(counties$filter(
  ee$Filter$eq("NAME", "Santa Clara")
)$first())

Map$addLayer(
  eeObject = santaClara,
  visParams = list(palette = "yellow"),
  name = "Santa Clara"
)

# Get the mean of NDVI in the region.
mean <- meanImage$select("nd_mean")$reduceRegion(
  reducer = ee$Reducer$mean(),
  geometry = santaClara$geometry(),
  scale = 30
)

# Print mean NDVI for the region.
cat("Santa Clara spring mean NDVI:", mean$get("nd_mean")$getInfo())

## Santa Clara spring mean NDVI: 0.4650797