# Function to download daily met data, format it and save it on the disk
#### AEMET
downloadAEMETcurrentday <- function(api, daily = TRUE, verbose=TRUE){
# # Utilitary functions
# nonUTF8 = "\u00D1\u00C0\u00C1\u00C8\u00C9\u00D2\u00D3\u00CC\u00CD\u00DC\u00CF"
# cname.func <- function(x){
# regmatches(x,gregexpr('(?<=\\n\\s{2}\\")[[:print:]]+(?=\\"\\s\\:)', x, perl = T))[[1]]
# }
# value.func <- function(x){
# value <- regmatches(x,gregexpr(paste0("(?<=\\:\\s)([[:print:]]|[",nonUTF8,"])*(?=\\n)"), x, perl = T))[[1]]
# value <- gsub('\\",', "", value)
# value <- gsub('\\"', "", value)
# value <- gsub(',', "", value)
# }
#
#
# # set options
# h = new_handle()
# handle_setheaders(h, "Cache-Control" = "no-cache", api_key=api)
# handle_setopt(h, ssl_verifypeer=FALSE)
# handle_setopt(h, customrequest="GET")
#
# url <- "https://opendata.aemet.es/opendata/api/observacion/convencional/todas"
# # get data url
# urldata_raw <- curl_fetch_memory(url, h)$content
# urldata_string <- value.func(rawToChar(urldata_raw))
#
# if(urldata_string[2]=="401"){
# stop("Invalid API key. (API keys are valid for 3 months.)")
# }
#
# urldata <- urldata_string[3]
#
# if(verbose)cat("Downloading hourly data from all available stations")
# data_raw <- curl_fetch_memory(urldata, h)$content
#
# if(verbose)cat("\nFormating data")
# data_string <- rawToChar(data_raw)
# #Add local encoding information to data_string
# # Encoding(data_string) <-"latin1"
# enclocal <- l10n_info()
# if(enclocal[[2]]) Encoding(data_string) <-"UTF-8"
# else if(enclocal[[3]]) Encoding(data_string) <-"latin1"
# data_string <- strsplit(data_string,"}\\,\\s{1}\\{")[[1]]
# cname <- lapply(data_string,FUN = cname.func)
# value <- lapply(data_string,FUN = value.func)
# value <- mapply(FUN = function(x,y){names(x) <- y;return(x)}, x = value,y = cname)
# unique_cname <- cname[[which.max(sapply(cname,FUN = length))]]
# value <- mapply(FUN = function(x,y){x[y]}, x = value, y = list(unique_cname))
#
# data_df <- data.frame(matrix(t(value), ncol = length(unique_cname), dimnames = list(NULL, unique_cname)),
# stringsAsFactors = F)
if(verbose)cat("Downloading hourly data from all available stations")
apidest = "/api/observacion/convencional/todas"
data_df = .get_data_aemet(apidest, api)
if(verbose)cat("\nFormating data")
varnames <-c("idema", "lon","lat", "ubi", "alt", "fint", "ta", "tamin", "tamax", "prec", "hr", "dv", "vv")
data_df <- data_df[,varnames]
numvar <- c("lon","lat","alt","ta", "tamin", "tamax", "prec", "hr", "dv", "vv")
data_df[,numvar] <- sapply(data_df[,numvar],as.numeric)
data_df$fint <- as.POSIXlt(sub("T", " ",data_df$fint), format = "%Y-%m-%d %H:%M:%S")
if(daily){
if(verbose)cat("\nAggregating hourly data to 24h-scale\n")
options(warn=-1)
data_agg <- aggregate(data_df[,numvar],list(idema = data_df$idema, ubi = data_df$ubi),
function(x){mean<-mean(x,na.rm=T);min<-min(x,na.rm=T);max<-max(x,na.rm=T);sum<-sum(x,na.rm=T)
return(c(mean=mean,min=min,max=max,sum=sum))})
# wind direction
dv_agg <- aggregate(list(dv = data_df$dv),list(idema = data_df$idema, ubi = data_df$ubi),
function(dvvec){
y = sum(cos(dvvec*pi/180), na.rm=TRUE)/length(dvvec)
x = sum(sin(dvvec*pi/180), na.rm=TRUE)/length(dvvec)
dv = (180/pi)*atan(y/x)
dv[dv<0] <- dv[dv<0]+360
return(dv)
})
options(warn=0)
data_df <- data.frame(ID = as.character(data_agg$idema), name = data_agg$ubi,
long = data_agg$lon[,"mean"],lat = data_agg$lat[,"mean"], elevation = data_agg$alt[,"mean"],
MeanTemperature = data_agg$ta[,"mean"], MinTemperature = data_agg$ta[,"min"], MaxTemperature = data_agg$ta[,"max"],
Precipitation = data_agg$prec[,"sum"], WindSpeed = data_agg$vv[,"mean"], WindDirection = dv_agg$dv,
MeanRelativeHumidity = data_agg$hr[,"mean"], MinRelativeHumidity = data_agg$hr[,"min"], MaxRelativeHumidity = data_agg$hr[,"max"])
data_df <- as.data.frame(lapply(data_df,function(x){
x. <- x
if(is.numeric(x.))x.[is.nan(x.)|is.infinite(x.)] <- NA
return(x.)
}))
data_sp <- SpatialPointsDataFrame(coords = data_df[,c("long", "lat")],
data = data_df[,which(!colnames(data_df) %in% c("long", "lat", "name", "ID"))],
proj4string = CRS("+proj=longlat"))
row.names(data_sp) <- data_df$ID
return(data_sp)
}else{
if(verbose)cat("\nHourly results are returned\n")
colnames(data_df) <- c("ID", "long", "lat", "name", "elevation", "date",
"MeanTemperature", "MinTemperature", "MaxTemperature",
"Precipitation", "MeanRelativeHumidity", "WindDirection", "WindSpeed")
return(data_df)
}
}
#### SMC
# download the variables metadata
downloadSMCvarmetadata <- function(api, type = "current"){
type = match.arg(type, c("current", "historical"))
if(type=="current") {
apidest <- "/variables/mesurades/metadades"
} else if(type == "historical") {
apidest <- "/variables/estadistics/diaris/metadades"
}
data <- .get_data_smc(apidest,api)
rownames(data) <- data$codi
return(data)
}
# download the met data
downloadSMCcurrentday <- function(api, daily_meteoland=TRUE, variable_code=NULL, station_id=NULL, date = Sys.Date(), verbose=TRUE){
if(daily_meteoland) variable_code <- c(32, 33, 35, 36, 46, 47)
else if(is.null(variable_code)) stop("variable_code must be specified")
variable_names = SMCvarcodes[as.character(variable_code),"nom"]
if(verbose)cat("Downloading hourly data from all available stations\n")
date_split <- strsplit(as.character(date), split = "-")[[1]]
# download variable per variable
for(i in 1:length(variable_code)){
apidest <- paste("/variables/mesurades", variable_code[i], date_split[1], date_split[2], date_split[3], sep = "/")
if(!is.null(station_id)){apidest<-paste0(apidest,"?codiEstacio=", station_id)}
data_list <- .get_data_smc(apidest, api)
data_list$variables <- sapply(data_list$variables, FUN = function(x)x$lectures)
data_i <- data.frame()
for(j in 1:length(data_list$codi)){
data_j <- data_list$variables[[j]][,c("data", "valor")]
data_j$ID <- data_list$codi[[j]]
data_i <- rbind(data_i,data_j)
}
colnames(data_i) <- c("date", as.character(variable_code[i]), "ID")
if(i == 1) {data <- data_i}else{data <- merge(data, data_i, all=T)}
}
if(verbose)cat("\nFormating data\n")
colsel <- !colnames(data) %in% c("date", "ID")
colnames(data)[colsel] <- SMCvarcodes[colnames(data)[colsel], "nom"]
data$date <- sub("T", " ", data$date)
data$date <- sub("Z", "", data$date)
data$date <- as.POSIXlt(sub("T", " ",data$date), format = "%Y-%m-%d %H:%M")
if(daily_meteoland){
if(verbose)cat("\nDownloading station info\n")
SMCstation_sp = downloadSMCstationlist(api, date = date)
if(verbose)cat("\nAggregating hourly data to 24h-scale\n")
options(warn=-1)
data$date <- as.Date(data$date)
numvar <- !colnames(data) %in% c("date", "date", "ID")
data_agg <- aggregate(data[,numvar],list(ID = data$ID),
function(x){mean<-mean(x,na.rm=T);min<-min(x,na.rm=T);max<-max(x,na.rm=T);sum<-sum(x,na.rm=T)
return(c(mean=mean,min=min,max=max,sum=sum))})
# wind direction
dv_agg <- aggregate(list(dv = data[,variable_names[6]]),
list(ID = data$ID),
function(dvvec){
y = sum(cos(dvvec*pi/180), na.rm=TRUE)/length(dvvec)
x = sum(sin(dvvec*pi/180), na.rm=TRUE)/length(dvvec)
dv = (180/pi)*atan(y/x)
dv[dv<0] <- dv[dv<0]+360
return(dv)
})
options(warn=0)
data_df <- data.frame(ID = data_agg$ID, name = SMCstation_sp@data[data_agg$ID,"name"],
long = SMCstation_sp@coords[data_agg$ID,"long"],
lat = SMCstation_sp@coords[data_agg$ID,"lat"],
elevation = SMCstation_sp@data[data_agg$ID,"elevation"],
MeanTemperature = data_agg[[variable_names[1]]][ ,"mean"],
MinTemperature = data_agg[[variable_names[1]]][ ,"min"],
MaxTemperature = data_agg[[variable_names[1]]][ ,"max"],
Precipitation = data_agg[[variable_names[3]]][ ,"sum"],
WindSpeed = data_agg[[variable_names[5]]][ ,"mean"],
WindDirection = dv_agg$dv,
MeanRelativeHumidity = data_agg[[variable_names[2]]][ ,"mean"],
MinRelativeHumidity = data_agg[[variable_names[2]]][ ,"min"],
MaxRelativeHumidity = data_agg[[variable_names[2]]][ ,"max"],
Radiation = data_agg[[variable_names[4]]][ ,"mean"])
data_df <- as.data.frame(lapply(data_df,function(x){
x. <- x
if(is.numeric(x.))x.[is.nan(x.)|is.infinite(x.)] <- NA
return(x.)
}))
data_df$Radiation = data_df$Radiation*(3600*24)/1e6 # From W/m2 to MJ/m2
data_sp <- SpatialPointsDataFrame(coords = data_df[,c("long", "lat")],
data = data_df[,which(!colnames(data_df) %in% c("long", "lat", "name", "ID"))],
proj4string = CRS("+proj=longlat"))
row.names(data_sp) <- data_df$ID
return(data_sp)
}else{
if(verbose)cat("\nHourly results are returned\n")
return(data)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.