# In which region does my (arbitrary) point lie?

1. Mar 8, 2013

### billiards

My problem is that I have a bunch of regions, in this case they are tectonic plates. Each region has a complex boundary, described by polygons, in a text file with co-ordinates.

I have a long list of points (>100,000) distributed globally and I need to know in which region each point lies.

Basically I need an algorithm, that can be easily coded up in fortran.

Any helpers?

Last edited: Mar 8, 2013
2. Mar 8, 2013

### Staff: Mentor

I don't know how quick that is, but for 100 000 points I think you can spend some time on preprocessing: Divide the map into parts which are easy to handle (like "rectangles" in latitude/longitude), and check which plates and which edges are within a specific rectangle first. This can be done with a loop over the edges, I think. It determines all rectangles with edges, the other rectangles are easy once this is done.
Afterwards, for each point, get the corresponding rectangle (trivial). If it has more than one plate, check the borders in this rectangle.

3. Mar 8, 2013

### billiards

Thanks for your reply mfb. I was considering that kind of approach but really I am looking for something that is more transportable, so the code can be used for example with other plate boundary models, or with any arbitrary data set.

I think I may have stumbled upon the solution.

http://en.wikipedia.org/wiki/Point_in_polygon

In fact the code looks like it is available here:

http://rosettacode.org/wiki/Ray-casting_algorithm

I'll have a go and report back.

Cheers

4. Mar 8, 2013

### Staff: Mentor

Well, that needs a lot of calculations, as you have to check every edge for every point.

It can, you just have to run the preprocessing again.

5. Mar 8, 2013

### billiards

Yes but that's what computers are for :-)

6. Mar 8, 2013

### Staff: Mentor

Well, if computing time is not an issue, sure....

7. Mar 9, 2013

### billiards

Just in case anyone wwas wondering. I got my code to work.

If anone ever tries to do this and stumbles upon this thread, the above code is good: I warn you, though: Think carefully about what happens about the Greenwich Meridian (lon=0/360).

ps, I haven't run on my 100 000+ points yet but will do on Monday so I'll let ya''ll know how long it takes. cos I'm sure you're dying to kknow