Datos para Junin

library(sf)
library(tidyverse)

mapaperu <- read_sf("D:/UC/CURSOS UC - 2019 I/Seminario de investigacion Departamental/MAPAS QGIS/PERU/Basica/limites_politicos/BAS_LIM_DEPARTAMENTO.shp")
plot(mapaperu["NOMBDEP"])

mapaperu <- mapaperu %>% mutate(km2 = HECTARES/100)

Junin <- mapaperu %>% dplyr::filter(NOMBDEP %in% ("JUNIN"))

ggplot() + geom_sf(data = Junin) + theme_bw()

library(leaflet)
Junin <- Junin %>% st_transform (crs = 4326)
leaflet(Junin) %>% addTiles() %>% addPolygons()

## Para obtener una imagen mas nitida, aqui nos dara la alongitud y latitud, lo cual es reemplazada en la funcion getData. 
st_centroid(Junin)

## Simple feature collection with 1 feature and 5 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -74.87737 ymin: -11.54012 xmax: -74.87737 ymax: -11.54012
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs
## # A tibble: 1 x 6
##   NOMBDEP COUNT FIRST_IDDP HECTARES    km2              geometry
##   <chr>   <dbl> <chr>         <dbl>  <dbl>          <POINT [°]>
## 1 JUNIN     123 12         4399729. 43997. (-74.87737 -11.54012)

library(raster)
library(rgdal)
Prec <- getData("worldclim", var = "prec", res=0.5, lon=-74.8773, lat=-11.54012)
Prec_Junin <- crop(Prec, Junin)
Prec_Junin <- Prec_Junin <- mask(Prec_Junin, Junin)
PPAnual_Junin<- do.call("sum", unstack(Prec_Junin))
plot(PPAnual_Junin)

Futuro_CC_2070_85 <- getData("CMIP5", var = "prec", res = 0.5, lon=-74.8773, lat=-11.54012, year= 70, model= "CC", rcp = 85)
Futuro_PJunin <- crop(Futuro_CC_2070_85, Junin)
Futuro_PJunin <- Futuro_PJunin <- mask(Futuro_PJunin, Junin)
PPFuturo_Junin<- do.call("sum", unstack(Futuro_PJunin))
plot(PPFuturo_Junin)

Diferencia <- (PPFuturo_Junin - PPAnual_Junin)                      
plot(Diferencia)

APPeru <- read_sf("D:/UC/CURSOS UC - 2019 I/Seminario de investigacion Departamental/MAPAS QGIS/PERU/Areas_protegidas/ap_peru.shp") 

APPeru$Hectareas = as.numeric(st_area(APPeru)/10000)
plot(APPeru["Hectareas"])

APJunin = st_intersection(APPeru, Junin)
plot(APJunin["NAME"])

APJunin2 = st_intersection(Junin, APPeru)
plot(APJunin2["NAME"])

ggplot() + geom_sf(data = Junin) + geom_sf(data = APJunin, aes(fill = IUCN_CAT))

AHJunin = st_intersection(APPeru, Junin)
plot(APJunin["REP_AREA"])

ggplot() + geom_sf(data = Junin) + geom_sf(data = AHJunin, aes(fill = IUCN_CAT))