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

Apologies for the off-topic solution, but PostGIS, GRASS & R are pretty interoperable :slight_smile:

I understand it may be a problem changing tools at this stage, but PostGIS can do this sort of operation very easily, if I understand you correctly.

If your data is in shapefiles, shape2pgsql will generate SQL commands to load the shapefile geometries & attribute data into PostGIS tables, and there are a number of spatial relationship query commands available to overlay the species points with polygons. This sort of query I run frequently, as I have used PostGIS for managing species occurrence data & all sorts of political/physical/management polygons for several years now. Even some of our Arc users are starting to use this, as it offers some advantages over Arc.

I appreciate the time constraints you are under, but if you can get me the shapefiles, I'm pretty sure I can provide a PostGIS based solution in an hour or two at most. And PostGIS is somewhat interoperable with Arc as well, which may be more politically acceptable in your situation.

Of course, a spatial data management/query tool like PostGIS does not provide graphical output, though the PostGIS tables can be viewed & queried on screen using several tools, QGIS, JUMP, OpenJUMP, uDIG, gvSIG spring to mind as options. The data can also be extracted & displayed or analysed further with GRASS.

My preferred FOSS solution is PostGIS for data management & some analysis, R for modeling/stats, GMT for surface modeling/gridding/publication quality cartography & QGIS for simple data viewing. Sometimes GRASS for some operations, though this is increasingly uncommon, mainly due to what I find an unduly restrictive data management capability of GRASS. But there are still some things that GRASS is the best FOSS tool for :slight_smile:

Corrado <ct529@york.ac.uk> 05/03/08 8:07 AM >>>

Dear friends,

Platform: Ubuntu 7.10, Grass 6.2.3,PostgresSQL 8.2,R 2.6.2

Yes, I am in a hurry. If I do not solve it by tomorrow, then it is back to
ARCGIS I am afraid. It is not my will, of course, but the project's leader is
already not keen on switching to GRASS, so she is pressing for me to solve it
using ARCGIS. :frowning:

Best,

On Tuesday 04 March 2008 18:33:26 Roger Bivand wrote:

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

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

Best Regards
--
Corrado Topi

Global Climate Change and Biodiversity
Area 18,Department of Biology
University of York, York, YO10 5YW, UK
Phone: + 44 (0) 1904 328645, E-mail: ct529@york.ac.uk
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user