<%
library(intamap)
library(Cairo)
library(sp)
library(automap)
library(rgdal)
library(maptools)
library(gstat)
impdata <- read.table('/var/www/R/Rtmp/ClimaticDataMK.txt', header = T, sep = '\t')
# variable names
datatype <- "Precipitation"
pltlabel <- "Precipitation (mm)"
mname <- "January"
pval <- 0
if(!is.null(POST)) {
datatype <- as.character(POST$Data)
if(datatype == "Precipitation") {
pltlabel <- "Precipitation (mm)"
} else if (datatype == "Temperature") {
pltlabel <- expression(paste("Temperature (",degree,"C)"))
} else if (datatype == "Relative Humidity") {
pltlabel <- "Relative Humidity (%)"
} else {
pltlabel <- "Wind Speed (m/s)"
}
mname <- as.character(POST$Month)
}
QuantileType <- function (v, p) {
v = sort(v)
m = 0
n = length(v)
j = floor((n * p) + m)
g = (n * p) + m - j
y = ifelse (g == 0, 0.5, 1)
((1 - y) * v[j]) + (y * v[j+1])
}
mkdata <- subset(impdata , impdata$DataType == datatype)
mkdata$value = mkdata[[mname]]
pval <- QuantileType(mkdata[[mname]], 0.75)
coordinates(mkdata) = ~X+Y
## Assign a Coordinate Reference System using a PROJ.4 string
prj <- CRS("+init=epsg:31277")
proj4string(mkdata) <- prj
## mkbnd <- readOGR(dsn = "https://geogis.bgsu.edu/R/mkclimate/polys", layer="mk_boundary")
mkbnd <- readShapePoly("/var/www/R/mkclimate/polys/mk_boundary.shp")
proj4string(mkbnd) <- prj
#generate some random values that serve as fixed points
value_points <- spsample(mkbnd, type="stratified", n = 400)
values <- data.frame(value = rnorm(dim(coordinates(value_points))[1], 0 ,1))
value_df <- SpatialPointsDataFrame(value_points, values)
#generate a grid that can be estimated from the fixed points
grd = spsample(mkbnd, type = "regular", n = 4000)
kr <- autoKrige(value~1, value_df, grd)
mkdata.grid = as.data.frame(kr$krige_output)
gridded(mkdata.grid)= ~x1+x2
# spplot(mkdata,"value",col.regions = bpy.colors())
#output = interpolate(observations = mkdata,
#predictionLocations = mkdata.grid,
#outputWhat = list(mean=T, variance=T, excprob = pval),
#maximumTime = 30,
#methodName = "automatic",
#optList = list(), optList = list())
#output = interpolate(observations = mkdata, mkdata.grid,
#outputWhat = list(mean = T, variance = T, excprob = pval),
#methodName = "automatic")
output = autoKrige(value~1, mkdata, mkdata.grid)
# plot(output)
# dev.off()
fnamepng <- paste(tempfile(pattern = "png_", tmpdir='/var/www/tmp'), '.png', sep='')
if (!file.exists(fnamepng)){
CairoPNG(filename=fnamepng,width=600,height=600)
par(mar=c(1,1,1,1), oma=c(3,1,0,0))
plot(output)
dev.off()
}
setContentType("text/html; charset=utf-8")
fname1 <- paste("https://geogis.bgsu.edu/tmp/", '', basename(fnamepng), sep='')
%>
Kriging
