pyqgis_workshop_foss4g_2023.py

from csv import DictReader
from datetime import datetime, date
from functools import reduce
from os.path import join

work_dir = "/Users/gis/Desktop/pyqgis_workshop_foss4g2023"

gpkg_path = join(work_dir, "naturalearth.gpkg")
layer_name = "ne_10m_admin_1_states_provinces"

layer_path = f"{gpkg_path}|layername={layer_name}"   
layer = QgsVectorLayer(layer_path, layer_name, "ogr")
layer.setSubsetString("iso_a2='IT'")

# Iterating over features.

prj_layers = QgsProject.instance().mapLayersByName(layer_name)
if len(prj_layers) == 0:
    QgsProject.instance().addMapLayer(layer)

features = layer.getFeatures()

region_dict = [{"name": f.attribute("region"), "geom": f.geometry()}
                for f in features]
region_names = set([f["name"] for f in region_dict])

regions_grouped = {n: [f["geom"] for f in region_dict if f["name"] == n] for n in region_names}

regions_merged = {k: reduce(lambda a, b: a.combine(b), v, v[0]) for k, v in regions_grouped.items()}

with open(join(work_dir, "dpc-covid19-ita-regioni.csv"), "r") as file:
    reader = DictReader(file)

    covid_rows = [{"date": datetime.fromisoformat(f["\ufeffdata"]).date(),
                  "total": int(f["totale_casi"]),
                  "geom": QgsGeometry(QgsPoint(float(f["long"]), float(f["lat"]))),
                  }
                  for f in reader]


days = set([d["date"] for d in covid_rows])

feature_days_dict = {}
for day in list(days)[:10]:
    filter_rows = [r for r in covid_rows if r["date"] == day]

    features = []
    for region_name, region_geom in regions_merged.items():
        cases_filtered = [f for f in filter_rows if region_geom.intersects(f["geom"])]
        
        if len(cases_filtered) > 1:
            print(region_name, cases_filtered)

        # Provisional fix.
        total_cases_filtered = sum([f["total"] for f in cases_filtered])
        # Create feature.
        feature = QgsFeature()
        feature.setGeometry(region_geom)
        feature.setAttributes([region_name, day.isoformat(), total_cases_filtered])

        features.append(feature)
    feature_days_dict[day.isoformat()] = features


for day, features_list in feature_days_dict.items():
    new_layer_name = "covid_italy"

    output_layer = QgsVectorLayer("Polygon?crs=EPSG:4326&field=day:string&field=region_name:string&field=total_cases:int", new_layer_name, "memory")
    output_layer.dataProvider().addFeatures(features_list)
    
    render_attributes = [
        [0, 1000, "yellow"],
        [1001, 3000, "green"],
        [3001, 10000, "blue"],
        [10001, 40000, "orange"],
        [40001, 1000000, "red"],
        [1000001, 10000000, "black"],
    ]

    renderer = QgsGraduatedSymbolRenderer('total_cases')

    for (min, max, color) in render_attributes:
        classification_range = QgsClassificationRange(f'{min} <= covid_cases <= {max}', min, max)
        properties = {
            "color": color,
            "outline_color": "white",
            "outline_width": 1,
        }
        symbol = QgsFillSymbol.createSimple(properties)
        render_range = QgsRendererRange(classification_range, symbol)
        renderer.addClassRange(render_range)

    output_layer.setRenderer(renderer)
    QgsProject.instance().addMapLayer(output_layer)
    
    layout = QgsPrintLayout(QgsProject.instance())
    layout.initializeDefaults()
    
    map = QgsLayoutItemMap(layout)
    map.attemptMove(QgsLayoutPoint(5, 25, QgsUnitTypes.LayoutMillimeters))
    map.attemptResize(QgsLayoutSize(285, 180, QgsUnitTypes.LayoutMillimeters))
    
    map.setExtent(QgsRectangle(0, 0, 285, 180))
    map.zoomToExtent(output_layer.extent())
    
    map.setFrameEnabled(True)
    
    layout.addLayoutItem(map)
    
    exporter = QgsLayoutExporter(layout)

    path = join(work_dir, f"output_{day}.png")
    exporter.exportToImage(path, QgsLayoutExporter.ImageExportSettings())