brightway-lca / brightway2-regional

Regionalized LCA calculations
https://docs.brightway.dev/projects/bw-regional/en/latest/index.html
BSD 3-Clause "New" or "Revised" License
6 stars 3 forks source link

calculate_intersection with non-remote engines #14

Open mfastudillo opened 1 month ago

mfastudillo commented 1 month ago

sometimes pandarus_remote seems to be unavailable so I was trying to run the calculations locally. Specifically I was trying to replace the function create_regionalized_intersections in bw2_lcimpact

I started by the intersections, that function calls

remote.intersection("world","watersheds-eq-sw-certain")

I replaced that by

bwr.gis_tasks.calculate_intersection("world", "watersheds-eq-sw-certain",engine='geopandas')

but I get a "KeyError: 'id'" ... which is strange because the bwr.geocollections["watersheds-eq-sw-certain"]['field'] is "id", and it seems to be the right column in the gpkg. Changing the engine to pandarus gives an AttributeError

mfastudillo commented 1 month ago

... when it does the intersection it ends up with a id_1 and id_2 columns... The problem is that at least one of the fields is also in the other dataframe. Both the "world" and "whatersheds" have columns named id

A quick fix could be to rename the columns of the other dataframe...

mfastudillo commented 1 month ago

something like this may work, although it is not obvious how to write the tests

https://github.com/brightway-lca/brightway2-regional/compare/main...mfastudillo:brightway2-regional:iss14_gpd_intersection

mfastudillo commented 1 month ago

the second part of create_regionalized_intersections in bw2_lcimpact calls

bwr.pandarus_remote.intersection_as_new_geocollection('world', 'watersheds-hh','world-wathersheds-hh')

which doesn't work at the moment either... I guess to replace this by a local calculation we would need to do something like

bwr.geocollections['world-wathersheds-hh'] = {
'filepath': 'some_path',
'field':'a_field',
'is_intersection':True,
"first": bwr.geocollections['world'], # this is a guess
"second": bwr.geocollections['watersheds-hh'],
}

but I don't know what those fields should be.

After "fixing" bwr.gis_tasks.calculate_intersection I can see in bwr.intersections[('watersheds-hh','world')] an abbriviation, but I don't know where that intersection has been stored... @cmutel :)

nrogy commented 2 weeks ago

Dear @mfastudillo,

Intersection

You can find the intersections files in the "processed" and "intermediate" folders of your current project ("C:\Users\ "XX"\AppData\Local\pylca\Brightway3\"project_folder").

You can also look at the data directly in bw after performing the lca:

import bw2regional as bwr

lcaX = bwr.TwoSpatialScalesLCA()
lcaX.lci()
lcaX.lcia()

You can see the intersection with the command : lcaX.geo_transform_mm.groups[X].calculate()

You will obtain something like:

(array([  1,   1,   1, ..., 151, 151, 151]), --> Rows of the geotransform matrix
 array([ 1090,  3425,  3695, ..., 10547, 10645, 10714]), --> Columns of the geotransform matrix
 array([7.34311723e+10, 4.41004588e+08, 5.88101152e+09, ...,
        1.45890146e+09, 9.82730954e+10, 4.70504117e+09])) --> data for pairs of row/column

However, you have to know what are the locations corresponding to the rows and columns of the geo transform matrix.

Mapping columns/rows of the geo transform matrix to locations

To do that, first, you use:

lcaX.geo_transform_mm.groups[X].row_mapper.to_dict() --> mapping of the Inventory spatial units
lcaX.geo_transform_mm.groups[X].col_mapper.to_dict() --> mapping of the LCIA spatial units

Second, you build a reverse dictionnary of the locations :

import geopandas as gp
import bw2regional as bwr

geom_mapping = {}
for gc in lca_petrol.get_inventory_geocollections():
    field = bwr.geocollections[gc]["field"]
    gdf = gp.read_file(bwr.geocollections[gc]["filepath"])
    for _, row in gdf.iterrows():
        if gc != "world":
            geom_mapping[(gc, row[field])] = row.geometry
        else:
            geom_mapping[row[field]] = row.geometry
    reversed_geomapping = {v: k for k, v in bd.geomapping.items()}

After that you will be able to see what location correspond to which rows and columns and you will be able to see what is the shared area between the inventory and the LCIA locations.