[R-sig-Geo] Exporting spatstat 'pixel image' to Arcview 9.x as an e00 or ascii
M.Stevenson at massey.ac.nz
M.Stevenson at massey.ac.nz
Tue Jul 12 12:29:57 CEST 2011
This might be of some use. Here (among other things) we kernel smooth a
point pattern using spatstat then export the raster image as an ASCII grid
ready for loading into a GIS.
Regards,
Mark Stevenson
------------------------------------------------------------------------
library(epiR); library(maptools); library(rgdal); library(shapefiles);
library(sp); library(spatstat); library(spatialkernel);
library(RColorBrewer)
data(epi.incin); dat <- epi.incin
xlim <- range(dat$xcoord); ylim <- range(dat$ycoord)
sigma <- 1200; dimyx = c(200, 200)
# Create an observation window and two ppp objects:
dat.w <- convexhull.xy(x = dat[,1], y = dat[,2])
cas.ppp <- ppp(x = dat$xcoord[dat$status == 1], y = dat$ycoord[dat$status
== 1], window = dat.w)
ctl.ppp <- ppp(x = dat$xcoord[dat$status == 0], y = dat$ycoord[dat$status
== 0], window = dat.w)
cas.den <- density(cas.ppp, sigma = sigma, dimyx = dimyx)
ctl.den <- density(ctl.ppp, sigma = sigma, dimyx = dimyx)
# Image plot showing the spatial distribution of laryngeal cancer risk:
breaks <- seq(from = 0, to = 0.4, length = 5)
col <- brewer.pal(n = 4, name = "Blues")
par(pty = "s")
plot(xlim, ylim, type = "n", xlab = "Easting (m)", ylab = "Northing (m)",
xlim = xlim, ylim = ylim)
image(x = cas.den$xcol, y = cas.den$yrow, z = (t(cas.den$v) /
t(ctl.den$v)), zlim = c(0, 0.4), col = col, breaks = breaks, add = TRUE)
points(ctl.ppp, pch = 1, col = "grey")
points(cas.ppp, pch = 16)
plot(dat.w, add = TRUE)
metre(xl = xlim[1], yb = ylim[1], xr = xlim[1] + 500, yt = ylim[1] + 5000,
lab = breaks, cols = col, shift = 0, cex = 0.75)
# Add the appropriate risk (really, an odds) estimate to each household
location. Convert the risk estimates to a SpatialGridDataFrame:
risk <- cas.den
risk$v <- (cas.den$v) / (ctl.den$v)
risk.sgdf <- as.SpatialGridDataFrame.im(risk)
# Convert the case locations to a SpatialPoints object:
cas.sp <- SpatialPoints(coords = cbind(cas.ppp$x, cas.ppp$y))
# Do a transect and add the values to the cas.sp table:
tmp.trans <- overlay(risk.sgdf, cas.sp)
cas <- data.frame(x = cas.ppp$x, y = cas.ppp$y, risk = tmp.trans$v)
# Plot the data adding the numeric estimates of risk to check:
par(pty = "s")
plot(xlim, ylim, type = "n", xlab = "Easting (m)", ylab = "Northing (m)",
xlim = xlim, ylim = ylim)
image(x = cas.den$xcol, y = cas.den$yrow, z = (t(cas.den$v) /
t(ctl.den$v)), zlim = c(0, 0.4), col = col, breaks = breaks, add = TRUE)
text(x = cas$x, cas$y, labels = round(cas$risk, digits = 2), cex = 0.50)
plot(dat.w, add = TRUE)
metre(xl = xlim[1], yb = ylim[1], xr = xlim[1] + 500, yt = ylim[1] + 5000,
lab = breaks, cols = col, shift = 0, cex = 0.75)
# As a further check, export all of the R objects you've created into a
format suitable for display in a GIS.
# First of all, export the boundary window created earlier as a shape file
using the shapefiles package:
Id <- rep(1, times = length(dat.w$bdry[[1]]$x))
X <- dat.w$bdry[[1]]$x
Y <- dat.w$bdry[[1]]$y
shpTable <- data.frame(Id, X, Y)
attTable <- data.frame(Id = 1, Name = c("VAR1"))
epiincin <- convert.to.shapefile(shpTable = shpTable, attTable = attTable,
field = "Id", type = 5)
write.shapefile(epiincin, "epiincin", arcgis = TRUE)
# Export the raster image in ESRI's asciigrid format:
xllcorner = epiincin$shp$shp[[1]]$box[1]
yllcorner = epiincin$shp$shp[[1]]$box[2]
cellsize.h <- (epiincin$shp$shp[[1]]$box[3] -
epiincin$shp$shp[[1]]$box[1]) / dimyx[1]
cellsize.v <- (epiincin$shp$shp[[1]]$box[4] -
epiincin$shp$shp[[1]]$box[2]) / dimyx[2]
risk <- t((cas.den$v) / (ctl.den$v))
ncol <- dim(risk)[1]
risk <- (risk[,ncol:1])
grd <- GridTopology(cellcentre.offset = c(xllcorner,yllcorner), cellsize =
c(cellsize.h, cellsize.v), cells.dim = dimyx)
crds <- coordinates(grd)
risk.sgdf <- SpatialGridDataFrame(grd, data = data.frame(z =
as.vector(risk)))
spplot(risk.sgdf)
writeGDAL(risk.sgdf, "risk_sgdf.asc")
# Export the household data (including the risk estimates):
write.table(cas, "cases.csv", sep = ",", row.names = FALSE, col.names = TRUE)
# These three objects can now be loaded and queried in a GIS.
------------------------------------------------------------------------
> @Rolf
>
> You mentioned:
> xy <- as.matrix(as.data.frame(raster.xy(as.owin(Im))))
> v <- lookup.im(Im,xy[,1],xy[,2])
>
> Just for information, can you tell me in which package the functions
> lookup.im and raster.xy can be found?
>
> Thanks!
>
> 2011/7/12 Rolf Turner <r.turner at auckland.ac.nz>
>
>>
>> This is definitely a case of the blind leading the blind, but I thought
>> I'd
>> chip
>> in anyway.
>>
>> Let "Im" be your pixel image. Do:
>>
>> xy <- as.matrix(as.data.frame(**raster.xy(as.owin(Im))))
>> v <- lookup.im(Im,xy[,1],xy[,2])
>> A <- SpatialPixelsDataFrame(points=**xy,data=data.frame(v))
>> B <- as(A,"SpatialGridDataFrame")
>>
>> Then "B" is an object of class "SpatialGridDataFrame" which ***may*** be
>> what
>> you need. I can't really tell, since I don't understand
>> SpatialGridDataFrames.
>>
>> The foregoing is kludgey, almost surely sub-optimal, and quite possibly
>> wrong .....
>> But there's a chance it will be of some use to you.
>>
>> cheers,
>>
>> Rolf Turner
>>
>>
>> On 12/07/11 10:54, Honey-Marie de la Giroday wrote:
>>
>>> Dear All,
>>>
>>> I have a pixel image (im) I received as a .RData file. I need to
>>> export this pixel image from R in a file format that can be imported
>>> into Arcview 9.x (e.g., .e00 or ascii).
>>>
>>> Below is some of the code that I have been trying to use. I have
>>> reviewed the manuals for rgdal, maptools, sp, and spatstat packages
>>> but have been unable to find another approach to converting the pixel
>>> image to an Arcview raster dataset. Please note that I have also
>>> searched the R-sig-geo archives and was unable to find this covered
>>> elsewhere.
>>>
>>> "#load spatst package
>>> library(spatstat)
>>>
>>> #load file with model data
>>> load(file="C:/bestmodel2004.**RData")
>>>
>>> #Plot best model data for 2004
>>> plot(bestmodel2004)
>>>
>>> #Examine attributes of pixel image
>>> attributes(bestmodel2004)
>>> $names
>>> [1] "v" "dim" "xrange" "yrange" "xstep" "ystep" "xcol"
>>> "yrow"
>>> [9] "type" "units"
>>>
>>> $class
>>> [1] "im"
>>>
>>> #Load sp package
>>> library(sp)
>>>
>>> #Coerce bestmodel2004 (pixel image) to SpatialGridDataFrame
>>> #......how do I do this?......
>>>
>>> #Output to Acrview = write.asciigrid(x, fname, attr=1,
>>> na.value=-9999.....
>>> #x=SpatialGridDataFrame
>>> #fname=filename
>>> write.asciigrid"
>>>
>>> Thank you for your time, and please let me know if I need to provide
>>> further information.
>>>
>>> Sincerely,
>>>
>>>
>>>
>>> H.M.C de la Giroday
>>>
>>> ______________________________**_________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo at r-project.org
>>> https://stat.ethz.ch/mailman/**listinfo/r-sig-geo<https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
>>>
>>
>> ______________________________**_________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/**listinfo/r-sig-geo<https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
>>
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
More information about the R-sig-Geo
mailing list