I am trying to calculate an intersectLCA, which requires the (ecoinvent, world) intersection in bw2regional (bwr), but the geo_transform matrix is empty, because the intersection is empty.
I use bw2r version 0.5.2, compatible with brightway2, but to calculate intersections:
- the pandarus remote server is not available
- the pandarus library says the two geocollections are overlapping
Newer version of bwr (0.7.DEV2) also provides geopandas as a way to calculate the intersections. So I replicated the function to calculate intersection with geopandas, which provides intersections data. The intersection is then writen. However, when loading it during the intersectLCA computation, the intersection array is only full of 4294967295 values, which are kind of error values. Therefore, the array is empty, and the geotransform matrix is null.
Why do I get an array filled whith when loading the intersection, altough the intersection is not empty when writing it ?
Here is the intersection creating code:
df1 = gpd.read_file(bwr.geocollections['world']["filepath"])
df2 = gpd.read_file(bwr.geocollections['ecoinvent']["filepath"])
id1 = bwr.geocollections['world']["field"]+'_1'
id2 = bwr.geocollections['ecoinvent']["field"]+'_2'
# creating the intersection
intersection = gpd.overlay(df1, df2, how='intersection', keep_geom_type=False)
# cleaning where no iso code
intersection = intersection[intersection['isotwolettercode_1'].notna()]
# filtering for weird geometry
intersection = intersection[(intersection.geom_type=='Polygon')|(intersection.geom_type=='MultiPolygon')]
# getting area
areas = intersection.to_crs("esri:54009").area # World Mollweidge in square meters
# writing intersection in bwr format
data = []
for index, feature in intersection.iterrows():
data.append((('world', feature[id1]), ('ecoinvent', feature[id2]), areas[index]))
obj = bwr.Intersection(('world', 'ecoinvent'))
obj.write(data)
obj.create_reversed_intersection()
when displaying data, intersection has values:
data[:10]
Returns:
[(('world', 'NF'), ('ecoinvent', 'UN-AUSTRALIANZ'), 41293322.44087714),
(('world', 'NZ'), ('ecoinvent', 'UN-AUSTRALIANZ'), 268812699946.50467),
(('world', 'TK'), ('ecoinvent', 'UN-AUSTRALIANZ'), 5481892.042473819),
(('world', 'AU'), ('ecoinvent', 'UN-AUSTRALIANZ'), 7724172328042.095),
(('world', 'NF'), ('ecoinvent', 'UN-OCEANIA'), 41293322.44087714),
(('world', 'NU'), ('ecoinvent', 'UN-OCEANIA'), 221754487.87432534),
(('world', 'VU'), ('ecoinvent', 'UN-OCEANIA'), 12366728898.030743),
(('world', 'KI'), ('ecoinvent', 'UN-OCEANIA'), 985794360.8650136),
(('world', 'PF'), ('ecoinvent', 'UN-OCEANIA'), 3356087824.517741),
(('world', 'NZ'), ('ecoinvent', 'UN-OCEANIA'), 268812699946.50467)]
But when loading intersections values (used in matrix computation after), I only get the error values 4294967295 in row and columns, which leads then to an empty matrix when passing the array to the matrix builder.
bc.utils.load_arrays([bwr.Intersection(('ecoinvent', 'world')).filepath_processed()])
Returns:
array([(569, 1869, 4294967295, 4294967295, 0, 1.2212809e+04, 1.2212809e+04, nan, nan, nan, nan, False),
(570, 1869, 4294967295, 4294967295, 0, 1.2212809e+04, 1.2212809e+04, nan, nan, nan, nan, False),
(571, 1869, 4294967295, 4294967295, 0, 1.2212809e+04, 1.2212809e+04, nan, nan, nan, nan, False),
...,
(587, 1777, 4294967295, 4294967295, 0, 1.6897256e+13, 1.6897256e+13, nan, nan, nan, nan, False),
(672, 1777, 4294967295, 4294967295, 0, 1.6897256e+13, 1.6897256e+13, nan, nan, nan, nan, False),
(673, 1777, 4294967295, 4294967295, 0, 1.6897256e+13, 1.6897256e+13, nan, nan, nan, nan, False)],
dtype=[('geo_inv', '<u4'), ('geo_ia', '<u4'), ('row', '<u4'), ('col', '<u4'), ('uncertainty_type', 'u1'), ('amount', '<f4'), ('loc', '<f4'), ('scale', '<f4'), ('shape', '<f4'), ('minimum', '<f4'), ('maximum', '<f4'), ('negative', '?')])
From what I understand, row and col should not be 4294967295, otherwise the matrix build from this array is empty
Here is the code of the LCA I am performing. It involves the Ecoinvent and Agribalyse Database. Prior to this code, the databases activity have been labeled, and the row locations have been taken care of with the rower datapackage.
# Naming databases
bio3 = bd.Database('biosphere3')
bioagb31 = bd.Database('agribalyse3.1 biosphere flows')
eidb = bw.Database('ecoinvent 3.8_cutoff_ecoSpold02')
agbalyse = bw.Database('agribalyse3.1')
# New database
my_db = bd.Database("Regionalised Bread")
# Inventory (a bread)
wheat = my_db.new_activity(code="Wheat",name="Wheat",location='FR',unit='kg')
wheat.save()
wheat.new_exchange(input=agbalyse.search('Soft wheat grain organic moisture')[0], amount=1, type="technosphere").save()
flour = my_db.new_activity(code="Flour",name="Flour",location='FR', unit='kg')
flour.save()
flour.new_exchange(input=eidb.search('market for electricity medium voltage', filter={'location':'FR'})[0], amount=0.114, type="technosphere").save()
flour.new_exchange(input =wheat, amount =1.25, type="technosphere").save()
breadbaked = my_db.new_activity(code="Bread at plant",name="Bread at plant",location='FR', unit='kg')
breadbaked.save()
breadbaked.new_exchange(input=eidb.search('market group for tap water', filter={'location':'RER'})[0], amount=1.69, type="technosphere").save()
breadbaked.new_exchange(input =agbalyse.search('Heat, district or industrial, natural gas {FR}')[0], amount =30.55, type="technosphere").save()
breadbaked.new_exchange(input = flour, amount=0.77, type="technosphere").save()
breadsaled = my_db.new_activity(code="Bread at sale",name="Bread at sale",location='FR', unit='kg')
breadsaled.save()
breadsaled.new_exchange(input=eidb.search('transport freight light commercial vehicle', filter ={'location': 'Europe'})[0], amount=0.03, type="technosphere").save()
breadsaled.new_exchange(input =breadbaked, amount =1, type="technosphere").save()
# labelling new db
bwr.label_activity_geocollections('Regionalised Bread')
The impact assessment method has already been created, it is regionalised at the country level, and looks like this:
landsrmethod = bw.Method(('Resource Supply Risk Potential (JRC)','Land Supply Risk Potential (JRC)'))
landsrmethod.load()[:5]
[(('biosphere3', 'e489cce4-a80f-417d-9ae6-9fc14cc7dd49'),
0.0604505136445917,
'AT'),
(('biosphere3', 'e489cce4-a80f-417d-9ae6-9fc14cc7dd49'),
0.8142923129937788,
'BE'),
(('biosphere3', 'e489cce4-a80f-417d-9ae6-9fc14cc7dd49'),
2.1836511415932987,
'CZ'),
(('biosphere3', 'e489cce4-a80f-417d-9ae6-9fc14cc7dd49'),
0.0827477341881703,
'DE'),
(('biosphere3', 'e489cce4-a80f-417d-9ae6-9fc14cc7dd49'),
0.5652521448442248,
'DK')]
landsrmethod.metadata
{'unit': 'SR',
'description': 'Land Supply Risk from the JRC criticality method, regionalised at country level',
'geocollections': ['world'],
'abbreviation': 'resource-supply-risk-potential-jrcl.11b595f16c6df637798224663df145f3',
'num_cfs': 1752}
One spatial scale LCA works:
lca = bwr.OneSpatialScaleLCA({breadsaled: 1}, ('Resource Supply Risk Potential (JRC)','Land Supply Risk Potential (JRC)'))
lca.lci()
lca.lcia()
lca.score
0.18872489371142512
For intersect LCA, intersections are created with code described above. Three intersections are created (and reversed) : (world, ecoinvent), (world Row), (world, world)
bwr.intersections
Intersections dictionary with 5 object(s):
('RoW', 'world')
('ecoinvent', 'world')
('world', 'RoW')
('world', 'ecoinvent')
('world', 'world')
But the intersectLCA returns a zero:
# Two spatial scale
lca = bwr.TwoSpatialScalesLCA({breadsaled: 1}, ('Resource Supply Risk Potential (JRC)','Land Supply Risk Potential (JRC)'))
lca.lci()
lca.lcia()
lca.score
0
When looking into the results, this is due to a null geo_transform_matrix:
lca.geo_transform_matrix
<599x24 sparse matrix of type '<class 'numpy.float64'>'
with 0 stored elements in Compressed Sparse Row format>
4294967295 is the maximum unsigned 32 bit integer - this is supposed to be there when you load the intersections file, as it is overwritten when building the actual matrix.
BTW, excellent work getting through a package which is not well documented. Very impressive!
Can you provide the LCA construction code you are using, and the point where you notice the actual error?