Python_scripts/cross_MNT.py

89 lines
2.2 KiB
Python

#!/usr/bin/env python
# -*- coding: UTF-8 -*-
# from binascii import Incomplete
from warnings import filterwarnings
import geopandas as gpd
import pandas as pd
from sqlalchemy import create_engine
from sqlalchemy.engine import URL
# from shapely import geometry
# from shapely.geometry import MultiPolygon, MultiLineString, MultiPoint, Polygon, LineString #, collection
# from shapely import ops
from os import listdir, chdir
from pathlib import Path
from zipfile import ZipFile
from rasterstats import zonal_stats
# Test de LiveShare DM
path0 = '/home/colas/Documents/9_PROJETS/3_PGZH/'
path = path0 + 'SIG/'
p_mltifct = 'multi_fonctions/IGN - BD Alti 25M/'
# Parametres bdd IN
user = 'cen_admin'
pwd = '#CEN38@venir'
adr = '192.168.0.189'
base = 'bd-cen-38'
url = URL.create('postgresql+psycopg2',
username=user,
password=pwd,
host=adr,
database=base,
)
con = create_engine(url)
df = gpd.read_postgis(
sql = 'SELECT * FROM pelouse_seche."cr_PS_CBNA_habitats_aggreg_06_2020"',
con = con
)
ddf = df[['id','geom']].copy()
ddf.set_geometry('geom',inplace=True)
home = str(Path.home())
chdir(path+p_mltifct)
Dir = listdir()
Dir = [x for x in Dir if '.zip' in x]
# stats = pd.DataFrame()
for i, d in enumerate(Dir):
zip = ZipFile(d).namelist()
z = [z for z in zip if 'MNT' in z][0]
print(z)
tmp = zonal_stats(ddf.geom,'/vsizip/{zip}/{mnt}'.format(zip=d,mnt=z),stats="max")
tmp = pd.DataFrame(tmp)
tmp.columns = ['max%s'%i]
if i == 0 :
stats = tmp
else:
stats = pd.merge(stats,tmp,how='left',left_index=True,right_index=True)
stats['bilan'] = stats.max(axis=1)
tmp = pd.merge(df,stats['bilan'],left_index=True,right_index=True)
tmp.bilan = tmp.bilan.round()
PSinfeq1200 = tmp[tmp.bilan <= 1200].copy()
PSsup1200 = tmp[tmp.bilan > 1200].copy()
PSinfeq1200.rename(columns={'bilan':'max_alti'}, inplace=True)
PSsup1200.rename(columns={'bilan':'max_alti'}, inplace=True)
PSinfeq1200.to_postgis(
name='cr_PS_CBNA_habitats_inf-eq_1200m',
con=con,
schema='pelouse_seche',
if_exists='replace',
index=False,
index_label='id',
geom_col='geom'
)
PSsup1200.to_postgis(
name='cr_PS_CBNA_habitats_sup_1200m',
con=con,
schema='pelouse_seche',
if_exists='replace',
index=False,
index_label='id',
geom_col='geom'
)