ST_Intersects slow across two tables

Question

I have a set of municipality maps, with an ID for each municipality, and a Multipolygon for the geometry. In PostGIS, I've put these data into tables. Table 1 has the municipalities in year 2019 (each row is 1 muni with ID and geometry), and Table 2 has the municipalities in year 2005. Spatial indices for both tables. About 400 rows in each table.

Here are some stats on this data: All in the same projection: WGS 84

NPoints | Avg   | Max    | Min
2019:   | 14378 | 379875 | 9
2005:   | 24675 | 556334 | 5

Num Geometries | Avg   | Max    | Min
2019:          |  63   | 4402   | 1
2005:          | 131   | 9480   | 1

ST_Valid: 4 invalid geometries out of about 400 in each set.

Basically what I want to do is look at the intersection and see how similar the two geographies are across time. So, I wrote the query below.

SELECT m1.munid, ST_Intersection(m1.wkb_geometry, m2.wkb_geometry)
INTO intersect_table
FROM m2019 m1
LEFT JOIN m2005 m2 on m1.munid = m2.munid

Essentially I join based on the muni IDs (which don't change across time), then create the intersection. But this query basically doesn't finish, I've left it for 30 minutes and I can't imagine it should take longer. Am I being unreasonable or am I going about this incorrectly?

Answer 1

So I did a redesign in a way, as suggested by @geozelot. Basically, I processed the data using a cursor, row by row, in Python using sql/geoalchemy, and dumped the result into another table, row by row. In particular, when doing the processing, I processed from least number of geometries to most (see the starred line below). So I had most of my answers within minutes, it's just took a few hours to run with the big geometries. But I was able to have incremental results and examine those while I let the computer process in the background.

Not sure how to achieve this in any kind of straight up SQL query that really does all the processing at once.

from sqlalchemy import create_engine
from sqlalchemy.ext.declarative import declarative_base
from sqlalchemy import Column, Integer, String, Float
from geoalchemy2 import Geometry
Base = declarative_base()

engine = create_engine('postgresql://NAME:PASSWORD@localhost/DB_NAME')

from sqlalchemy.orm import sessionmaker
Session = sessionmaker(bind=engine)

class Muni2005(Base):
    __tablename__ = 'm2005r'
    id = Column(Integer, primary_key=True)
    munid = Column(Integer)
    name = Column(String)
    wkb_geometry = Column(Geometry(geometry_type='MULTIPOLYGON', srid=4326))
                       
class Muni2019(Base):
    __tablename__ = 'm2019r'
    id = Column(Integer, primary_key=True)
    munid = Column(Integer)
    name = Column(String)
    wkb_geometry = Column(Geometry(geometry_type='MULTIPOLYGON', srid=4326))

class Inters(Base):
    __tablename__ = 'inters'
    id = Column(Integer, primary_key=True)
    id_2019 = Column(Integer)
    munid = Column(Integer)
    name = Column(String)
    inter = Column(Geometry(srid=4326))
    g2005 = Column(Geometry(srid=4326))
    g2019 = Column(Geometry(srid=4326))


conn = engine.connect()
import time
from sqlalchemy.orm import sessionmaker
from sqlalchemy import func,text

Inters.__table__.drop(engine)
Inters.__table__.create(engine)

Session = sessionmaker(bind=engine)
session = Session()
**munis = session.query(Muni2019).order_by(text("ST_NumGeometries(wkb_geometry)"))**
for muni in munis:
    print muni.name, muni.munid
    match2005 = session.query(Muni2005).filter(Muni2005.munid==muni.munid)
    for muni2005 in match2005:
        i = Inters()
        i.id_2019 = muni.id
        i.munid = muni.munid
        i.name = muni.name
        i.g2005 = muni2005.wkb_geometry
        i.g2019 = muni.wkb_geometry
        i.inter = session.query(func.ST_Intersection(i.g2005, i.g2019))
        session.add(i)
    session.commit()