MySQL and Grass

Hi,

I have a bunch (500k) of line segments(related to a 'polygon' table) and
around 70k points, both in a MySQL database. I need to do a
point-in-polygon query. i.e. I want to know which polygons contain each
of the 70k points. I decided to use MySQL because it is fast and
Postgres would not handle the size of polygons I have. (Among some other
apparently built-in limits of Postgres)

Anyway, I have been trying to do this query and it is eating my lunch.
This is what I tried...

I have a polygon table that contains the name of the entity that has a
one to many relationship with the segment table.

So trying the "vertical line to infinity and count the intersections
approach" I figured I could query for the segments that bounded the x
coord of the point where the segment was "above" the point. This gives
in SQL

"select count(*) from segment where segment.x <= x and segment.x >= x
where segment.y >=y group by segment.polygon_name"

I do get back the count of the "crossings" but I am not sure it is
correct. My polygons are supposed to be distinct and I get multiple
polygons reporting to contain the point.

Can anyone offer advice? Any good books on these types of algorithms or
code?

I would love a lot of GIS related tools as Perl modules that connect to
MySQL, shapelib, and GRASS. But if I can't get this thing going I think
I will save my pennies for ArcInfo.

Thanks in advance.
Robb

On Wed, 10 Nov 1999, Robb Hill wrote:

Can anyone offer advice? Any good books on these types of algorithms or
code?

Robb,

  That's an excellent question! I, too, would like to have such a reference
available. The last I heard -- and this goes back a dozen years when I first
worked with ARC/Info -- is that the algorithms are rather closely held by
those (like Jack Dagerman) who developed them. After all, the
point-in-polygon and similar problems are why ESRI is able to charge the
extremely high prices they do. It's also the reason why there are more
raster-based analytical GISes than vector-based analytical offerings.

  Please let me know what you find out. Have you posted to GIS-L?

Rich

Dr. Richard B. Shepard, President

                       Applied Ecosystem Services, Inc. (TM)
              Making environmentally-responsible mining happen. (SM)
                       --------------------------------
            2404 SW 22nd Street | Troutdale, OR 97060-1247 | U.S.A.
+ 1 503-667-4517 (voice) | + 1 503-667-8863 (fax) | rshepard@appl-ecosys.com

Robb Hill wrote:

Hi,

I have a bunch (500k) of line segments(related to a 'polygon' table) and
around 70k points, both in a MySQL database. I need to do a
point-in-polygon query. i.e. I want to know which polygons contain each
of the 70k points.

---------snip--------
Robb,
I have been looking at the GRASS Programmer's Manual, and the Vector
Library does include routines which might be helpful for this task, but
they then need to be called within appropriate code. Some GRASS programmers
out there may be able to offer more useful advice as to whether what you
need has already been done or can in fact be done.
The routines I refer to are:
dig_point_in_area (Map,x,y,pa) is point in area?
dig_point_to_area (Map,x,y) find which area point is in

This doesn't actually solve your problem, but it may help you get there :slight_smile:

Cheers,

  --
Terry Duell, Senior Mobility Engineer
Army Engineering Agency
Maribyrnong, Victoria, Australia
ph:61-3-93195837 fax:61-3-93195830

Hi all,

I have a set of spot height data of a site and its DEM.
I wish to perform cut & fill calculation based on a few proposed final level
using r.volume.
Is there anyone who can help me to locate some tutorial for the above
subject?

Thank you
Khin Fah

On Wed, 10 Nov 1999, Robb Hill wrote:

Can anyone offer advice? Any good books on these types of algorithms or
code?

Robb,

That's an excellent question! I, too, would like to have such a reference
available. The last I heard -- and this goes back a dozen years when I first
worked with ARC/Info -- is that the algorithms are rather closely held by
those (like Jack Dagerman) who developed them. After all, the
point-in-polygon and similar problems are why ESRI is able to charge the
extremely high prices they do. It's also the reason why there are more
raster-based analytical GISes than vector-based analytical offerings.

Please let me know what you find out. Have you posted to GIS-L?

Rich

This excercise seems to have confounded numerous folks in numerous areas of
endeavour it seems...I've appended two possible solutions (headers from
code) that I came across while searching for an answer some time ago. The
first is a simple minimalist solution, the second is a much more
sophisticated solution (can handle an array of polygons for the p_in_poly
test)...but will need more work to customise. If anybody else would like
the code (which can be freely and publicly distributed) let me know and
I'll forward same.

Regards, Roderick

/***************************************************************************
* *
* INPOLY.C *
* *
* Copyright (c) 1995-1996 Galacticomm, Inc. Freeware source code. *
* *
* Please feel free to use this source code for any purpose, commercial *
* or otherwise, as long as you don't restrict anyone else's use of *
* this source code. Please give credit where credit is due. *
* *
* Point-in-polygon algorithm, created especially for World-Wide Web *
* servers to process image maps with mouse-clickable regions. *
* *
* Home for this file: http://www.gcomm.com/develop/inpoly.c *
* *
* 6/19/95 - Bob Stein & Craig Yap *
* stein@gcomm.com *
* craig@cse.fau.edu *
* *
***************************************************************************/

/*************************************************************

Copyright (c) 1995 Raimund Seidel
All rights reserved.

This program is distributed WITHOUT ANY WARRANTY;
without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE.

**************************************************************

  A RANDOMIZED PLANAR POINT LOCATION STRUCTURE

  by R. Seidel

  last edited: 90 Feb 10

This file contains several routines for dealing with a general version
of the so-called planar point location problem, which can be described
as follows: one is to preprocess a set of non-vertical straight line
segments into a data structure so that for any given query point p one
can determine quickly which segment is the first one directly below p.
It is assumed that no two segments have interior points in common.

The method employed is a randomized one that I have so far not gotten
around to write down formally in a paper. It has the property that
a search requires fewer than 4*ln(n) comparisons in expectation. It has
O(nlogn) expected preprocessing time and needs O(n) expected space.
All expectations are over "coin flips" in the algorithm and not over
assumed input distributions.

The Test Program allows for two usages: either the input is a
set of segments as described above and for query points the next
segment below are returned, or the input is the sequence of corners
counter clockwise around a simple polygon and for query points it
is determined whether they lie inside this polygon or not.

Just compile this file and run it. It should be fairly self-explanatory.

Some warnings: There is only minimal checking whether the input
  is correct (whether there are improper intersections).
  Although the program should not crash, the answers will not
  be meaningful.
  There is also no overflow checking. Currently everything
  is done in terms of integer coordinates. Keeping their
  magnitudes smaller than 2^15 will definitely avoid overflow.
  Using a coordinate larger than `infinity' can cause a crash.

Remarks:
   The coordinate type can be changed from integer to something
   else quite easily by changing the first 5 definitions. A change
   to floating point will of course open the usual Pandora's box
   of precision problems. However, they are refined to the two
   macros point_above_segment and point_below_segment.

   These two macros along with segment_above_segment are also
   the only things that would need to be adjusted if one wanted to
   use this algorithm for segments that are not straight, but
   still x-monotone.

For any remarks, complaints, etc. send email to

  seidel@lyra.berkeley.edu

_____________________________________________________________________

Roderick Brown Tele: (Int+ 61 3) 9344-9868
School of Earth Sciences Dept.: (Int+ 61 3) 9344-7675
The University of Melbourne Facsimile: (Int+ 61 3) 9344-7761
Parkville, 3052 Home Telephone: (Int+ 61 3) 9397-8848
Australia Email: r.brown@earthsci.unimelb.edu.au

______________________________________________________________________

Hi,

Thanks to everyone for the response. I think I can now get MySQL to
work. Also I have learned something more about GRASS.

Thanks Again,
Robb