[STATSGRASS] Re: [GRASS5] r.moran

[Roger Bivand]

On Thu, 23 Mar 2006, Wolf Bergenheim wrote:

> On Wed, 22 Mar 2006, Roger Bivand wrote:
> 
> > On Wed, 22 Mar 2006, Wolf Bergenheim wrote:
> >
> >> On Wed, 22 Mar 2006, Roger Bivand wrote:
> >>
> > Fine. Can I suggest that you make a personal directory for the additional
> > contributed packages, then you don't need to change who you are too often?
> > In my .bashrc I have:
> >
> > R_LIBS="$HOME/lib/R" ; export R_LIBS
> 
> good idea. (: done.
> 
> > which R picks up, and uses in addition to the system-wide library. Did you
> > install R from source, or otherwise? Source is nice, because then you know
> 
> Just plain debian binary packages.
> 
> > that it figures out your machine and gets dependencies right straight
> > away, but the debian binary builders are very good too.
> >
> > Under Linux (and other platforms), the easiest way to install contributed
> > packages is from the R prompt, as you've seen. Some packages depend on
> > other R packages, but some have external dependencies, which also need
> > meeting.
> 
> This is absolutely the coolest plugin system I have ever seen! It simply 
> rocks!
> 
> >
> > You'll need the spdep package and all its dependencies, so:
> >
> > install.packages("spdep", depencencies=TRUE)
> 
> for future reference there is a typo in that. it should naturally be
> install.packages("spdep", dependencies=TRUE)
>                                 ^
>                                 +-- notice this 'd'
> 
> > does that, and the packages end up in your own library.
> >
> > The easiest way for now to import the shapefile is to do v.out.ogr from
> > GRASS, and
> 
> Ok.
> 
> >
> > library(maptools)
> > my_shape <- readShapePoly("my_shape.shp")
> > plot(my_shape)
> > names(my_shape)
> >
> > maybe
> >
> > summary(my_shape)
> 
> 
> >
> > If you are running R within GRASS (start R from the GRASS shell prompt),
> > you can give my_shape its coordinate reference system by:
> >
> > library(spgrass6)
> > meta <- gmeta6()
> > proj4string(my_shape) <- CRS(meta$proj4)
> >
> 
> I didn't do this, I'll see if I need it later..
> 
> > in case that is needed.
> >
> > Being brave:
> >
> > library(classInt)
> 
> I had to install classInp package, no biggie (:
> 
> > var_class_intervals <- classIntervals(my_shape$variable_of_interest,
> >  style="fisher")
> > cols <- findColours(var_class_intervals, pal=c("blue", "red"))
> 
> so far so good.
> 
> > plot(my_shape, cols)
> 
> my grass wanted this to be
> 
> plot(my_shape, col=cols)
> 
> >
> > and hope the plot shows the data in similar places to the PNGs.
> >
> 
> It does!
> 
> > Shall we try to get that far first - I'm not too sure how many polygons
> > there are, but there seem to be a fair number?
> 
> reading them took a while, but other then that R seems to be able to 
> handle this just fine.

Yes, each polygon - and here there are quite a few - gets checked fairly 
thoroughly, though not as thoroughly or topologically as v.clean, gets a 
centroid, an area, and a plot order to try to stop included small polygons 
being overpainted by large including ones. It could be faster, but once 
it's done, the object is there, and can be save()ed out as a 
SpatialPolygonsDataFrame object.

> 
> >
> > The next bits are simple:
> >
> > my_nb <- poly2nb(my_shape)
> 
> R says: Error: couldn't find function "poly2nb"
> 
> perhaps something is missing?
> 
> Ahh... I found it. help.search("poly2nb") talks about spdep.
> so I had to:
> 
> library(spdep)

No, I didn't do that on purpose just to see how you did! Very typically 
users and developers work in an edit/run/edit/run/... mode, I tend to use 
Vim to edit and source() to run (or copy&paste), very many power users use 
Emacs Speaks Statistics which has a run buffer, and the same facilities 
are available in RWinEdt and other editors. So a lot of things get 
forgotten first time through ...

> 
> This is simply just great. I really wish I'd have had more time to study 
> R. I guess now is as good a time as any.

Being driven by pressing need is always good motivation, IMHO!

> 
> >
> > maybe setting queen=FALSE to only choose neighbours sharing more than
> > just
> > a corner point.
> 
> Hmm it complains that my_shape isn't a polygon list... ):
> What is the correct spell?

Another slip on my part - I should have cast to a SpatialPolygons object:

my_nb <- poly2nb(as(my_shape, "SpatialPolygons"))

which will also run for a fair while.

> 
> >
> > print(my_nb)
> > plot(my_nb, coordinates(my_shape))
> >
> > display the neighbours.
> >
> > From there it is just
> >
> > moran.test(my_shape$variable_of_interest, listw=nb2listw(my_nb, style="W",
> > zero.policy=TRUE), zero.policy=TRUE)
> >
> > roughly, depending on whether you have cells that have zero neighbours.
> > But we can get back to this once the shapefile is on board.
> 
> OK. Now I have the shapefile locked and loaded, now to just pull the 
> moran numbers out of it (: and I will be solid.
> 

moran.test() returns an object of class "htest", which has a print method. 
So what you see is the output of that. How do you need the output? The

res <- moran.test(my_shape$variable_of_interest, listw=nb2listw(my_nb, 
  style="W", zero.policy=TRUE), zero.policy=TRUE)
res # just to have a look
sink("tmp/output_file.txt")
res
sink()

makes a text file copy. You can get it as HTML formatted too:

library(R2HTML)
dir.create("tmp/morantest")
HTMLStart("tmp/morantest")
HTML(res)
HTMLStop()

but you'll need to edit it a bit afterwards (or copy & paste from a 
browser to say oowriter). There is pretty good LaTeX support too, as 
you'd expect with statisticians.

By the way, if names(my_shape) says for example:

"children", "totalPop"

and you want the percentage,

my_shape$percent_children <- (100*my_shape$children)/my_shape$totalPop

does the stuff.

Or if you want the local Moran's I in the same object, do str() on the
output object of localmoran() (a matrix, nb. different access syntax than
data frames) to work our which column you want - probably Z.Ii, put it in
my_shape, and do:

writePolyShape(my_shape, "new_my_shape")

to write out to a new shapefile. At this point simple working out the 
appropriate row ID and writing out using functions in the RSQLite package 
might be the most natural for you? 

One of the ideas Markus and Helena mention in their book is writing small 
shell scripts from GRASS passing some values to R through environmental 
variables to be used in an R script with the R interface hidden from the 
user, because R can be run with a script from the command line. 

So "get the local moran values for my polygons" could be run from the
GRASS prompt if enough of the options and other user decisions can be
handled sensibly. R can write to DBF too, and most external databases, but
the very range of choices makes writing a robust script hard. In the 
polygon contiguity neighbour definition, we can actually pass out a 
neighbour structure directly from GRASS, and R can read the GRASS 
attribute databases directly, so it could be made quite "tight".

Just thinking aloud ...

Roger


> The ever grateful
> Wolf
> 
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no