[GRASS-user] Checking if a point (X, Y) is inside the, polygons of a shape file using command line. Solution using R

Corrado:

This is a good task for the R programming environment, which has strong geospatial tools and
integrates with GRASS:

Check out the following example, which solves the point-in-polygon problem for shape files:

http://nceas.ucsb.edu/scicomp/GISSeminar/UseCases/PointInPolygonAnalysis/point_in_poly.html

Some other good examples on this site as well...
Regards,
Rick Reeves

Rick Reeves wrote:

This is a good task for the R programming environment, which has
strong geospatial tools and integrates with GRASS:

FYI, in GRASS's vector library (available outside of C with e.g.
SWIG-Python) there a number of x-in-polygon functions found in
lib/vector/Vlib/poly.c. v.buffer uses such a test to check that its
centroid is inside the new area, you could trace that back if you
wanted to see the method in use.

from GRASS modules there is v.select and 'v.distance dmax=0' (?) and
probably other ways too.

many roads,
Hamish

      ____________________________________________________________________________________
Be a better friend, newshound, and
know-it-all with Yahoo! Mobile. Try it now. http://mobile.yahoo.com/;_ylt=Ahu06i62sR8HDtDypao8Wcj9tAcJ

Dear John, Rick, Hamish,

thanks for your kind answers.

I explain the problem:

I have a file with some hundreds of species and 58,000 observations. For each
species, I have the coordinates (X,Y) of the sites where they were observed.

I also have a shape file, with polygons representing conservation areas, like
natural reserves.

I need to know for each species, how many point fall inside conservation areas
and how many points fall outside conservation areas. Also, I need to run the
same test on the conservation areas with a buffer of 100 m and with a buffer
of 500 m.

This is what I thought: I test if each point (that is site) is inside the
polygon (that is conservation area), and I have an answer like yes or no (or
0 or 1). I write a routine that does test all the points for each species.

I add the buffer (I do not know how to do it, but I think it is possible). I
run the test again.

What do you think of this approach?

If I use v.select, would that return me 0 or 1 (or yes or no), when I try to
overlap one point with the vector map?

Or do you think it would be better to use R?

Or would it be better to prepare a vector map for each species with all the
sites, and overlap species by species instead of point by point?

I apologise for my questions, but I am trying to use GRASS for the first time.
I would like to switch from ArcGIS and convince also the other people in the
research group (there is at least two of us campaigning now).

Thanks!

On Tuesday 04 March 2008 00:38:26 Hamish wrote:

Rick Reeves wrote:
> This is a good task for the R programming environment, which has
> strong geospatial tools and integrates with GRASS:

FYI, in GRASS's vector library (available outside of C with e.g.
SWIG-Python) there a number of x-in-polygon functions found in
lib/vector/Vlib/poly.c. v.buffer uses such a test to check that its
centroid is inside the new area, you could trace that back if you
wanted to see the method in use.

from GRASS modules there is v.select and 'v.distance dmax=0' (?) and
probably other ways too.

many roads,
Hamish

___________________________________________________________________________
_________ Be a better friend, newshound, and
know-it-all with Yahoo! Mobile. Try it now.
http://mobile.yahoo.com/;_ylt=Ahu06i62sR8HDtDypao8Wcj9tAcJ

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

Hamish:

Good point. For a large dataset such as Corrado's, it would be better to use GRASS,
once Corrado surmounts the GRASS learning curve. R is a good alternative choice
in that it is simpler to use.

Regards, RR

Hamish wrote:

Rick Reeves wrote:
  

This is a good task for the R programming environment, which has
strong geospatial tools and integrates with GRASS:
    
FYI, in GRASS's vector library (available outside of C with e.g.
SWIG-Python) there a number of x-in-polygon functions found in
lib/vector/Vlib/poly.c. v.buffer uses such a test to check that its
centroid is inside the new area, you could trace that back if you
wanted to see the method in use.

from GRASS modules there is v.select and 'v.distance dmax=0' (?) and
probably other ways too.

many roads,
Hamish

      ____________________________________________________________________________________
Be a better friend, newshound, and know-it-all with Yahoo! Mobile. Try it now. http://mobile.yahoo.com/;_ylt=Ahu06i62sR8HDtDypao8Wcj9tAcJ

I think your approach with two runs is adequate, although I am not sure which tool you would use to extend your polygon as I have not looked into that issue recently. I have similar questions that I address. In order to get names of the polygons (conservation areas for you), I go beyond v.select to populate a field with the names of the areas being queried. I then dump to a text file that I can open in a spreadsheet application for summarizing. My example assumes you have a field called 'name' in your input polygon file:

v.select ainput=in_points binput=conservation_polygons out=points_in_polygons
v.db.addcol points_in_polygons col="polygon_name VARCHAR(40)"
v.distance from= points_in_polygons to=conservation_polygons dmax=0 upload=to_attr to_column=name col=polygon_name
v.db.select points_in_polygons > points_in_polygons.txt

Hope that helps. Maybe others will have similar or improved solutions, but this has worked for me.

John

On Mar 4, 2008, at 3:00 AM, Corrado wrote:

Dear John, Rick, Hamish,

thanks for your kind answers.

I explain the problem:

I have a file with some hundreds of species and 58,000 observations. For each
species, I have the coordinates (X,Y) of the sites where they were observed.

I also have a shape file, with polygons representing conservation areas, like
natural reserves.

I need to know for each species, how many point fall inside conservation areas
and how many points fall outside conservation areas. Also, I need to run the
same test on the conservation areas with a buffer of 100 m and with a buffer
of 500 m.

This is what I thought: I test if each point (that is site) is inside the
polygon (that is conservation area), and I have an answer like yes or no (or
0 or 1). I write a routine that does test all the points for each species.

I add the buffer (I do not know how to do it, but I think it is possible). I
run the test again.

What do you think of this approach?

If I use v.select, would that return me 0 or 1 (or yes or no), when I try to
overlap one point with the vector map?

Or do you think it would be better to use R?

Or would it be better to prepare a vector map for each species with all the
sites, and overlap species by species instead of point by point?

I apologise for my questions, but I am trying to use GRASS for the first time.
I would like to switch from ArcGIS and convince also the other people in the
research group (there is at least two of us campaigning now).

Thanks!

Dear friends,

Where do I find these files, like poly.c? I cannot find them on my
installation. I have 6.2.3 on Ubuntu 7.10, with the grass-dev package
installed.

Do I have to use python?

I thought I could use a bash script.

On Tuesday 04 March 2008 16:45:47 Rick Reeves wrote:

Hamish:

Good point. For a large dataset such as Corrado's, it would be better to
use GRASS,
once Corrado surmounts the GRASS learning curve. R is a good alternative
choice
in that it is simpler to use.

Regards, RR

Hamish wrote:
> Rick Reeves wrote:
>> This is a good task for the R programming environment, which has
>> strong geospatial tools and integrates with GRASS:
>
> FYI, in GRASS's vector library (available outside of C with e.g.
> SWIG-Python) there a number of x-in-polygon functions found in
> lib/vector/Vlib/poly.c. v.buffer uses such a test to check that its
> centroid is inside the new area, you could trace that back if you
> wanted to see the method in use.
>
> from GRASS modules there is v.select and 'v.distance dmax=0' (?) and
> probably other ways too.
>
>
> many roads,
> Hamish
>
>
>
>
> _________________________________________________________________________
>___________ Be a better friend, newshound, and
> know-it-all with Yahoo! Mobile. Try it now.
> http://mobile.yahoo.com/;_ylt=Ahu06i62sR8HDtDypao8Wcj9tAcJ

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

On Tue, March 4, 2008 18:04, John C. Tull wrote:

I think your approach with two runs is adequate, although I am not
sure which tool you would use to extend your polygon as I have not
looked into that issue recently. I have similar questions that I
address. In order to get names of the polygons (conservation areas for
you), I go beyond v.select to populate a field with the names of the
areas being queried. I then dump to a text file that I can open in a
spreadsheet application for summarizing. My example assumes you have a
field called 'name' in your input polygon file:

v.select ainput=in_points binput=conservation_polygons
out=points_in_polygons
v.db.addcol points_in_polygons col="polygon_name VARCHAR(40)"
v.distance from= points_in_polygons to=conservation_polygons dmax=0
upload=to_attr to_column=name col=polygon_name
v.db.select points_in_polygons > points_in_polygons.txt

Just out of curiosity: do you really need the v.select ? Wouldn't this work:

v.db.addcol in_points col="polygon_name VARCHAR(40)"
v.distance from=in_points to=conservation_polygons dmax=0 /
   upload=to_attr to_column=name col=polygon_name

?

Moritz

On Mar 4, 2008, at 1:08 PM, Moritz Lennert wrote:

On Tue, March 4, 2008 18:04, John C. Tull wrote:

I think your approach with two runs is adequate, although I am not
sure which tool you would use to extend your polygon as I have not
looked into that issue recently. I have similar questions that I
address. In order to get names of the polygons (conservation areas for
you), I go beyond v.select to populate a field with the names of the
areas being queried. I then dump to a text file that I can open in a
spreadsheet application for summarizing. My example assumes you have a
field called 'name' in your input polygon file:

v.select ainput=in_points binput=conservation_polygons
out=points_in_polygons
v.db.addcol points_in_polygons col="polygon_name VARCHAR(40)"
v.distance from= points_in_polygons to=conservation_polygons dmax=0
upload=to_attr to_column=name col=polygon_name
v.db.select points_in_polygons > points_in_polygons.txt

Just out of curiosity: do you really need the v.select ? Wouldn't this work:

v.db.addcol in_points col="polygon_name VARCHAR(40)"
v.distance from=in_points to=conservation_polygons dmax=0 /
upload=to_attr to_column=name col=polygon_name

?

Moritz

I suppose you might be able to do that. I used the v.select to first subsample a larger dataset to create a file of points that fall within the polygons. Your approach would append polygon_name to the entire points file, which may be perfectly suitable for some peoples needs. Unless I am missing something, a definite possibility.

John

Corrado wrote:

Where do I find these files, like poly.c? I cannot find them on my
installation. I have 6.2.3 on Ubuntu 7.10, with the grass-dev package
installed.

Do I have to use python?

I thought I could use a bash script.

I was meaning fundamentally those functions are present in the C source
code if you want to have a peek under the bonnet:
  http://trac.osgeo.org/grass/browser/grass/trunk/lib/vector/Vlib/

If that's a bit dense you can check out what the Programmer's module:
  http://download.osgeo.org/grass/grass6_progman/

But typically you should be able to get at whatever functionality you
need from the modules (command line / GUI interfaces). There isn't much
you can't do with the modules, it may not be obvious which module
contains the functionality, but if you keep searching+asking eventually
a solution comes to light.

If you must have access to a function not covered by a module or e.g.
want to write a high performance custom model into a new module, you
can get access to those C library functions using the SWIG Python
interface without having to mess with knowing C.

Hamish

      ____________________________________________________________________________________
Be a better friend, newshound, and
know-it-all with Yahoo! Mobile. Try it now. http://mobile.yahoo.com/;_ylt=Ahu06i62sR8HDtDypao8Wcj9tAcJ

On Tue, Mar 4, 2008 at 12:00 PM, Corrado <ct529@york.ac.uk> wrote:

Dear John, Rick, Hamish,

thanks for your kind answers.

I explain the problem:

I have a file with some hundreds of species and 58,000 observations. For each
species, I have the coordinates (X,Y) of the sites where they were observed.

I also have a shape file, with polygons representing conservation areas, like
natural reserves.

I need to know for each species, how many point fall inside conservation areas
and how many points fall outside conservation areas.

(while digging around myself for a solution, I found one and thought to
post it here as followup:)

http://grass.osgeo.org/wiki/Count_points_in_polygon

Perhaps helpful?

Markus