0

How can I simplify this for loop with a list comprehension and get more speed?

import numpy as np
import shapely

#buffer_original_points --> shapely polygons in a list
#raster --> shapely polygons in a list


raster_array = np.zeros((len(raster)))

for i in range(0, len(buffer_original_points)):
    for j in range(0, len(raster)):
        if buffer_original_points[i].intersects(raster[j]) or raster[j].contains(buffer_original_points[i]):
            raster_array[j] += 1

raster_array = raster_array.reshape(rows, columns)

I tried this, but the raster_array result was wrong and is filled with ones:

raster_array=[raster_array[j]+1 for i in range(0, len(buffer_original_points)) for j in range(0, len(raster)) if (buffer_original_points[i].intersects(raster[j]) or raster[j].contains(buffer_original_points[i]))]
Georgy
  • 12,464
  • 7
  • 65
  • 73
brez_123
  • 63
  • 1
  • 7
  • This comprehension just returns values and doesn't change the values in `raster_array`. Since you wan't to change a variable inside the loop, a loop seems to be the right approach. – Aguy Aug 14 '20 at 12:15
  • @Georgy Thanks for your support. Could you please give me some advise to implement Mike T 's solution in my case? I am a little bit overwhelmed with the topic Rtree and do not understand how to use it in my case. My aim is to build a matrix where the cell value depends on the number of points which lay in a raster cell. – brez_123 Aug 14 '20 at 21:23
  • I suggest trying the [second answer](https://stackoverflow.com/a/43105613/7851470) instead, so you won't have to install other libraries. Initialize the `STRtree` by the list of raster polygons, then iterate over `buffer_original_points` and use the `query` method to get those raster polygons that intersect with the polygon on the current iteration. To keep the track of counts you could also use `collections.Counter` with `raster`s as keys, which you would convert back to an array of counts in the end. – Georgy Aug 14 '20 at 22:09
  • @Georgy please correct me if I am wrong. The STRtree gets initialized with `s=STRtree(raster)`. After that I queried the results with a list-comprehension and stored the results in a list `result=[s.query(buffer_original_points[i] for i in range(0,len(buffer_original_points)`. Then I tested to see if any result gets true with `[raster[i] in result for i in range(0,len(raster))` but the whole list is **false**. Sorry for my misunderstanding but I am a novice. – brez_123 Aug 15 '20 at 16:14
  • The first list comprehension doesn't look right, there are missing parentheses. The second one won't work as expected because `result` is (most probably) a list of lists of polygons, but you check for the presence of `raster[i]` in it as it was a flat list of polygons. I sketched out a notebook with an example of how you would do it. Take a look: https://nbviewer.jupyter.org/gist/LostFan123/f6c46d304681b0262823637dd5e3ea5d – Georgy Aug 15 '20 at 17:28
  • 1
    Big thanks @Georgy. Your notebook helped me out enormously! From two minutes runtime to <1 sec. Thanks! – brez_123 Aug 16 '20 at 09:47

1 Answers1

0

You can probably try the following:

raster_array = [sum([buffer_original_points[i].intersects(raster[j]) or raster[j].contains(buffer_original_points[i])) for i in range(0, len(buffer_original_points))] for j in range(0, len(raster))]

But I'm not sure this is not a good coding practice due to the reduced clarity.

Aguy
  • 7,851
  • 5
  • 31
  • 58
  • Thanks for your help, but your code snippet returns me a syntax error and I can't figure out where the problem occurs. – brez_123 Aug 14 '20 at 20:53