library(rgdal)
## Loading required package: sp
## rgdal: version: 0.8-14, (SVN revision 496)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 1.10.0, released 2013/04/24
## Path to GDAL shared files: /usr/share/gdal/1.10
## Loaded PROJ.4 runtime: Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]
## Path to PROJ.4 shared files: (autodetected)
library(sp)
setwd("/home/bartosz/R/skryptydoR/gis/")
a=read.csv("temp.txt",sep=" ", header=F)
a=a[which(a$V7>13.01 & a$V7<24.99 & a$V6>48.01 & a$V7<55.99),] # przycinamy mniej wiecej do obszaru rastra
xym=data.frame(a$V7,a$V6,a$V8,a$V5)
names(xym)=c("lon","lat","t2m","code")
rm(a)
coordinates(xym)=~lon+lat
### KONIEC CZESCI Z USTAWIANIEM KOORDYNAT
srtm <- read.asciigrid("elevation.asc", as.image=FALSE, plot.image=TRUE)
names(srtm)="elevation"
srtm$lon=coordinates(srtm)[,1]
srtm$lat=coordinates(srtm)[,2]
srtm.ov = overlay(srtm, xym) # create grid-points overlay
## Warning: '.local' is deprecated.
## Use 'over' instead.
## See help("Deprecated") and help("sp-deprecated").
xym$srtm.ov =srtm.ov$elevation # copy the slope values
lm.depth <- lm(t2m~srtm.ov+lon+lat, xym)
tmp=summary(lm.depth)
#plot(t2m~srtm.ov+lon+lat, as.data.frame(xym))
library(gstat)
# temp
xym=as.data.frame(xym)
xym=xym[which(is.na(xym$srtm.ov)==F),]
names(xym)[5]="elevation"
coordinates(xym)=~lon+lat
atmin=xym
null.vgm <- vgm(var(atmin$t2m), "Sph", sqrt(areaSpatialGrid(srtm))/4, nugget=1) # initial parameters
vgm_depth_r <- fit.variogram(variogram(t2m~elevation+lon+lat, xym), model=null.vgm)
## Warning: singular model in variogram fit
plot(variogram(t2m~elevation+lon+lat,atmin), vgm_depth_r, main="fitted by gstat")
tmin=krige(t2m~elevation+lon+lat, atmin, srtm)
## [ordinary or weighted least squares prediction]
spplot(tmin)
