I have an irregular grid of points in the form (X,Y,Z), where (X,Y) are coordinates on a plane (can be geographical longitude/latitude), Z is some property to interpolate between points. I use akima interpolation package in R. The data set can contain missed values, and akima package does not like it. This can be remedied by a complete.cases() directive reorganizing the data set. But there is a following issue. Certain points contain no data in the sense that the interpolated quality is absent there (NA in R). As a closest example, Z is a depth of a stratigraphic interval, for example, Quaternary deposits. In these places I need to have a "hole" in the interpolated grid, showing that this layer is absent here; meanwhile the algorithm simply interpolates between available points with data.
#small extract from data
mydf<-read.csv(text="lon, lat, Q
411,362,1300
377,395.5,1425
427,370,1800
435.5,352,
428,357,
390,423,1700")
library("akima")
bbb<-data.frame(lon=mydf$lon,lat=mydf$lat,H=mydf$Q)
ccc<-bbb[complete.cases(bbb),]
Q.int<-interp(ccc$lon,ccc$lat,ccc$H,linear=TRUE,
extrap=FALSE,nx=240,ny=240,duplicate="mean")
Then it can be visualized, for example, with image.plot() from fields package
library("fields")
image.plot(Q.int)
In this data set, points 4,5 are lacking. This can be either 1) lack of data on these points, or 2) indication, that the deposits are absent here. I can note it in data set explicitly, for example, with NA symbol. But I need the interpolation, which distinguishes these two cases. My solution was to interpolate "as is", and then to use a trick: declare that all interpolated values of a property Z on a grid <30 meters are effectively NA's, and then plot it:
Q.int$z.cut<-ifelse(Q.int$z<30,NA,Q.int$z)
This could reflect the geological situation, since layers with decreasing thickness indeed "fade out", and can be stripped on a map, but would it be possible to arrange this problem in a more elegant solution?
