Skip to content

Commit

Permalink
fix bug when calculating intersection of areas in gis module
Browse files Browse the repository at this point in the history
The current intersection of areas is using tree.query(o) but this provides the index of the elements in "tree" that overlaps with "o". This gives an error because it attempts to calculate the intersection of an index "o" with a geometry "d". The intersection must be calculated between the geometry corresponding to the index "o", that is "orig[o]" and the element in the tree "d"
  • Loading branch information
martavp committed Oct 12, 2023
1 parent 506f76c commit 5652445
Showing 1 changed file with 5 additions and 5 deletions.
10 changes: 5 additions & 5 deletions atlite/gis.py
Original file line number Diff line number Diff line change
Expand Up @@ -147,13 +147,13 @@ def compute_indicatormatrix(orig, dest, orig_crs=4326, dest_crs=4326):
indicator = sp.sparse.lil_matrix((len(dest), len(orig)), dtype=float)
tree = STRtree(orig)
idx = dict((id(o), i) for i, o in enumerate(orig))

for i, d in enumerate(dest):
for o in tree.query(d):
if o.intersects(d):
j = idx[id(o)]
area = d.intersection(o).area
indicator[i, j] = area / o.area
if orig[o].intersects(d):
j = idx[id(orig[o])]
area = d.intersection(orig[o]).area
indicator[i, j] = area / orig[o].area

return indicator

Expand Down

0 comments on commit 5652445

Please sign in to comment.