Point inside the polygon

Question

I'm working in Python with a .las file which I open with laspy and a .shp file which I open with the shapefile library.

I have multiple polygons in the shapefile and the task is for each polygon to take LAS points that are inside this polygon and calculate some statistics.

For now, I manually check every point from the .las file whether it is inside the polygon or not. Such a procedure is too much time-consuming. Is there a better way to deal with point extraction?

Answer 1

Leveraging one or more of the above suggested tools may well be the easiest and possibly most economical way to go. That's your call; however, if you are determined to "home-grow" the solution, the key to this problem is using an appropriate data structure to store your LAS points. As with so many spatial query problems like yours, an excellent data structure in general is a tree. Such structures spatially order or index your data, so that you can jump directly to a subset of data that makes sense for the query at hand. As I understand your problem, you want to calculate stats for the set of LAS points within each of a set of polygons. I agree that the brute force method of comparing all LAS points against each of the polygons is too time-consuming. True, brute force is often good enough for relatively small datasets and one-off analyses, but humongous data sets, as LAS data sets tend to be, make brute force an impractical solution. In general in a problem like yours, it makes no sense to compare a polygon with a point that does not fall within the bounding box of your polygon; that point has no chance of being within the footprint of the polygon.

How, then, does one determine which points fall within a polygon without going though all the points? The answer brings us back to the data structure, specifically a tree. The type of query in your case is generally called an "orthogonal range query" in computationl geometry. Typically, the data structure supporting such a query is a tree, and I specifically suggest a 2-dimensional kd-tree or a range tree. (There are plenty of other tree structures that could be used or adapted, each with its own advantages and disadvantages.) That said, one still has to build the tree and store the LAS data in it to make the data efficiently searchable, and this is the lunch that you have to pay for. It requires an algorithm, which itself requires time to develop and implement and time to run before any actual analysis is performed. Once the data are stored, though, the per polygon queries can be performed on just the sets of points that fall in the spatial ranges of each polygon. Compared to the time required initially to store the data in the tree, the time saved using a spatially filtered query vs. brute force query can be huge with the biggest advantages coming from processing large data sets. There are myriad computational geometry resources out there, but, as the one I most oftern start with sits literally a foot away from me in my small reference book stack, I will offer that one: de Berg et al. 2008. Computational Geometry (3rd ed.). Springer-Verlag, Berlin; 398 pp.