[GRASS-user] Re: Using R to check points in polygons does not work

[Roger Bivand]

Corrado <ct529 <at> york.ac.uk> writes:

> 
> Dear friends,
> 
> I attempted to solve the problem using R, as Hamish suggested.
> 
> Because there are some ring direction changes in the shape files, the 
> conversion to polylist returns "NA" in some fields. When I run inout() on 
> those list elements the function fails with error NA/NaN/Inf in foreign 
> function call (arg 4).
> 
> Is there another method to test if a point falls in any of the polygons 
> without using R? or, is there a way to solve this problem?

Are you in a hurry? You can certainly do what you need in R, in GRASS, or 
very possibly in a combination. But you need to try to break you task 
down into smaller pieces and make sure at each step that you know what 
is going on. What is your platform (OS)?

I suggest that you import your polygon shapefile into GRASS and clean the
topology, that is first establish an appropriate location, then v.in.ogr 
(or v.in.ogr location=), then v.clean (see flags and options), and check 
that you have areas. 

If the ring direction was changed, I think you were using the read.shape()
internal function in the maptools package, followed by Map2polylist - this 
is no longer supported code, use rgdal and readOGR() instead. The polylist
structure is also deprecated. Use SpatialPolygons, into which readOGR() 
puts data  without fuss. 

Indeed, if the ring direction changed, your shapefile cannot have been 
clean (this happens, even for shapefiles generated by ArcGIS, and there 
are lots of dirty shapefiles in the wild, dirty meaning not conforming 
to the ESRI specification). Then the point-in-polygon function inout() 
from the splancs package should fail, shouldn't it? By the way, inout() 
does not work as such with polylist objects.

Access to the GRASS book (Neteler & Mitasova 2008) will help with GRASS 6 
vector operations. Then read in the point locations (are there 58000
observations on 100+ variables, or three columns: (x, y, species)? Are 
they also in a shapefile?

If you are familiar with databases (or their view of the world), you can
continue within GRASS with the db.* commands as suggested. If your output 
tables will be happiest here, that's fine (Neteler & Mitasova 2008 may 
help).

If you are going to use R with GRASS, review the GFOSS 2006 course 
materials:

http://www.foss4g2006.org/contributionDisplay.py?
 contribId=46&sessionId=59&confId=1

(sorry, gmane didn't like the long line)

(slides and scripts towards the foot of the page).

R will also give you statistical graphics, so that you can visualise your
results if need be. Tallies for conservation areas can also be moved back 
to GRASS if need be.

Starting R from the GRASS prompt, say (roughly):

library(spgrass6)
polys <- readVECT6("polys") # correct GRASS vector name
pts <- readVECT6("pts") # ditto
res0 <- overlay(pts, polys)

# yes, it's a one-liner!

This returns a vector as long as the number of points, giving either the 
number of the polygon, or NA if the point does not fall into any polygon.

Assuming that pts$species are the species, some variant of:

tapply(is.na(res0), pts$species, table)

will give you: "I need to know for each species, how many point fall 
inside conservation areas and how many points fall outside conservation 
areas".

tapply(pts$species, res0, table) or some such should tabulate the species
observations by conservation area (but you may need to drop the NAs and the
corresponding rows of pts).

If there are 100 columns of presence/absence, variants of the *apply() 
family should deliver the goods, very possibly apply() across the species
columns of pts.

The procedure for buffers could be to use v.buffer in GRASS to create the
buffers, then do the overlaying and tallying in R.

Learning enough of GRASS, or GRASS and R to do this will take a little
time, but undoubtedly less than learning C or Python, with certainly 
fewer bugs in the code you use, and with full journals (history from 
the shell in GRASS and session history() to save in R -
savehistory("now.Rhistory") is a remarkably helpful tool!

Hope this helps,

Roger


> 
> Best Regards