Speed up intersection on large polygons with PyQGIS?

Question

Is there a way to speed up my script? I am new to PyQGIS. I am using QGIS 3.16.

What I need: I have 50 polygons (up to 1.675.000 hectares, mean 420.000 ha; one shapefile per polygon) and for each of these areas I need to know which percentage consists of roads, which percentage is water, forest, buildings, etc. The information about the area (roads, buildings, water, etc.) is available in separate layers (.gpkg files) for the entire country; multi-polygon layers.

My current script works for the smaller polygons, but for the larger ones it is just taking too much time to complete.

One idea I have which could help speed things up is to use "Simplify" on the area-information-files but I have no idea which simplification method would be better and what value of tolerance is acceptable?

And maybe there is a better way altogether, perhaps only calculating the total areas of roads, buildings, water etc. within the polygons without actually saving the geographical information? I really only need to know how many hectares within the polygons are of the type road, water, forest, etc.

list = 'f62_2020','m61_2019','f61_2021' # id's of polygons
onedrive = #location of folder on my computer

for b in list: 
    by = b
    output = #location of the map with temporary files
    output2 = #location of the map with output-files (.csv files to be used as input for R)
 
 ## WATER
    processing.run("native: intersection",
        {'INPUT': onedrive + 'Waterdeel_Limburg.gpkg|layername=Waterdeel_Limburg',
         #'OVERLAY': onedrive + 'Werkmap/MCP_def/MCP_' + by + '.shp',
         'OVERLAY': onedrive + 'Werkmap/MCP_gpkg/' + by + '.gpkg|layername=' + by,
         'INPUT_FIELDS': ['typewater'],   # I only need 1 of the columns
         'OVERLAY_FIELDS': ['area2'],
         'OVERLAY_FIELDS_PREFIX': '',
         'OUTPUT': output + 'KLAD_' + by + '_waterdeel.kml'})
    
    # re-calculating area in hectares: 
    processing.run("native: fieldcalculator",
        {'INPUT': output + 'KLAD_' + by + '_waterdeel.kml',
         'FIELD_NAME': 'area2',
         'FIELD_TYPE': 0,
         'FIELD_LENGTH': 14,
         'FIELD_PRECISION': 2,
         'FORMULA': '$area/10000',
         'OUTPUT': output2 + by + '_waterdeel.csv'})
 
    ## TERREIN
    processing.run("native: intersection",
        {'INPUT': onedrive + 'Terrein_Limburg.gpkg|layername=Terrein_Limburg',
         #'OVERLAY': onedrive + 'Werkmap/MCP_def/MCP_' + by + '.shp',
         'OVERLAY': onedrive + 'Werkmap/MCP_gpkg/' + by + '.gpkg|layername=' + by,
         'INPUT_FIELDS': ['typelandgebruik'],
         'OVERLAY_FIELDS': ['area2'],
         'OVERLAY_FIELDS_PREFIX': '',
         'OUTPUT': output + 'KLAD_' + by + '_terrein.kml'})

    processing.run("native: fieldcalculator",
        {'INPUT': output + 'KLAD_' + by + '_terrein.kml',
         'FIELD_NAME': 'area2',
         'FIELD_TYPE': 0,
         'FIELD_LENGTH': 14,
         'FIELD_PRECISION': 2,
         'FORMULA': '$area/10000',
         'OUTPUT': output2 + by + '_terrein.csv'})

# ...same for the other layers with information about the area

Things I have already tried to speed the process up:

• I have created a spatial index to the polygons and saved the shapefile as a geopackage file. The area-information files already had a spatial index.
• For each of the area-information files I made a copy of only part of the country, saved this as a .gpkg and clicked on "Update Extents" on the Source-tab under Properties. And tried using these subset-layers as input for only part of the polygons that lay within this smaller part of the country.
• I used "Clip" first but read that "Intersection" should run faster.
• I have paused synchronisation from my files to the cloud while running my script.