R/interpolateDailyWeather.R

# fields, Tools for spatial data
# Copyright 2004-2011, Institute for Mathematics Applied Geosciences
# University Corporation for Atmospheric Research
# Licensed under the GPL -- www.gpl.org/licenses/gpl.html
# Modified from a version sent to me in an email by Doug Nychka (nychka@ucar.edu)
.TpsLonLat <- function(x, Y, Z) {
    x <- as.matrix(x)
    d <- ncol(x)
	m <- max(c(2, ceiling(d/2 + 0.1)))
    p <- (2 * m - d)
    if (p <= 0) {stop(" m is too small  you must have 2*m -d >0")}
    Tpscall <- match.call()
    Tpscall$cov.function <-
        "Thin plate spline radial basis functions (RadialBasis.cov) using great circle distance "
    Krig(x=x, Y=Y, Z=Z, cov.function = stationary.cov, m = m, scale.type = "range",
        outputcall = Tpscall, GCV = TRUE,
         cov.args=list(Covariance="RadialBasis",M=m, dimension=2, Distance="rdist.earth",
            Dist.args=list(miles=TRUE)) )
}

interpolateDailyWeather <- function(tableGSOD, locations, startDate, endDate, vars=c("TEMP","MAX","MIN","PRCP"), covars="ALT", stations, sqrtTr=("PRCP"), silent=TRUE) #covars should be present in the stations file
{
	require(fields)

	#checking and formatting of incoming objects
	if(is.null(colnames(locations))) {stop("object locations has no column names.")}
	if(!all(c("ID","LON","LAT","ALT") %in% colnames(locations))) stop("object locations has non-standard column names (standard is ID,LON,LAT,ALT).")

	locations <- locations[c("ID","LON","LAT","ALT")]
	
	if(!is.numeric(locations$LON)){locations$LON <- as.numeric(locations$LON); warning("locations column LON forced to numeric")}
	if(!is.numeric(locations$LAT)){locations$LAT <- as.numeric(locations$LAT); warning("locations column LAT forced to numeric")}
	if(!is.numeric(locations$ALT)){locations$ALT <- as.numeric(locations$ALT); warning("locations column ALT forced to numeric")}
	#TODO more controls
	
	#transform some variables
	if(!is.null(sqrtTr)) {tableGSOD[sqrtTr] <- sqrt(tableGSOD[sqrtTr])} #TO DO

	#prepare vectors of the days
	days <- as.Date(startDate):as.Date(endDate)
	days <- as.Date(days, origin="1970-01-01")
	monthDays <- paste(substring(as.character(days),6,7),substring(as.character(days),9,10),sep="")
	years <- as.integer(substring(as.character(days),1,4))

	#make the final table as a list: station name, year, date and weather vars
	inDaWe <- vector(mode="list",length=length(vars)+3)
	names(inDaWe) <- c("ID", "Year", "Date", vars)
	inDaWe$ID <- rep(locations$ID,times=length(days))
	inDaWe$Year <- rep(years,each=dim(locations)[1])
	inDaWe$Date <- rep(days,each=dim(locations)[1])

	ll <- length(locations[,1])
	for(i in 4:length(inDaWe)) inDaWe[[i]] <- rep(NA, times=ll*length(days))

	#loop through days and variables
	for(i in 1:length(years))
	{
		ssub <- tableGSOD[tableGSOD["MODA"]==monthDays[i] & tableGSOD["YEAR"]== years[i],]
		if(dim(ssub)[1]>5)
		{
			#add coordinates and altitude
			ssub <- cbind(ssub, stations[match(ssub$ID, stations$ID),c("LON","LAT","ALT")])
			for(Yvar in vars)
			{
				ssubV <- ssub[c(Yvar, "LON","LAT","ALT")]
				ssubV <- na.omit(ssubV)
				if(dim(ssubV)[1]>5)
				{
					#make daily interpolation and predict for each location
					model <- try(.TpsLonLat(Y=ssubV[Yvar], x=ssubV[c("LON","LAT")], Z=unlist(ssubV["ALT"]))) 
					if(!silent) surface(model)
					#model <- fastTps(Y=ssub[Yvar], x=ssub[c("LON","LAT")], Z=ssub["ALT"], lon.lat=TRUE, theta=theta)
					if(!inherits(model, "try-error")) {inDaWe[[Yvar]][(ll*(i-1)+1):(ll*i)] <- predict(model, x=locations[,c("LON","LAT")], Z=locations[,"ALT"])}
				}
			}
		}
		
		#if there are 3 or 4 data points I could do IDW for points within the convex hull of the coords
		
	}
	result <- as.data.frame(inDaWe)
	if(!is.null(sqrtTr)) 
	{
		result[sqrtTr] <- pmax(rep(0, times=dim(result)[1]), as.vector(unlist(result[sqrtTr])))
		result[sqrtTr] <- as.vector(result[sqrtTr]^2)
	}
	return(result)
}

Try the weatherData package in your browser

Any scripts or data that you put into this service are public.

weatherData documentation built on May 2, 2019, 5:47 p.m.