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

[Brent Wood]

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

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 :-)






>>> Corrado <ct529 at 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. :(

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 at 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 at york.ac.uk
_______________________________________________
grass-user mailing list
grass-user at lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user