library(raster)
library(rgdal)
library(fields)
library(RColorBrewer)
#Directorio de trabajo
setwd("D:/Documents/Goes_R/netcdf")
# Shapefile
shp1 = readOGR("D:/Documents/GIS DataBase/Shape_WGS84/departamentos.shp")
shp2 = readOGR("D:/Documents/GIS DataBase/Shape_WGS84/Centro_America.shp")
# Proyeccion de goes
geost = CRS("+proj=geos +a=6378137. +b=6356752.31414 +lon_0=-75 +f=0.00335281068119356027 +h=35786023. +sweep=x +no_defs")
# Coordenadas geograficas
geo = CRS("+proj=longlat +datum=WGS84")
# Leer archivo netcdf y voltearlo verticalmente
goes_16 = flip(raster("OR_ABI-L2-CMIPF-M3C14_G16_s20181742345453_e20181742356220_c20181742356317.nc"),direction = "y")
crs(goes_16) = geost # Asignar la proyeccion del satelite geoestacionario
extent(goes_16) = c(-5434894.885056, 5434894.885056, -5434894.885056, 5434894.885056) # Asignarle la extension
# Recortar el area de interes (Centro america)
cgoes_16 = crop(goes_16,extent(c(-3424153, -231218.2, 791968.4, 3216755)))
# Reproyectar de GRS80 a coordenadas WGS84
cgoes_16 = projectRaster(cgoes_16, crs=geo)
# Plotear y agregar el mapa de centro america y Departamentos Guatemala
plot(cgoes_16 - 273.15, xlab="Longitud", ylab="Latitud", main="Banda 14")
plot(cgoes_16-273.15, xlab="Longitud", ylab="Latitud", main="Banda 14", xlim=c(-105, -67), ylim=c(8,30), col=rainbow(20, alpha = .6))
lines(shp2)
lines(shp1)
