rosa de viento excel

Extraccion lluvia diario chirps para una cuenca

Leer los archivos de CMORPH utilizando R

CMORPH (técnica de CPC MORPHing) este análisis de precipitación global con una resolución espacial y temporal muy alta. Esta técnica utiliza estimaciones de precipitación derivadas exclusivamente de observaciones de satélites de microondas de órbita baja, y cuyas características se transportan a través de información de propagación espacial que se obtiene completamente de los datos de IR de los satélites geoestacionarios. 

En la actualidad se incorporan las estimaciones de precipitación derivadas de los microondas pasivos a bordo del DMSP 13, 14 y 15 (SSM / I), NOAA-15, 16, 17 y 18 (AMSU-B) y AMSR-E y TMI a bordo del Aqua de la NASA. y la nave espacial TRMM, respectivamente. Estas estimaciones son generadas por los algoritmos de Ferraro (1997) para SSM / I, Ferraro et al. (2000) para AMSU-B y Kummerow et al. (2001) para TMI. Tenga en cuenta que esta técnica no es un algoritmo de estimación de precipitación, sino un medio por el cual se pueden combinar las estimaciones de los algoritmos de lluvia de microondas existentes. Por lo tanto, este método es extremadamente flexible, de modo que se pueden incorporar las estimaciones de precipitación de cualquier fuente de satélite de microondas.

Las estimaciones de precipitación están disponibles en una cuadrícula con un espaciado de 8 km (en el ecuador), la resolución de las estimaciones individuales derivadas de satélites menor que eso, y oscilan entre los 12 x 15 km más o menos, para obtener una resolución más fina,  se realizan interpolaciones de la estimaciones.

Los datos los puede descargar de:

ftp://ftp.cpc.ncep.noaa.gov/precip/CMORPH_V0.x/RAW/8km-30min/2018/

los registros históricos desde enero de 1998 a julio de 2017, los puede descargar de:

ftp://ftp.cpc.ncep.noaa.gov/precip/CMORPH_V1.0/RAW/8km-30min/

Los datos los pueden visualizar en R y guardar en Geotiff, para luego porcesar en Qgis o cualquier GIS

Comparación de la lluvia diaria de las estaciones y CMORPH

Añadir leyenda

Grado de asociación entre la lluvia diaria y CMORPH para las estaciones de la red de INSIVUMEH, durante el periodo 2001 a 2015

############################################################################
#
#       Leer archivos de CMORPH con R
#       Walter Bardales (bardaleswa@gmail.com)
#
#

# Infomracion del archivo CTL
# TITLE  Precipitation estimates
# XDEF 4948 LINEAR   0.036378335 0.072756669
# YDEF 1649 LINEAR -59.963614    0.072771377
# ZDEF   01 LEVELS 1
# TDEF 999999 LINEAR  00z01jan2005 30mn
# VARS 1
# cmorph   1  99  hourly cmorph [ mm/hr ]

# Cargar librerias
library(raster)
library(tibble)
library(stats)
library(rgdal)

# Definir directorio de trabajo
setwd("D:/Documents/chirps")

# Leer archivo binario de cmorph
cmorph.file = file("CMORPH_V1.0_RAW_8km-30min_1998010105", "rb")

# Cerrar
on.exit(close(cmorph.file))

# Cantidad de bites XDEF * YDEF * 2
bites = (4948*1649*2)

# Coordenadas iniciales CMORPH, coordenada inicial, cantidad de pixeles y tamaño de pixel
longitud = seq(0.036378335, length.out = 4948, by=0.072756669)
latitud = seq(-59.963614, length.out = 1649, by=0.072771377)

# Leer datos binarios
temporal = readBin(cmorph.file, numeric(), n = bites, size = 4, endian = "little")
temporal = temporal[seq_len(bites/2)] * 0.1

# hacer data.frame
datacmorph = as_data_frame(
  setNames(
    cbind(expand.grid(longitud, latitud), temporal),
    c("lon", "lat", "lluvia")
  )
)

# Convertir en puntos
coordinates(datacmorph) = ~lon+lat

# Crear una grilla de puntos
gridded(datacmorph) = TRUE

# Crear un raster a partir de los puntos y poner los valores como NA si son menores a 0
rascmorph = raster(datacmorph)
rascmorph[rascmorph<0 0="" na="" p="">

<0 na="" p=""># Convertir longitude 0-360 a longitude -180-180
rascmorph = rotate(rascmorph)

# ASignarle la proyección al raster
crs(rascmorph) = "+proj=longlat +datum=WGS84"

# Guardar el raster
writeRaster(rascmorph, "cmorph.tif", format="GTiff", overwrite=TRUE)

Leer los archivos del Hidroestimador Nesdis usando R

La técnica Hidroestimador es empleada para estimar lluvia a través de satélites, esta técnica fue desarrollada por la Administración Nacional Oceánica y Atmosférica / Servicio Nacional de Satélites Ambientales, Datos e Información. (NOAA/NESDIS).  

NESDIS trabaja con las imágenes del canal infrarrojo del satélite GOES e información de variables meteorológicas pronosticadas por el modelo NAM.  

• Los píxeles de lluvia y no lluvia, se separan de acuerdo a un valor construido con la media y la desviación estándar de la temperatura de brillo en un área centrada en el píxel de interés. 
• Las nubes en un determinado píxel, producen precipitación si poseen topes más fríos que la media de los píxeles circundantes. 
• La tasa de precipitación, se ajusta teniendo en cuenta la humedad del entorno, a partir de datos de humedad relativa y agua precipitable de los modelos de pronóstico numérico. 

Los datos para la región centroamericana están disponibles  en:
https://www.star.nesdis.noaa.gov/smcd/emb/ff/HEcentralAmerica.php

Esta página proporciona estimaciones de la tasa de lluvia en tiempo real:

• Los promedios de la tasa de lluvia de 1 hora y las estimaciones acumuladas de 3 horas están disponibles cada hora a partir de las 00:00 h con actualizaciones cada 15 minutos después de la hora. 
• Las estimaciones de 6 horas acumuladas se actualizan a las 00: 15Z, 06: 15Z, 12: 15Z y 18: 15Z y la acumulación de 24 horas cada 12: 15Z. 
• Las estimaciones pueden ser demasiado altas para las nubes frías y demasiado bajas para las nubes cálidas. El uso de este producto es a discreción del usuario.

Esta información puede servir para generar información en regiones donde no se tiene el monitoreo o para mejorar lo métodos de interpolación empleados en la actualidad. 

Al leer los datos en R y guardarlos como raster en formato geotif se pueden abrir en qgis u otro software gis.

############################################################################
#
#       Leer archivos de hydroestimador Nesdis con R
#       Walter Bardales (bardaleswa@gmail.com)
#
#


# Configuración del archivo Nesdis para Centro América
#DSET K20053141609
#UNDEF -1.0
#TITLE Rain Rate Estimates (mm)
#XDEF 1464 LINEAR 245 0.0382513
#YDEF 667 LINEAR 6 0.0359477
#ZDEF 01 LEVELS 1
#TDEF 1 LINEAR 16:09Z10NOV2005 1hr
#VARS 1
#rr 0 99 Goes Data
#ENDVARS

# Cargar librerias
library(raster)
library(tibble)
library(stats)

# Directorio de trabajo
setwd("D:/Documents/Nesdis")

# Leer archivo binario y almacenarlo en objeto
archivo.nesdis = file("N20182631200", "rb")

# Salir
on.exit(close(archivo.nesdis))


# Cantidad de bits del archivo XDEF * YDEF * 2 bits
bites = (667*1464*2)


# Coordenadas de inicio, informacion extraidad de XDEF y YDEF
longitud = seq(-115, length.out = 1464, by=0.0382513)
latitud = seq(6, length.out = 667, by=0.0359477)


# Leer datos binarios
temporal = readBin(archivo.nesdis, numeric(), n = bites, size = 4, endian = "little")
temporal = temporal[seq_len(bites/2)] * 0.1

# hacer data.frame de las coordenadas y cantidad de lluvia
datanesdis = as_data_frame(
  setNames(
    cbind(expand.grid(longitud, latitud), temporal),
    c('lon', 'lat', 'lluvia_mm')
  )
)


# Convertir el dataframe en puntos
coordinates(datanesdis) = ~lon+lat

# Asignar el valor True para realizar una grilla con los puntos
gridded(datanesdis) = TRUE

# Convertir los puntos a raster, valores negativos como NA
rasnesdis = raster(datanesdis)
rasnesdis[rasnesdis < 0] = NA

# ASignarle la proyección al raster
crs(rasnesdis) = "+proj=longlat +datum=WGS84"

# Guardar el raster
writeRaster(rasnesdis, "nesdis.tif", format="GTiff", overwrite=TRUE)

Interpolar lluvia usando de covariable elevación con el método Thin Plate Spline

Thin Plate Splines (Tps)

La interpolación de Thin Plate Spline corresponde a las técnicas de interpolación n-dimensional a partir de datos dispersos. La función de interpolación correspondiente es una función continuamente diferenciable (función C1) que tiene la energía (norma L2) de su segundo derivado mínimo.

El residuo de interpolación con respecto a esta interpolación puede ser 0 ó, lo que es más interesante, se puede hacer una compensación de regularidad entre el error de interpolación y la energía de la curvatura de la función.

Wahba (1990) propone un método basado en la validación cruzada para determinar la compensación óptima.

• Thin Plate Spline adapta los datos espaciados irregularmente a una superficie. 
• El parámetro de suavizado se elige mediante validación cruzada generalizada. 
• El modelo asumido es el aditivo Y = f (X) + e donde f (X) es una superficie dimensional. 
• Este es un caso especial de la estimación espacial.

Thin Plate Spline es el resultado de minimizar la suma residual de cuadrados sujetos a una restricción de que la función tiene un cierto nivel de suavidad (o penalización de rugosidad). En donde la rugosidad se cuantifica por la integral de los derivados de orden m cuadrado al cuadrado.

Para una dimensión y m = 2, la penalización de rugosidad es el cuadrado integrado de la segunda derivada de la función. Para dos dimensiones, la penalización por rugosidad es la integral de

(Dxx (f))^2 + 2 (Dxy (f)) ^2 + (Dyy (f)) ^2

Además de controlar el orden de las derivadas, el valor de m también determina el polinomio base que se ajusta a los datos. El grado de este polinomio es (m-1).

El parámetro de suavizado controla la cantidad de suavización de los datos. En la forma habitual, esto se denota con lambda, el multiplicador de Lagrange del problema de minimización. Aunque esta es una escala incómoda, lambda = 0 no corresponde con restricciones de suavidad y los datos se interpolan. lambda = infinito corresponde a simplemente ajustar el modelo de base polinomial por mínimos cuadrados ordinarios.

Este estimador se implementa pasando la función de covarianza generalizada correcta basada en funciones de base radial a la función más general Krig.  Una ventaja de esta implementación es que una vez que se crea un objeto Tps / Krig, el estimador se puede encontrar rápidamente para otros datos y parámetros de suavizado siempre que las ubicaciones permanezcan sin cambios. Esto hace que la simulación dentro de R sea eficiente.  Tps no admite actualmente el argumento de los nudos donde se puede usar un conjunto reducido de funciones básicas. Esto es principalmente para simplificar y una buena alternativa al usar nudos sería usar una covarianza válida de la familia Matern y un parámetro de gran rango.

Código utilizado en R

##################################################################

#

#       Walter Arnoldo Bardales Espinoza (bardaleswa@gmail.com)

#       julio de 2018

#

##################################################################

# Cargar librerias

library(fields)

library(raster)

library(rgdal)

library(lattice)

library(stringr)

# Cargar modelo de elevación

dem = raster("D:/Documents/GIS DataBase/DEM_CA/geotif/dem_90_gtm.tif")

# Extraer las coordenadas del DEM

coords = xyFromCell(dem,1:ncell(dem))

# Extraer los factores de las coordenadas x e y

coord_x = as.numeric(levels(factor(coords[,"x"])))

coord_y = as.numeric(levels(factor(coords[,"y"])))

# Convertir el dem en matriz y transponerlo

z = t(as.matrix(dem))

# Corregir la matriz de elevacion debido a que esta rotada 180 grados

coord_z = matrix(, nrow=length(coord_x), ncol=length(coord_y))

for(i in 1:length(coord_y)){

  j = length(coord_y)+1-i

  coord_z[,j] = z[,i]

}

#crear lista de elevación a partir del dem

elevacion = list(x=coord_x,y=coord_y,z=coord_z)

# Evaluar la lista

str(elevacion)

# Cargar los datos de lluvia

datospp = read.csv("D:/Documents/Lluvia/Lluviaxy_1980_2016.csv",sep=",")

datos = datospp

# Convertir en puntos geograficos y asignarle las coordenadas GTM

coordinates(datos) = ~x + y

gtm = CRS("+proj=tmerc +lat_0=0 +lon_0=-90.5 +y_0=0 +x_0=500000 +k=0.9998  +datum=WGS84 +units=m  +no_defs +ellps=WGS84 +towgs84=0,0,0")

proj4string(datos) = gtm

# Extraer los valores de elevacion de las estaciones a partir del DEM

altura = extract(dem,datos)

# Crear lista con los datos x, y, z y variable

xydatos = data.frame(lon=datospp$x,lat=datospp$y)

zdatos = altura

# Antes de continuar verificar que no hallan datos faltantes en las series

# Nombrar columnas

mes = c(0,0,0,1:13)

meses = c("","","","enero", "febrero", "marzo", "abril", "mayo", "junio", "julio", "agosto", "septiembre", "octubre", "noviembre", "diciembre", "Anual")

for(k in 4:16){

  ydatos = datospp[,k]

  

  precip = list(x=xydatos,elev=zdatos,y=ydatos)

  

  # Graficar las estaciones y valores de lluvia

  quilt.plot(precip$x,precip$y)

  

  # Crear el modelo de lluvia elevacion

  fit = Tps(precip$x, precip$y, Z=precip$elev)

  

  # Resumen del modelo lluvia elevacion con TPS

  summary(fit)

  

  # Graficar los resultados

  set.panel(2,2)

  plot(fit)

  

  # Intrpolar superficie del modelo creado

  set.panel()

  # Crear la lista de coordenadas a interpolar

  grid.list = list(x=elevacion$x, y=elevacion$y)

  

  #Interpolacion del modelo usando de covariable la elevación

  fit.full = predictSurface(fit, grid.list, extrap=TRUE, ZGrid = elevacion)

  # Transponer la matriz de resultados

  pp_z1 = t(fit.full$z)

  

  # Crear el raster interpolado

  ras = raster(pp_z1, xmn=extent(dem)[1],xmx=extent(dem)[2],ymn=extent(dem)[3],ymx=extent(dem)[4])

  ras = flip(ras,direction="y")

  crs(ras) = gtm

  ras[ras<0 0="" p="">

  

  # Guardaer el raster

  destino5 = paste("D:/Documents/Consultorias_2018/balance_hidrico/Raster/lluvia/prueb/prec_",str_pad(mes[k],width =2, pad="0"),".tif", sep ="")

  writeRaster(ras,filename=destino5, overwrite=TRUE)

}

Resultados de la interpolación 

Hacer un Rscript para QGIS

Vincular R con QGIS

Estimar parámetros de la curva de calibración de caudales o descarga con R

#################################################################################
#
#
# Estimación de los parámetros de l curva de descarga
# Ajusta y plotea la grafica de curva descarga siguiendo una funcion exponencial
#
#  Walter Arnoldo Bardales Espinoza
#
#################################################################################

# remover todos los objetos de almacenamiento
rm(list = ls())

# cargar libreria
require(Hmisc)
library(minpack.lm)

# leer datos de aforos y niveles
dato.aforos = read.table("D:/Documents/Curva_calibracion_de_caudales/est_Matucuy.txt", header = TRUE, sep = "\t",na.strings = "-99.9")


plot(Caudal~Nivel)

# Para usar los nombres de los campos de las columnas
attach(dato.aforos)


# Estima los valores de los parametros iniciales de la curva
# Estableciendo h0 a un valor por debajo del minimo observado
rango=range(Nivel)
h.min=rango[1]
h.max=rango[2]
ho.inicial =  h.min - 0.1*(h.max - h.min)

# modelo lineal para estimar los valores inicial de k y n
modelo.log = lm( log10(Caudal) ~ log10(Nivel - ho.inicial) )
a.inicial = 10^modelo.log$coefficients[[1]]
b.inicial = modelo.log$coefficients[[2]]

# Utilizar los coeficientes calculados inicialmente para estimar los parametros de la regresion no lineal
curva.descarga = nlsLM( Caudal ~ k*(Nivel-ho)^n, start = list(k=a.inicial,n=b.inicial, ho=ho.inicial), control = list(maxiter=500))

# Resumen de los estadisticos de la curva de descarga
summary(curva.descarga)

# Grafico de la curva de descarga o calibracion
xmin = min(Nivel)
xmax = max(Nivel)
dx = 0.01*(xmax - xmin)
x = seq(xmin,xmax,dx)
coef.model = coef(curva.descarga)
qest = coef.model[1]*(x - coef.model[3])^coef.model[2]


plot(Nivel, Caudal,
     ylim = c(0, max(Caudal)+0.2*(max(Caudal)-min(Caudal))),
     xlim = c(min(Nivel)-0.1*(rango[2]-rango[1]), max(Nivel)+0.2*(rango[2]-rango[1])),
     xlab = "Nivel (m)",
     ylab = expression("Q ("*m^3*s^{-1}*")"),
     pch = 21, bg = "black",
     main = "Curva de calibracion de caudales de la estacion Matucuy, Periodo 2002 - 2015"
)
lines(x, qest, col = "red")
text(min(Nivel)+0.1*(max(Nivel)-min(Nivel)),max(Caudal), paste0("Q = ",round(coef.model[1],3),"*(H-(",round(coef.model[3],3),"))^",round(coef.model[2],3)))
text(min(Nivel)+0.1*(max(Nivel)-min(Nivel)),0.95*max(Caudal),paste0("RMS =", round(mean((Caudal-(coef.model[1]*(Nivel - coef.model[3])^coef.model[2]))^2))))

#################################################################### Fin código

Descargar datos y scripts R

Relleno de datos de lluvia mensual de una estación con CHIRPS

Anteriormente, generamos una función para descargar los raster de lluvia estimada a partir de satélite (CHIRPS V 2.0), ahora lo que vamos hacer es extraer los datos para un punto en especifico o la ubicación de una estación climática, el periodo de tiempo de la serie de datos chirps debe ser la misma que la serie de datos de la estación (no importa que existan datos faltantes, solo especificar el datos faltante como -99.9) para que sean comparables.

Luego se creara un modelo de ajuste entre los datos estimados y los datos reales, esto con el fin de que los datos a rellenar sean acordes a la tendencia de los datos observados o pluviometricos.

La importancia de completar series de datos radica darle continuidad a la serie y que esto me permita realizar mejores análisis, crear modelos hidrológicos estables, balances hídricos, etc.

####################################################################

library(raster)

# Ruta del directorio de trabajo
setwd("D:/Documents/chirps/")

# Directorio de los datos de chirp V-2.0
ruta_chirps = "D:/Documents/chirps/raster/"

# Periodo de datos a analizar
fecha_inicial = c("1981/1/1")
fecha_final = c("2018/5/1")

# Coordenada geografica de las estacion (longitud, Latitud)
coor.x = c(-89.8106)
coor.y = c(15.6083)

# Convertir las coordenadas de las estaciones a puntos
ubicacion = data.frame(lon=coor.x, lat=coor.y)
coordinates(ubicacion) = c("lon", "lat")
proj4string(ubicacion) = CRS("+init=epsg:4326") # WGS 84

# Cargar los datos raster
lista_chirps = paste0(ruta_chirps,"chirps-v2.0.", format(seq(as.Date(fecha_inicial), as.Date(fecha_final), "month"), "%Y.%m"), ".tif")
ras_chirps = stack(lista_chirps)

# Extraer datos de lluvia de la estacion
chirps = as.vector((extract(ras_chirps, ubicacion)))

# Cargar datos pluviometricos (reales u observados) de la estacion
estacion = scan("estaciones/Cahabon_1_1981_5_2018.txt", sep = "\t", na.strings = "-99.9")
serie = estacion


# Modelo de ajuste de datos chirps a pluviometicos
model = lm(estacion~chirps)

# estadisticos del modelo de ajuste lineal
summary(model)

# coeficientes del modelo de ajuste lineal
alfa = coef(model)[[2]]
beta = coef(model)[[1]]

# Grafica
plot(estacion~chirps)
abline(b=1,a=1, col="red")

# posiciones de datos faltantes pluviometricos
faltantes = which(is.na(serie))

# Datos chirps ajustados para relleno
relleno = round(alfa*chirps[faltantes]+beta,1)

# Completacion de la serie real
serie[faltantes] = relleno

# Grafico de serie completa (negro) y serie incompleta (rojo)
plot(serie,type="l", xlab="Lluvia [mm]", main="Comparacion de serie completa e incompleta de lluvia")
lines(estacion, col="red")

# Crear tabla de fechas y serie completa
fecha = format(seq(as.Date(fecha_inicial), as.Date(fecha_final), "month"), "%Y.%m")
serie.completa = data.frame(fecha,serie)

write.csv(serie.completa, "Cahabon_rellenada.csv", row.names = FALSE)

#############################################################

Descargar
Archivo_R

Delimitacion de cuencas con qgis y grass

https://youtu.be/ytsEYNy4Xn4

Extraer datos de lluvia en una cuenca usando R

La falta de datos de lluvia en una cuenca, no es impedimento para poder realizar un balance hídrico hoy en día, ya que existen datos de hidroestimación o datos de lluvia estimada a partir de satelites o radares.  Estos datos tienen la ventaja que proporcionan la distribución espacial de la lluvia, pero su desventaja es que tienen a subestimar o sobrestimar en algunos eventos de lluvia y esto depende del tipo de lluvia; pero mediante un ajuste por regresión se puede reducir el error.
 También, pueden usarse para rellenar datos mensuales, trimestrales, anuales de lluvia en las series observadas, mediante un ajuste.

## Script para descargar datos chirps 2.0 y extraer la lluvia media de la cuenca
## Walter Bardales
## Fecha: 01/07/2018


#################################################################################

library(R.utils)

# Funcion de descarga Chirps V2.0
download.chirps2.0 = function(Ruta, fecha_inicial, fecha_final){
fecha=format(seq(as.Date(fecha_inicial), as.Date(fecha_final), "month"), format="%Y.%m")
for (i in 1:length(fecha)){
wurl=paste0("ftp://ftp.chg.ucsb.edu/pub/org/chg/products/CHIRPS-2.0/camer-carib_monthly/tifs/","chirps-v2.0.",fecha[i],".tif.gz")
download.file(wurl,paste0(Ruta,"chirps-v2.0.",fecha[i],".tif.gz"))
gunzip(paste0(Ruta,"chirps-v2.0.",fecha[i],".tif.gz"))
}
}

#################################################################################

# Cargar libreria de trabajo
library(raster)
library(rgdal)

# Directorio de trabajo
setwd("D:/Documents/chirps/")

# Cargar funcion de descarga chirps V2.0
source("Funcion_descarga_chirps.R")

# Directorio de la descarga
Ruta = "D:/Documents/chirps/raster/"

# Periodo de descarga
fecha_inicial = c("1983/5/1")
fecha_final = c("1983/8/1")

# Descargar archivos chirps y descomprimir
download.chirps2.0(Ruta,fecha_inicial, fecha_final)

# Desactivar el paquete R.utils por tener conflicto con el paquete Raster
detach("package:R.utils", unload=TRUE)

# Cargar los datos de chirps 2.0
lista_archivos = paste0("raster/","chirps-v2.0.",format(seq(as.Date(fecha_inicial), as.Date(fecha_final), "month"),"%Y.%m"),".tif")
chirps = stack(lista_archivos)

# Cargar el poligono de la cuenca en coordenadas geograficas
cuenca = readOGR("shapefile/La_Pasion.shp") # Indicar el directorio de la cuenca

# Extraer los datos
ppcuenca = t(extract(chirps, cuenca, fun=mean))
row.names(ppcuenca) = NULL # Crear serie de datos y guardar
lluvia_media = data.frame(fecha=format(seq(as.Date(fecha_inicial), as.Date(fecha_final), "month"),"%Y.%m"), Lluvia=ppcuenca)

write.csv(lluvia_media, "Lluvia_La_Pasion.csv", row.names = FALSE) # Cambiar el nombre del archivo de salida

#######################################################################

Descargar scripts

Archivo

Graficar imagenes goes con R

library(ncdf4)
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)