Skip to content

Instantly share code, notes, and snippets.

@viciana
Created April 4, 2017 16:30
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save viciana/5f1ec6c00354b212694e6d8659bfa1eb to your computer and use it in GitHub Desktop.
Save viciana/5f1ec6c00354b212694e6d8659bfa1eb to your computer and use it in GitHub Desktop.
Ejemplo de descarga de datos climaticos
# require(devtools)
# install_github('SevillaR/Andaclima')
require(Andaclima)
rm(list = ls())
stations <- getAndalusia_ACS()
metainfo <- getMetaData(provincia = stations$province.code,
estacion = stations$station.code ,
nombre_estacion = stations$station.name )
# save(stations, metainfo, file='010_listaEstaciones.RData')
load('010_listaEstaciones.RData')
#' Convierte cadena de caracter en grados, minutos y segundo
#' en grados decimales
Char2ll <- function(e) {
is.W <- gsub("^[1-90 º']+","",e) %in% c("W","N")
e <- gsub(" [NSEW]$",'',e)
ee <- strsplit(e, "[º' ]+")
ee<-as.numeric((ee[[1]]))
ee <- sum(ee * c(1,1/60,1/3600))
if (!is.W) ee<- ee*-1
return(ee)
}
Char2ll("30º 12' 15'' W")
## Coordenada UTM
metainfo$x <- as.numeric(metainfo$x)
metainfo$y <- as.numeric(metainfo$y)
## Coordenada Latitud y Longitud
metainfo$longitud <- sapply(metainfo$longitud, Char2ll )
metainfo$longitud <- sapply(metainfo$longitud, Char2ll )
start <- '15-07-2015'
end <- '17-07-2015'
dclima <- data.frame()
i <- 1
for (i in 1:nrow(stations)) {
cat(i, '/', nrow(stations),
': Downloading data from station ', stations$station.code[i],
' at province ', stations$province.code[i], '\n', sep ='')
dclima <- rbind(dclima, readIFAPA(stations$province.code[i],
stations$station.code[i],
start = start , end =end))
}
## Return data
# save(dclima, file='dclima')
load('dclima.RData')
### Une datos temperatura con metainfo
merge(dclima,metainfo, by=c("c_provincia","estacion")) -> dclima
###----
## Interpolación espacial ...
###
library(ggplot2)
library(gstat)
library(sp)
library(maptools)
subset(dclima, FECHA==as.Date('20150715', "%Y%m%d")) -> cdia
### Convierte en SpatialPoinsDataFrame
coordinates(cdia) = ~x + y
plot(cdia)
####--------------------
range(cdia$x) # 123271.7 608847.0 21.86411 28.03960
range(cdia$y) # 4019310 4263000 57.77708 59.50061
# Define the grid extent:
x.range <- c( 123271, 608847) # min/max longitude of the interpolation area
y.range <- c(4019310, 4263000) # min/max latitude of the interpolation area
# Create a data frame from all combinations of the supplied vectors or factors.
# See the description of the return value for precise details of the way this is done.
# Set spatial coordinates to create a Spatial object. Assign gridded structure:
grd <- expand.grid(x = seq(from = x.range[1], to = x.range[2], by = 10000),
y = seq(from = y.range[1],to = y.range[2], by =10000))
# expand points to grid
class(grd)
coordinates(grd) <- ~x + y
gridded(grd) <- TRUE # Transforma en un GRID
# Plot the weather station locations and interpolation grid:
plot(grd, cex = 1.5, col = "grey" )
points(cdia, pch = 19, col = "red", cex = .4)
## Interpolate surface and fix the output:
idw <- idw(formula = TMed ~ 1,
locations = cdia,
newdata = grd) # apply idw model for the data
idw.output = as.data.frame(idw) # output is defined as a data table
class(idw.output)
names(idw.output)[1:3] <- c("long", "lat", "var1.pred") # give names to the modelled variables
head(idw.output)
# Plot the results:
ggplot() + geom_tile(data = idw.output, aes(x = long, y = lat, fill = var1.pred)) +
geom_point(data = metainfo,
aes(x = x, y = y),
shape = 19, colour = "red", size=.5)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment