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
