[R-sig-Geo] Distmap: the way I could do it

[Agustin Lobo]

Finally, I've given up
using distmap on an object imported
from a shp file. Could not solve the casting
problem from the sp object to the owin object.

I've rasterized the vector prior to importing to R
and this is the way I calculate the distmap and export it.

refG25 <- readGDAL("GFOEOABUFActi25.tif")
a <- as.image.SpatialGridDataFrame(refG25) #here proj info is lost
ai <- as.im(a)
ai$v[ai$v==0]<-NA #0 values in the raster are the background
adist <- distmap(as.owin(ai))
adistlog <- adist
adistlog$v <- log(adistlog$v+1) #log transform

#Export works but Global Mapper claims there is no Datum
#TNT claims there is no Gereferenciation at all (not true!)
adistlogim <- as.SpatialGridDataFrame.im(adistlog)
adistlogim at proj4string <- refG25b at proj4string #Recover proj info
writeGDAL(adistlogim, "GFODISTLOG.tif", drivername = "GTiff", type = 
"Float32", mvFlag = NA, options=NULL)

#Another alternative, same result: Datum is lost
GFODISTLOG <- refG25
GFODISTLOG at data$band1 <- adistlogim at data$v
writeGDAL(GFODISTLOG, "GFODISTLOG.tif", drivername = "GTiff", type = 
"Float32", mvFlag = NA, options=NULL)

Hope this helps somebody else. Please let me know if you think of
a better way.

Agus
-- 
Dr. Agustin Lobo
Institut de Ciencies de la Terra "Jaume Almera" (CSIC)
LLuis Sole Sabaris s/n
08028 Barcelona
Spain
Tel. 34 934095410
Fax. 34 934110012
email: Agustin.Lobo at ija.csic.es
http://www.ija.csic.es/gt/obster