Export DVF geolocalisées en Geoparquet

export_dvf_to_geoparquet.py

import gzip
import shutil
import logging
from urllib.request import urlretrieve

import duckdb
from bs4 import BeautifulSoup
import requests

geo_dvf_entry = 'https://files.data.gouv.fr/geo-dvf/latest/csv/'

geo_dvf_entry_content = requests.get(geo_dvf_entry)
soup = BeautifulSoup(geo_dvf_entry_content.text, "html.parser")
years = [a.text.replace('/', '') for a in soup.find_all('a') if a.attrs.get('href') is not None and a.text.startswith('20')]

urls_geodvf_with_year = [[f'https://files.data.gouv.fr/geo-dvf/latest/csv/{year}/full.csv.gz', year] for year in years]

for url_geodvf_with_year in urls_geodvf_with_year:
   urlretrieve(url_geodvf_with_year[0], f'geodvf_{url_geodvf_with_year[1]}.csv.gz')

INSTALL_SPATIAL = "INSTALL spatial;"
LOAD_SPATIAL = "LOAD spatial;"
CREATE_GEODVF_TABLE = f"""CREATE OR REPLACE TABLE geodvf AS
SELECT id_mutation, date_mutation, numero_disposition, nature_mutation,
       valeur_fonciere, adresse_numero, adresse_suffixe, adresse_nom_voie,
       adresse_code_voie, code_postal, code_commune, nom_commune, code_departement,
       ancien_code_commune, ancien_nom_commune, id_parcelle, ancien_id_parcelle,
       numero_volume, lot1_numero, lot1_surface_carrez, lot2_numero,
       lot2_surface_carrez, lot3_numero, lot3_surface_carrez, lot4_numero,
       lot4_surface_carrez, lot5_numero, lot5_surface_carrez, nombre_lots,
       code_type_local, type_local, surface_reelle_bati, nombre_pieces_principales,
       code_nature_culture, nature_culture, code_nature_culture_speciale,
       nature_culture_speciale, surface_terrain, longitude, latitude,
       CAST(SUBSTRING(id_parcelle, 6, 5) AS VARCHAR) AS section_prefixe,
       ST_Point(longitude, latitude) AS geometry
FROM read_csv("geodvf_*.csv.gz",
  quote='"',
  names = ['id_mutation', 'date_mutation', 'numero_disposition', 'nature_mutation', 'valeur_fonciere', 'adresse_numero', 'adresse_suffixe', 'adresse_nom_voie', 'adresse_code_voie', 'code_postal', 'code_commune', 'nom_commune', 'code_departement', 'ancien_code_commune', 'ancien_nom_commune', 'id_parcelle', 'ancien_id_parcelle', 'numero_volume', 'lot1_numero', 'lot1_surface_carrez', 'lot2_numero', 'lot2_surface_carrez', 'lot3_numero', 'lot3_surface_carrez', 'lot4_numero', 'lot4_surface_carrez', 'lot5_numero', 'lot5_surface_carrez', 'nombre_lots', 'code_type_local', 'type_local', 'surface_reelle_bati', 'nombre_pieces_principales', 'code_nature_culture', 'nature_culture', 'code_nature_culture_speciale', 'nature_culture_speciale', 'surface_terrain', 'longitude', 'latitude'],
  types = {{'id_mutation': 'VARCHAR',
'date_mutation': 'DATE',
'numero_disposition': 'INTEGER',
'nature_mutation': 'VARCHAR',
'valeur_fonciere': 'DOUBLE',
'adresse_numero': 'INTEGER',
'adresse_suffixe': 'VARCHAR',
'adresse_nom_voie': 'VARCHAR',
'adresse_code_voie': 'VARCHAR',
'code_postal': 'VARCHAR',
'code_commune': 'VARCHAR',
'nom_commune': 'VARCHAR',
'code_departement': 'VARCHAR',
'ancien_code_commune': 'VARCHAR',
'ancien_nom_commune': 'VARCHAR',
'id_parcelle': 'VARCHAR',
'ancien_id_parcelle': 'VARCHAR',
'numero_volume': 'VARCHAR',
'lot1_numero': 'VARCHAR',
'lot1_surface_carrez': 'DOUBLE',
'lot2_numero': 'VARCHAR',
'lot2_surface_carrez': 'DOUBLE',
'lot3_numero': 'VARCHAR',
'lot3_surface_carrez': 'DOUBLE',
'lot4_numero': 'VARCHAR',
'lot4_surface_carrez': 'DOUBLE',
'lot5_numero': 'VARCHAR',
'lot5_surface_carrez': 'DOUBLE',
'nombre_lots': 'INTEGER',
'code_type_local': 'VARCHAR',
'type_local': 'VARCHAR',
'surface_reelle_bati': 'DOUBLE',
'nombre_pieces_principales': 'INTEGER',
'code_nature_culture': 'VARCHAR',
'nature_culture': 'VARCHAR',
'code_nature_culture_speciale': 'VARCHAR',
'nature_culture_speciale': 'VARCHAR',
'surface_terrain': 'DOUBLE',
'longitude': 'DOUBLE',
'latitude': 'DOUBLE'}}
)
"""

CREATE_GEODVF_INDEX = "CREATE INDEX geodvf_idx ON geodvf USING RTREE (geometry);"

EXPORT_DVF_TO_PARQUET = f"""COPY (
    WITH bbox AS (
        SELECT ST_Extent(ST_Extent_Agg(geometry))::BOX_2D AS b
        FROM geodvf
    )
    SELECT   t.*
    FROM     geodvf AS t
            CROSS JOIN bbox
    ORDER BY ST_Hilbert(t.geometry, bbox.b)
) TO 'geodvf.parquet' (FORMAT 'parquet', COMPRESSION 'zstd', ROW_GROUP_SIZE '100000');"""


output_duckdb_database_name = f"dvf.duckdb"

queries = [INSTALL_SPATIAL, LOAD_SPATIAL, CREATE_GEODVF_TABLE, CREATE_GEODVF_INDEX, EXPORT_DVF_TO_PARQUET]
with duckdb.connect(output_duckdb_database_name) as con:
    for query in queries:
        logging.info(query)
        con.sql(query)

# gdal vector convert geodvf.parquet out.parquet -lco "COMPRESSION=ZSTD" -lco "ROW_GROUP_SIZE=100000" -lco WRITE_COVERING_BBOX=YES
# ogr2ogr out.parquet -lco "COMPRESSION=ZSTD" -lco "ROW_GROUP_SIZE=100000" -lco WRITE_COVERING_BBOX=YES geodvf.parquet
# mv out.parquet geodvf.parquet

# with gzip.open('/tmp/full.csv.gz', 'rb') as f_in:
#     with open('/tmp/geodvf.csv', 'wb') as f_out:
#         shutil.copyfileobj(f_in, f_out)

# ogr2ogr -f Parquet -lco COMPRESSION=ZSTD -lco SORT_BY_BBOX=YES -lco WRITE_COVERING_BBOX=YES /tmp/geodvf_with_bbox.parquet /tmp/geodvf.parquet