118 lines
3.7 KiB
Python
118 lines
3.7 KiB
Python
#!/usr/bin/env python3
|
|
# -*- coding: UTF-8 -*-
|
|
#Nom : : extract_cadastre.py
|
|
#Description : Extraction du cadastre en fonction du des coordonnées NO,NE,SO,SE
|
|
# ou d'un polygon (.shp).
|
|
#Copyright : 2021, CEN38
|
|
#Auteur : Colas Geier
|
|
#Version : 1.0
|
|
|
|
import pandas as pd
|
|
import geopandas as gpd
|
|
from sqlalchemy import create_engine
|
|
from sqlalchemy.engine import URL
|
|
from geoalchemy2 import Geometry
|
|
# from shapely.geometry.multipolygon import MultiPolygon
|
|
from shapely.geometry.polygon import Polygon
|
|
from pyproj import Transformer
|
|
|
|
|
|
# import contextily as ctx
|
|
# import matplotlib.pyplot as plt
|
|
def geom_map(geom):
|
|
map1 = ctx.providers.OpenTopoMap
|
|
map2 = ctx.providers.Esri.WorldImagery
|
|
maps = geom.to_crs(epsg=3857).copy()
|
|
ax = maps.plot(alpha=0.5, edgecolor='k')
|
|
ctx.add_basemap(
|
|
ax,
|
|
attribution_size=5,
|
|
reset_extent=False,
|
|
source=map1,
|
|
zoom='auto'
|
|
# zoom=14
|
|
)
|
|
ax.set_axis_off()
|
|
|
|
output = '/home/colas/Documents/tmp/AUDE/cadastre_for_Aude.gpkg'
|
|
# Shape : chemin/nom.shp -- None
|
|
shp = None
|
|
# Coordonnées [[NO],[NE],[SE],[SO]] -- None
|
|
lat_point_list = [45.414249,45.414611,45.407813,45.407282] #,45.414249]
|
|
lon_point_list = [ 5.472428, 5.486879, 5.487019, 5.472701] #, 5.472428]
|
|
# bdd
|
|
bdd = True
|
|
code_site = 'CRAS'
|
|
type_zone = None
|
|
|
|
# Parametres bdd
|
|
user = 'cgeier'
|
|
pwd = 'adm1n*bdCen'
|
|
adr = '91.134.194.221'
|
|
base = 'bd_cen'
|
|
epsg = 2154
|
|
crs = 'EPSG:%s'%epsg
|
|
|
|
# Connexion bdd
|
|
url = URL.create('postgresql+psycopg2',
|
|
username=user,
|
|
password=pwd,
|
|
host=adr,
|
|
database=base,
|
|
)
|
|
con = create_engine(url)
|
|
con_open = con.connect()
|
|
|
|
# Shape : chemin/nom.shp -- None
|
|
# shp = None
|
|
# Coordonnées [[NO],[NE],[SE],[SO]] -- None
|
|
# lat_point_list = []
|
|
# lon_point_list = []
|
|
if shp:
|
|
print('IMPORT du mask : %s'%shp)
|
|
polygon = gpd.read_file(shp)
|
|
if polygon.crs.srs != crs.lower():
|
|
polygon.to_crs(epsg, inplace=True)
|
|
elif bdd:
|
|
from pycen import con_bdcen
|
|
gdf = gpd.read_postgis("SELECT * FROM sites.c_sites_zonages WHERE code_site = '%s' --and type_zonage = '%s'"%(code_site,type_zone),con_bdcen)
|
|
polygon = gdf[['geom']]
|
|
polygon.rename_geometry('geometry',inplace=True)
|
|
else:
|
|
if epsg != 2154:
|
|
from pyproj import Transformer
|
|
transformer = Transformer.from_crs(crs, "EPSG:2154", always_xy=True)
|
|
lon_point_list,lat_point_list = transformer.transform(lon_point_list,lat_point_list)
|
|
epsg = 2154
|
|
crs = 'EPSG:%s'%epsg
|
|
polygon_geom = Polygon(zip(lon_point_list, lat_point_list))
|
|
polygon = gpd.GeoDataFrame(index=[0], crs=crs, geometry=[polygon_geom])
|
|
|
|
sql = '''SELECT t1.par_id, t1.geom, t1.codcom, t4.libelle vl,
|
|
t1.ccopre, t1.ccosec, t1.dnupla, t1.dparpi, t5.typprop_lib,t5.typprop, t1.ccocomm,
|
|
t1.ccoprem, t1.ccosecm, t1.dnuplam, t1.type,
|
|
t8.ccodem,
|
|
--substring(t8.dnupro from 3) dnupro,
|
|
--substring(t9.dnuper from 3) dnuper,
|
|
t9.ddenom, t9.jdatnss, t9.dldnss, t9.dsglpm, t9.dlign3, t9.dlign4,
|
|
t9.dlign5, t9.dlign6, t9.dnatpr, t9.gtoper, t9.ccogrm
|
|
FROM cadastre.parcelles_38 t1
|
|
LEFT JOIN cadastre.vl_38 t4 ON (t1.vl_id = t4.vl_id)
|
|
LEFT JOIN cadastre.d_typprop t5 ON (t1.typprop_id = t5.typprop_id)
|
|
LEFT JOIN cadastre.lots_38 t2 ON (t1.par_id = t2.par_id)
|
|
LEFT JOIN cadastre.lots_natcult_38 t3 ON (t2.lot_id = t3.lot_id)
|
|
LEFT JOIN cadastre.cadastre_38 t6 ON (t2.lot_id = t6.lot_id)
|
|
LEFT JOIN cadastre.cptprop_38 t7 ON (t6.dnupro = t7.dnupro)
|
|
LEFT JOIN cadastre.r_prop_cptprop_38 t8 ON (t7.dnupro = t8.dnupro)
|
|
LEFT JOIN cadastre.proprios_38 t9 ON (t8.dnuper = t9.dnuper)
|
|
WHERE ST_Intersects (t1.geom, 'SRID={epsg};{poly}');'''.format(epsg=epsg,poly=polygon.geometry[0])
|
|
df = gpd.read_postgis(
|
|
sql=sql,
|
|
con=con,
|
|
)
|
|
for col in df.columns:
|
|
if (df[col].isna()).all():
|
|
del df[col]
|
|
|
|
df.to_file(output,driver='GPKG')
|