[R-sig-Geo] Cropping a interpolated map by a polygon
Cristian Osan
felsmechanik at email.de
Tue Jul 30 09:28:56 CEST 2013
Hello everybody,
i need to create a map showing the temperature distribution in a concrete dam, as measured by sensors installed. I can do a lot of calculation on the data, working just fine. Now i want to show a picture with the interpolated temperature distribution. Up to now, i did follwing:
library(raster) # For raster
library(tgp) # For interp.
library(fields) # For image.plot
library(akima) # For akima
library(sp) # For spatial representation
#================================================================
##definition of dam contour polygon
x<-c(35.0,46.5,51.5,99.0,152.7,167.2,147.5,147.5,35.0)
y<-c(9.25, 6.70, 6.70, 7.90, 1.50, 5.80, 32.20, 43.10, 43.10)
damcontour<-cbind(x,y)
dc<-owin(poly=damcontour)
#================================================================
# # Definition of irregularly spaced -z- values
#================================================================
#read from table exported from excel
#irregularily spaced temperature sensors
temp15july<- read.delim("temp15july0845.txt", dec=",")
xx<-temp15july$x
yy<-temp15july$y
zz<-temp15july$t1507
#================================================================
# Definition of the bounding box and related variables - taken from other maillist answer
xmin<-0
xmax<-200
ymin<-0
ymax<-200
dxy <- 0.5
nx <- (xmax-xmin)/dxy
ny <- (ymax-ymin)/dxy
xbox <- c(xmin, xmax, xmax, xmin)
ybox <- c(ymin, ymin, ymax, ymax)
#================================================================
coo<-cbind(xx,yy)
prds<-predict.surface(Tps(coo,zz)) #produces the expected results!
hmr<-raster(prds) ##creates raster
image.plot(hmr, nlevel=35)
contour(prds, nlevels=35, add=T)
#================================================================
##overlay dam contour
plot(dc,add=T)
# now crop the image on the dam contour!!
I found a method in a mailing list, but i need to transform the coordinates in LatLon or so, it appears to complicated for me. I use a local rectangular coordinate system. Data file is available on request.
Any thougths?
Thanks a lot in advance
--
Dipl. Ing. Cristian Osan VDI
Geotechniker
Stuttgart/Germany
More information about the R-sig-Geo
mailing list