In which region does my (arbitrary) point lie?

[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?

 

[mfb]

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.

 

[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

 

[mfb]

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

so the code can be used for example with other plate boundary models, or with any arbitrary data set.

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

 

[billiards]

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

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

 

[mfb]

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

 

[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