[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