From 565244559fa985849fcb08ab7ead26d8204731e0 Mon Sep 17 00:00:00 2001 From: martavp Date: Thu, 12 Oct 2023 20:43:32 +0200 Subject: [PATCH] fix bug when calculating intersection of areas in gis module 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" --- atlite/gis.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/atlite/gis.py b/atlite/gis.py index a92310c5..22632a2a 100644 --- a/atlite/gis.py +++ b/atlite/gis.py @@ -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