Primero algunos paquetes que tuve que usar para aprender o aplicar.
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(plyr)
## -------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## -------------------------------------------------------------------------
##
## Attaching package: 'plyr'
## The following objects are masked from 'package:dplyr':
##
## arrange, count, desc, failwith, id, mutate, rename, summarise,
## summarize
library(ggplot2)
library(survival)
library(HURDAT)
library(sf)
## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(raster)
## Loading required package: sp
##
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
##
## select
library(spData)
library(spDataLarge)
library(stringr)
library(tidyr)
##
## Attaching package: 'tidyr'
## The following object is masked from 'package:raster':
##
## extract
library(lubridate)
##
## Attaching package: 'lubridate'
## The following object is masked from 'package:plyr':
##
## here
## The following object is masked from 'package:base':
##
## date
library(tmap)
Usando HURDAT como la base de datos “base”. También se juntaron ambos océanos en una sola base. Se filtraron solo las variables y años de interés. Am corresponde al mapa de las americas, es dónde interactuan los huracanes.
setwd("C:/Users/USUARIO/Documents/R/Tesis/Americas")
Am <- st_read("Americas.shp", layer = "Americas")
## Reading layer `Americas' from data source `C:\Users\USUARIO\Documents\R\Tesis\Americas\Americas.shp' using driver `ESRI Shapefile'
## Simple feature collection with 53 features and 1 field
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -186.5396 ymin: -58.49861 xmax: -12.15764 ymax: 83.6236
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
AL <- get_hurdat(basin = "AL")
AL$Basin <- "AL"
EP <- get_hurdat(basin = "EP")
EP$Basin <- "EP"
hu1 <- rbind(AL,EP) %>% dplyr::select( Key, Name, DateTime, Record, Status, Lat, Lon, Wind, Basin)%>%filter(year(DateTime) > 1969,!is.na(Lat), !is.na(Lon)) %>% mutate(Year = year(DateTime))
Hurricane <- st_as_sf(hu1, coords = c("Lon", "Lat"), crs = 4326 )
Se tuvo que transformar las coordenadas en un objeto sf para poder usar los .shp que son los mapas en R.
Year <- c(1970:2018)
el_nino <- c(0, 0, 3, 0, 0, 0, 1, 1, 0, 1, 0, 0, 4, 0, 0, 0, 2, 3, 0, 0, 0, 3, 0, 0, 2, 0, 0, 4, 0, 0, 0, 0, 2, 0, 1, 0, 1, 0, 0, 2, 0, 0, 0, 0, 1, 4, 0, 0, 1)
la_nina <- c(2, 1, 0, 3, 1, 3, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 3, 0, 0, 0, 0, 0, 0, 2, 0, 0, 3, 3, 1, 0, 0, 0, 0, 1, 0, 3, 1, 0, 3, 2, 0, 0, 0, 0, 1, 1, 0)
enzo <- cbind(Year, el_nino, la_nina)
enzo
## Year el_nino la_nina
## [1,] 1970 0 2
## [2,] 1971 0 1
## [3,] 1972 3 0
## [4,] 1973 0 3
## [5,] 1974 0 1
## [6,] 1975 0 3
## [7,] 1976 1 0
## [8,] 1977 1 0
## [9,] 1978 0 0
## [10,] 1979 1 0
## [11,] 1980 0 0
## [12,] 1981 0 0
## [13,] 1982 4 0
## [14,] 1983 0 1
## [15,] 1984 0 1
## [16,] 1985 0 0
## [17,] 1986 2 0
## [18,] 1987 3 0
## [19,] 1988 0 3
## [20,] 1989 0 0
## [21,] 1990 0 0
## [22,] 1991 3 0
## [23,] 1992 0 0
## [24,] 1993 0 0
## [25,] 1994 2 0
## [26,] 1995 0 2
## [27,] 1996 0 0
## [28,] 1997 4 0
## [29,] 1998 0 3
## [30,] 1999 0 3
## [31,] 2000 0 1
## [32,] 2001 0 0
## [33,] 2002 2 0
## [34,] 2003 0 0
## [35,] 2004 1 0
## [36,] 2005 0 1
## [37,] 2006 1 0
## [38,] 2007 0 3
## [39,] 2008 0 1
## [40,] 2009 2 0
## [41,] 2010 0 3
## [42,] 2011 0 2
## [43,] 2012 0 0
## [44,] 2013 0 0
## [45,] 2014 1 0
## [46,] 2015 4 0
## [47,] 2016 0 1
## [48,] 2017 0 1
## [49,] 2018 1 0
db <- merge(Hurricane, enzo, by.x = "Year", by.y = "Year")
db
## Simple feature collection with 46112 features and 10 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -180 ymin: 0.4 xmax: 179.9 ymax: 70.7
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
## Year Key Name DateTime Record Status Wind Basin el_nino
## 1 1970 AL011970 ALMA 1970-05-17 18:00:00 <NA> TD 25 AL 0
## 2 1970 AL011970 ALMA 1970-05-18 00:00:00 <NA> TD 25 AL 0
## 3 1970 AL011970 ALMA 1970-05-18 06:00:00 <NA> TD 25 AL 0
## 4 1970 AL011970 ALMA 1970-05-18 12:00:00 <NA> TD 25 AL 0
## 5 1970 AL011970 ALMA 1970-05-18 18:00:00 <NA> TD 25 AL 0
## 6 1970 AL011970 ALMA 1970-05-19 00:00:00 <NA> TD 25 AL 0
## 7 1970 AL011970 ALMA 1970-05-19 06:00:00 <NA> TD 25 AL 0
## 8 1970 AL011970 ALMA 1970-05-19 12:00:00 <NA> TD 25 AL 0
## 9 1970 AL011970 ALMA 1970-05-19 18:00:00 <NA> TD 25 AL 0
## 10 1970 AL011970 ALMA 1970-05-20 00:00:00 <NA> TS 35 AL 0
## la_nina geometry
## 1 2 POINT (-79 11.5)
## 2 2 POINT (-79.2 11.7)
## 3 2 POINT (-79.7 12.1)
## 4 2 POINT (-80.1 12.3)
## 5 2 POINT (-80.5 12.5)
## 6 2 POINT (-81 13)
## 7 2 POINT (-81.5 13.5)
## 8 2 POINT (-82 14)
## 9 2 POINT (-82.5 14.5)
## 10 2 POINT (-82.5 15.5)
Se agregaron las columnas de los años dónde hubo el fenómeno del niño y la niña y sus intensidades (0 - No ocurrio, 1 - baja, 2 - media, 3 - fuerte, 4 - muy fuerte). Esto se refiere al cambio de temperaturas en la superficie del mar. El niño aumenta temperaturas y la niña las baja.
hur_join <- st_join(db, Am)
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
## Warning in st_is_longlat(x): bounding box has potentially an invalid value
## range for longlat data
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
hur_join
## Simple feature collection with 46112 features and 11 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -180 ymin: 0.4 xmax: 179.9 ymax: 70.7
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## First 10 features:
## Year Key Name DateTime Record Status Wind Basin el_nino
## 1 1970 AL011970 ALMA 1970-05-17 18:00:00 <NA> TD 25 AL 0
## 2 1970 AL011970 ALMA 1970-05-18 00:00:00 <NA> TD 25 AL 0
## 3 1970 AL011970 ALMA 1970-05-18 06:00:00 <NA> TD 25 AL 0
## 4 1970 AL011970 ALMA 1970-05-18 12:00:00 <NA> TD 25 AL 0
## 5 1970 AL011970 ALMA 1970-05-18 18:00:00 <NA> TD 25 AL 0
## 6 1970 AL011970 ALMA 1970-05-19 00:00:00 <NA> TD 25 AL 0
## 7 1970 AL011970 ALMA 1970-05-19 06:00:00 <NA> TD 25 AL 0
## 8 1970 AL011970 ALMA 1970-05-19 12:00:00 <NA> TD 25 AL 0
## 9 1970 AL011970 ALMA 1970-05-19 18:00:00 <NA> TD 25 AL 0
## 10 1970 AL011970 ALMA 1970-05-20 00:00:00 <NA> TS 35 AL 0
## la_nina COUNTRY geometry
## 1 2 <NA> POINT (-79 11.5)
## 2 2 <NA> POINT (-79.2 11.7)
## 3 2 <NA> POINT (-79.7 12.1)
## 4 2 <NA> POINT (-80.1 12.3)
## 5 2 <NA> POINT (-80.5 12.5)
## 6 2 <NA> POINT (-81 13)
## 7 2 <NA> POINT (-81.5 13.5)
## 8 2 <NA> POINT (-82 14)
## 9 2 <NA> POINT (-82.5 14.5)
## 10 2 <NA> POINT (-82.5 15.5)
Se juntaron las dos bases de datos, la de el mapa de las Américas y la de las coordenadas con tiempos, velocidad de viento, nombre, etc. Aquí podemos ver ya trayectorias como en el siguiente mapa.
am_map <- tm_shape(Am) + tm_polygons()
am_map + tm_shape(hur_join %>%
filter(Year == 1975) ) +
tm_dots(col = "Name", size = "Wind" , alpha =.5)
## Warning in st_is_longlat(x): bounding box has potentially an invalid value
## range for longlat data
## Warning in st_is_longlat(x): bounding box has potentially an invalid value
## range for longlat data
Como son tantos huracanes y tantos puntos este solo correponde a los huracanes del año 1975.