Informe 2 – Analisis de Raster
Carlos Luis Mora Cañas
2023-05-22
Taller 2
A continuación, se realizara el análisis del comportamiento de la temperatura y precipitación de todo el mundo, donde se realizara inicialmente la visualización por meses y al final se realizara la identificación de zonas donde es más optimo realizar siembras por la temperatura y precipitación acumulada en el año
Cargue de archivos por mes:
precipitacion_enero=raster('~/R/Analisis Geoespacial/Semana 2/YDRAY-data_global-copy/data_global copy/prec/current_prec_01.asc')
precipitacion_junio=raster('~/R/Analisis Geoespacial/Semana 2/YDRAY-data_global-copy/data_global copy/prec/current_prec_06.asc')
temp_enero=raster('~/R/Analisis Geoespacial/Semana 2/YDRAY-data_global-copy/data_global copy/tmean/current_tmean_01.asc')
temp_junio=raster('~/R/Analisis Geoespacial/Semana 2/YDRAY-data_global-copy/data_global copy/tmean/current_tmean_06.asc')Plot de informacion cargada:
# visualizacion de los meses cargados
plot(temp_enero)
plot(temp_junio)
#visualizacion de los dos meses cargados paratemperatura y precipitacion:
both=stack(temp_enero,temp_junio)
plot(both)
both1=stack(precipitacion_enero,precipitacion_junio)
plot(both1)
Carga en Batch de todos los Raster
archivos=list.files("~/R/Analisis Geoespacial/Semana 2/YDRAY-data_global-copy/data_global copy/prec/",full.names = TRUE)
precipitacion=stack(archivos)
precipitacion## class : RasterStack
## dimensions : 300, 720, 216000, 12 (nrow, ncol, ncell, nlayers)
## resolution : 0.5, 0.5 (x, y)
## extent : -180, 180, -60, 90 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : current_prec_01, current_prec_02, current_prec_03, current_prec_04, current_prec_05, current_prec_06, current_prec_07, current_prec_08, current_prec_09, current_prec_10, current_prec_11, current_prec_12archivos1=list.files("~/R/Analisis Geoespacial/Semana 2/YDRAY-data_global-copy/data_global copy/tmean/",full.names = TRUE)
temperaturas=stack(archivos1)
temperaturas## class : RasterStack
## dimensions : 300, 720, 216000, 12 (nrow, ncol, ncell, nlayers)
## resolution : 0.5, 0.5 (x, y)
## extent : -180, 180, -60, 90 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : current_tmean_01, current_tmean_02, current_tmean_03, current_tmean_04, current_tmean_05, current_tmean_06, current_tmean_07, current_tmean_08, current_tmean_09, current_tmean_10, current_tmean_11, current_tmean_12Visualizar y modificacion del mes de todos los Raster de temperaturas y precipitacion
names(temperaturas)=month.name
names(precipitacion)=month.name
plot(temperaturas)
plot(precipitacion)
Probemos ahora la visualización con la libreria rasterVis la cual estandariza la misma simbologia de los mapas lo cual es ideal para comparar los meses:
# visualizaciones:
levelplot(temp_enero)
levelplot(temp_enero, par.settings=BuRdTheme,
panel=panel.levelplot.raster, interpolate=TRUE)
#visualizacion y estandarizacion de todo el año de Raster Temperatura
levelplot(temperaturas, panel = panel.levelplot.raster,
col.regions = hcl.colors, cuts = 30, interpolate = TRUE)
#visualizacion y estandarizacion de todo el año de Raster Precipitacion
levelplot(precipitacion, par.settings=BTCTheme,
panel=panel.levelplot.raster, interpolate=TRUE)
#Visualizacion con Leaflet
#Temperatura
# Se evidencia que los Raster tienen una granularidad muy grande con respecto a la grilla que cubre.
crs(temp_enero)='+init=EPSG:4326'
leaflet() %>%addTiles() %>%addRasterImage(temp_enero,opacity = 0.5)#Precipitacion
crs(precipitacion_enero)='+init=EPSG:4326'
leaflet() %>%addTiles() %>%addRasterImage(precipitacion_enero,opacity = 0.5)Analisis y Limpieza de datos
#Primero que todo conocer como vienen los datos:
summary(precipitacion)## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (46.3% of all cells)## January February March April May June July August September
## Min. 0 0 0 0 0 0 0 0 0
## 1st Qu. 9 9 11 14 14 14 19 20 19
## Median 23 21 24 28 34 43 52 53 44
## 3rd Qu. 61 59 65 63 68 78 88 89 81
## Max. 828 698 569 707 759 1471 1728 1185 903
## NA's 150697 150697 150697 150697 150697 150697 150697 150697 150697
## October November December
## Min. 0 0 0
## 1st Qu. 17 13 11
## Median 36 30 25
## 3rd Qu. 72 68 64
## Max. 916 802 705
## NA's 150697 150697 150697summary(temperaturas)## Warning in .local(object, ...): summary is an estimate based on a sample of 1e+05 cells (46.3% of all cells)## January February March April May June July August September
## Min. -51 -46 -42 -34 -19 -11 -7 -9 -16
## 1st Qu. -21 -19 -12 -4 3 9 11 10 5
## Median -2 -1 3 8 13 16 18 17 14
## 3rd Qu. 21 22 22 23 23 24 24 24 24
## Max. 33 33 33 33 35 37 39 38 35
## NA's 150697 150697 150697 150697 150697 150697 150697 150697 150697
## October November December
## Min. -26 -40 -48
## 1st Qu. -2 -11 -18
## Median 9 3 0
## 3rd Qu. 23 22 21
## Max. 32 32 32
## NA's 150697 150697 150727# Se observa que existen muchos NA y se procede a eliminarlos
precipitacion1=na.omit(precipitacion)
precipitacion1## class : RasterStack
## dimensions : 300, 720, 216000, 12 (nrow, ncol, ncell, nlayers)
## resolution : 0.5, 0.5 (x, y)
## extent : -180, 180, -60, 90 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : January, February, March, April, May, June, July, August, September, October, November, Decembertemperaturas1=na.omit(temperaturas)
temperaturas1## class : RasterStack
## dimensions : 300, 720, 216000, 12 (nrow, ncol, ncell, nlayers)
## resolution : 0.5, 0.5 (x, y)
## extent : -180, 180, -60, 90 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : January, February, March, April, May, June, July, August, September, October, November, Decembernames(temperaturas1)=month.name
names(precipitacion1)=month.name
# R calcula estas estadisticas utilizando todas las celdas en el raster, esto puede tomar una muestra aleatorea y calcular a partir de esto. Para forzar el calculo en mas valores, o todos los valores. Pero se observa que al final no varian las estadisticas.
summary(precipitacion1, maxsamp = ncell(precipitacion1))## January February March April May June July August September
## Min. 0 0 0 0 0 0 0 0 0
## 1st Qu. 9 9 11 14 14 14 19 20 19
## Median 23 21 24 28 34 43 52 53 44
## 3rd Qu. 61 59 65 62 67 78 89 89 81
## Max. 828 698 609 707 759 1471 1728 1232 903
## NA's 151051 151051 151051 151051 151051 151051 151051 151051 151051
## October November December
## Min. 0 0 0
## 1st Qu. 17 13 11
## Median 36 30 25
## 3rd Qu. 72 68 64
## Max. 919 802 705
## NA's 151051 151051 151051summary(temperaturas1, maxsamp = ncell(temperaturas1))## January February March April May June July August September
## Min. -51 -46 -44 -35 -20 -11 -8 -9 -16
## 1st Qu. -21 -19 -12 -4 3 9 11 10 5
## Median -2 -1 3 8 13 16 18 17 14
## 3rd Qu. 21 22 23 23 23 24 24 24 24
## Max. 33 33 33 33 35 38 39 38 35
## NA's 151051 151051 151051 151051 151051 151051 151051 151051 151051
## October November December
## Min. -27 -40 -48
## 1st Qu. -2 -11 -18
## Median 9 3 0
## 3rd Qu. 23 22 21
## Max. 32 32 32
## NA's 151051 151051 151076================================================================================
Tambien pasamos los Raster a Dataframe para realizar mas analisis
# pasamos el Raster a Dataframe
df <- raster::as.data.frame(precipitacion, xy = TRUE)
str(df)## 'data.frame': 216000 obs. of 14 variables:
## $ x : num -180 -179 -179 -178 -178 ...
## $ y : num 89.8 89.8 89.8 89.8 89.8 ...
## $ January : num NA NA NA NA NA NA NA NA NA NA ...
## $ February : num NA NA NA NA NA NA NA NA NA NA ...
## $ March : num NA NA NA NA NA NA NA NA NA NA ...
## $ April : num NA NA NA NA NA NA NA NA NA NA ...
## $ May : num NA NA NA NA NA NA NA NA NA NA ...
## $ June : num NA NA NA NA NA NA NA NA NA NA ...
## $ July : num NA NA NA NA NA NA NA NA NA NA ...
## $ August : num NA NA NA NA NA NA NA NA NA NA ...
## $ September: num NA NA NA NA NA NA NA NA NA NA ...
## $ October : num NA NA NA NA NA NA NA NA NA NA ...
## $ November : num NA NA NA NA NA NA NA NA NA NA ...
## $ December : num NA NA NA NA NA NA NA NA NA NA ...# Se observan muchos NA y se eliminan
df1=na.omit(df)
head(df1)## x y January February March April May June July August September
## 9569 -75.75 83.25 3 3 3 3 4 6 14 13 11
## 9570 -75.25 83.25 3 3 3 3 4 6 14 13 11
## 9571 -74.75 83.25 3 3 3 3 4 6 14 13 11
## 9572 -74.25 83.25 3 3 3 3 4 6 14 13 11
## 9573 -73.75 83.25 3 3 3 3 4 6 14 13 11
## 9574 -73.25 83.25 3 3 3 3 4 6 14 14 11
## October November December
## 9569 8 5 4
## 9570 8 5 4
## 9571 8 5 4
## 9572 8 5 4
## 9573 8 5 4
## 9574 8 5 4#volver un dataframe un spatial dataframe
df1_map=SpatialPointsDataFrame(coords = df1[,1:2],data =df1 )
df1_map## class : SpatialPointsDataFrame
## features : 64949
## extent : -179.75, 179.75, -56.75, 83.25 (xmin, xmax, ymin, ymax)
## crs : NA
## variables : 14
## names : x, y, January, February, March, April, May, June, July, August, September, October, November, December
## min values : -179.75, -56.75, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
## max values : 179.75, 83.25, 828, 698, 609, 707, 759, 1471, 1728, 1232, 903, 919, 802, 705# ahora se le agrega el sistema geodésico de coordenadas geográficas
crs(df1_map)='+init=EPSG:4326'
#grafica interactiva con leaflet
leaflet(df1_map) %>%addTiles() %>%addCircleMarkers(radius = 0.1,clusterOptions = markerClusterOptions())#Mapa visualizado de otra manera
spplot(df1_map[,5])
# visualizacion con ggplot2
ggplot(df1, aes(x = x, y = y, fill = January)) +
geom_raster() +
scale_fill_gradientn(colors = terrain.colors(10))## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.
Analisis por histogramas para validar cual es la mejor forma de identificar las zonas mas propensas para cultimos
#histograma todo el año - no se evidencia bien entonces se realiza a nivel de mes
suppressWarnings(hist(precipitacion, breaks = 30, main = "Histograma Precipitacion", xlab = "Valores", ylab = "Frecuencia", col = "blue"))
suppressWarnings(hist(temperaturas, breaks = 30, main = "Histograma Temperatura", xlab = "Valores", ylab = "Frecuencia", col = "blue"))
## otra manera de visualizar el histograma de enero
ggplot(df1, aes(x = January)) +
geom_histogram(fill = "steelblue", color = "white", bins = 50) +
labs(title = "Histograma de datos", x = "Valores", y = "Frecuencia")
# con lo anterior podemos evidenciar y sacar concluciones de cuales serian las metricas y sus respectivos valores.como primer analisis se identifico la en la precipitacion la tendencia se mantiene a lo largo del año, por eso decidimos tomar las precipitaciones mayores a 150 y menores a 300 suponiendo que estos son buenos rangos de agua acumulada durante un mes.
Por otro lado con respecto a la temperatura, esta si tiene variaciones durante el año, por lo cual se estimo que la temperatura promedio del todo el año estaria entre 20 a 30 grados centigrados.
================================================================================
ahora vamos a empezar a jugar con todas las posibles ajustes para realizar diferentes indicadores, unos van a ser individuales y otros agrupados para ver cual seria la mejor opcion.
#condicionales
# filtrar rangos y sacarlo en binario para temperatura
binario_temperatura = temperaturas>20&temperaturas<30
plot(binario_temperatura[[1:4]])
# ahora sumamos para identificar donde hay mejores temperaturas del año
ind1= sum(binario_temperatura)
plot(ind1)
# ahora lo pasamos a tasa para tener mejores rangos de 0 a 100
ind1_temp=(ind1*100)/12
plot(ind1_temp)
# filtrar rangos y sacarlo en binario para precipitacion
binario_precipitacion = precipitacion>150&precipitacion<300
plot(binario_precipitacion[[1]])
# ahora sumamos para identificar donde hay mejores precipitacion del año
ind2= sum(binario_precipitacion)
plot(ind2)
# ahora lo pasamos a tasa para tener mejores rangos de 0 a 100
ind2_pre=(ind2*100)/12
plot(ind2_pre)
# intento combinado filtrar rangos y sacarlo en binario para precipitacion y temperatura
binario_combinado= precipitacion>150&precipitacion<300&temperaturas>=20&temperaturas<30
plot(binario_combinado[[1:4]])
# ahora vamos a sumar el binario combinado para ver la distribucion espacial:
ind3= sum(binario_combinado)
plot(ind3)
# ahora lo pasamos a tasa para tener mejores rangos de 0 a 100
ind3_comb=(ind2*100)/12
plot(ind3_comb)
# visualizacion agrupada
indicadores_bi=stack(ind1_temp,ind2_pre,ind3_comb)
plot(indicadores_bi)
# dejamos solo donde se quiere realizar el analisis
# mejores zonas con temperatura
ind1_tempa=ind1_temp
ind1_tempa[which(ind1_tempa[]<70)]=NA
plot(ind1_tempa)
crs(ind1_tempa)='+init=EPSG:4326'
leaflet() %>%addTiles() %>%addRasterImage(ind1_tempa,opacity = 0.5)# mejores zonas con precipitacion
ind2_prea=ind2_pre
ind2_prea[which(ind2_prea[]<60)]=NA
plot(ind2_prea)
crs(ind2_prea)='+init=EPSG:4326'
leaflet() %>%addTiles() %>%addRasterImage(ind2_prea,opacity = 0.5)bin_com_f=sum(ind1_tempa,ind2_prea)
bin_com_final=(bin_com_f/200)
plot(bin_com_final)
crs(bin_com_final)='+init=EPSG:4326'
leaflet() %>%addTiles() %>%addRasterImage(bin_com_final,opacity = 0.5)# intento combiar de manera simultanea los rangos para sacarlo en binario para precipitacion y temperatura
binario_combinado_f= ind2_pre>70&ind1_temp>70
plot(binario_combinado_f)
crs(binario_combinado_f)='+init=EPSG:4326'
leaflet() %>%addTiles() %>%addRasterImage(binario_combinado_f,opacity = 0.5)# Se evidencia que aunque dan diferentes espacalas las zonas se mantienen como las mejores para realizar siembras
indicadores_f=stack(bin_com_final,binario_combinado_f)
plot(indicadores_f)
Concluciones
library(jpeg)
imagen <- readJPEG('~/R/Mapamundia.jpg')
# Mostrar la imagen
plot(0, 0, type = "n", xlim = c(0, 1), ylim = c(0, 1), xlab = "", ylab = "", asp = 1)
rasterImage(imagen, 0, 0, 1, 1)
Concluciones:
Se identifica que las zonas con mejores temperaturas y precipitaciones durante todo el año seria en la amazonia entre colombia y brasil en latinoamerica y la otra parte seria en Indonesia