#!/usr/bin/env python # -*- coding: UTF-8 -*- from statistics import geometric_mean import geopandas as gpd import pandas as pd from sqlalchemy.engine import URL from sqlalchemy import create_engine # from shapely.geometry import Polygon import pycen from os import listdir, chdir from pathlib import Path from zipfile import ZipFile from rasterstats import zonal_stats # Path MNT path0 = '/home/colas/Documents/9_PROJETS/2_PS/' path0_mnt = '/home/colas/Documents/9_PROJETS/3_PGZH/' path = path0_mnt + 'SIG/' p_mltifct = 'multi_fonctions/IGN - BD Alti 25M/' path2_mnt = path0+'MNT/' file_mnt25 = 'MNT_25m.tif' file_mnt5 = 'MNT_5m.tif' # 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) # con = create_engine('postgresql+psycopg2://{0}:{1}@{2}/{3}'.format(user,pwd,adr,base), echo=False) def keepgeom_indept(df): dnat = pycen.ref.get_districtNat().geom.unary_union df['in_isere'] = df.intersects(dnat) df = df[df.in_isere].copy() df.dropna(axis=1, how='all', inplace=True) del df['in_isere'] return df # df = gpd.read_postgis( # sql = 'SELECT * FROM pelouse_seche."cr_PS_CBNA_habitats_aggreg_06_2020"', # con = con # ) lst_tab = ['"cr_VERCORS_habitats_CBNA_1999-2007"','"cr_TRIEVES+VERCORS_habitats_CBNA_2014"','"cr_ECRIN_habitats_CBNA_2014"','"cr_CHARTREUSE_habitats_CBNA_2000-2012"',] table = lst_tab[3] df = gpd.read_postgis( sql = 'SELECT * FROM habitat.%s'%table, con = con ) df.set_index('idfinal', inplace=True) # df.set_index('data_id', inplace=True) colnhab = df.columns[df.columns.str.startswith('n_hab')] nhab = df[colnhab].copy() nhab['SUMM'] = nhab.sum(axis=1) nhab['SUMM'].unique() lst_hab = ['34.1','34.3','34.4','34.7','35.2','62.3','64.1'] df2 = pd.DataFrame() for h in lst_hab: tmp = df[ (df.code_hab1.str.contains(h,na=False))| (df.code_hab2.str.contains(h,na=False)&(df.n_hab2 >= 20))| (df.code_hab3.str.contains(h,na=False)&(df.n_hab3 >= 20))| (df.code_hab4.str.contains(h,na=False)&(df.n_hab4 >= 20)) ] df2 = pd.concat([df2,tmp]) df3 = df2.copy() # Before 1385 / After 1422 lst_hab2 = ['41','42','43','44','81','82','83'] for h2 in lst_hab2: tmp2 = df2[ (df2.code_hab1.str.contains('&'+h2,na=False)&(df2.n_hab1 >= 60)) ] df2.drop(labels=tmp2.index,inplace=True) for h2 in lst_hab2: tmp2 = df2[ (df2.code_hab1.str.contains(h2+'.',na=False)&(df2.n_hab1 >= 60)) ] df2.drop(labels=tmp2.index,inplace=True) for h2 in lst_hab2: tmp3 = df3[ (df3.code_hab1.str.contains(h2,na=False)&(df3.n_hab1 >= 60)) ] df3.drop(labels=tmp3.index,inplace=True) df2 = keepgeom_indept(df2) df3 = keepgeom_indept(df3) df = df2.copy() df.reset_index(inplace=True,drop=False) # ps = pycen.ps.get_sitesGeom() sql = 'SELECT * FROM ps.v_pelouseseches_all' PS = gpd.read_postgis(sql,pycen.con) ps = PS[PS.source==table.replace('"','')].copy() # 1385 not_in_azalee = df2[~df2.index.isin(df3.index)].copy() in_azalee = df2[df2.index.isin(df3.index)].copy() not_in_azalee.to_file(path0+'CBNA/MISSING_DATA_INDB/chartreuse.gpkg',layer='not_in_azalee') in_azalee.to_file(path0+'CBNA/MISSING_DATA_INDB/chartreuse.gpkg',layer='in_azalee') ps.columns # Drop poly who intersect poly non-CBNA # ps_autre = ps[~ps.auteur_geom.str.contains('CBNA')] # intersect = gpd.sjoin(df,ps_autre,op='intersects').idfinal.tolist() # drop_index = df[df.idfinal.isin(intersect)].index # df.drop(drop_index,inplace=True) # Drop poly who intersect poly CBNA # ps_cbna = ps[ps.auteur_geom.str.contains('CBNA')] # intersect = gpd.sjoin( # df, # gpd.GeoDataFrame( geometry=ps_cbna.centroid), # op='intersects').idfinal.tolist() # drop_index = df[df.idfinal.isin(intersect)].index # df.drop(drop_index,inplace=True) from os.path import join PATH = '/home/colas/Documents/9_PROJETS/2_PS/CBNA/MISSING_DATA_INDB' FILE = 'CHARTREUSE.gpkg' df = gpd.read_file(join(PATH,FILE)) if 'date' not in df.columns and 'annee' in df.columns: df.loc[~df.annee.isna(),'date'] = df.loc[~df.annee.isna(),'annee'].astype(str) + '-01-01' df['date'] = pd.to_datetime(df['date']) colnhab2 = df2.columns[df2.columns.str.startswith('n_hab')] nhab2 = df2[colnhab2].copy() nhab2['SUMM'] = nhab2.sum(axis=1) nhab2['SUMM'].unique() if table == '"cr_CHARTREUSE_habitats_CBNA_2000-2012"': dic = { 'JCV':'VILLARET Jean-charles', 'GP' :'PACHE Gilles', 'AM' :'MIKOLAJCZAK Alexis', 'TL' :'LEGLAND Thomas', 'HM' :'MERLE Hugues', 'CC' :'Crassous Claire', # 'AL' # 'LA' # 'JC' } df.obs = df.obs.replace(dic) df.loc[~df.obs.isin(['ONF38',*dic.values()]),'obs'] = 'CBNA' df.loc[df.structur=='ONF 38','structur'] = 'ONF38' if table == '"cr_ECRIN_habitats_CBNA_2014"': df['auteur'] = 'PnE' df['structure'] = 'PnE' df['date'] = '2014-01-01' df['date'] = pd.to_datetime(df['date']) if table == '"cr_TRIEVES+VERCORS_habitats_CBNA_2014"': df['auteur'] = 'CBNA' df['structure'] = 'CBNA' df['date'] = '2014-01-01' df['date'] = pd.to_datetime(df['date']) if table == '"cr_VERCORS_habitats_CBNA_1999-2007"': df['structure'] = 'CBNA' if table in [ '"cr_VERCORS_habitats_CBNA_1999-2007"','"cr_TRIEVES+VERCORS_habitats_CBNA_2014"' ]: df.loc[df.statut=="Habitat d'intérêt communautaire",'statut'] = \ "Communautaire" df.loc[df.statut=="Habitat d'intérêt communautaire retenu prioritaire",'statut'] = \ "Prioritaire" df.loc[df.statut=="Habitat d'intérêt communautaire, retenu prioritaire",'statut'] = \ "Prioritaire" df.loc[ df.statut=="""Habitat d'intérêt communautaire, retenu prioritaire pour les sites riches en orchidées""", 'statut'] = "Prioritaire" df.loc[ df.statut=="""Habitat d'intérêt communautaire retenu prioritaire pour les sites riches en orchidées""", 'statut'] = "Prioritaire" df.loc[ df.statut=="""Habitat communautaire, retenu prioritaire pour les sites riches en orchidées""", 'statut'] = "Prioritaire" df.loc[df.statut=="Habitat d'intérêt communautaire retenu prioritaire",'statut'] = \ "Prioritaire" if table in [ '"cr_CHARTREUSE_habitats_CBNA_2000-2012"', '"cr_ECRIN_habitats_CBNA_2014"']: dic = {'PR':'Prioritaire','IC':'Communautaire','NC':'Non communautaire'} df.statut.replace(dic,inplace=True) df.reset_index(inplace=True, drop=True) df.to_postgis( 'ps'+table[3:-1],pycen.con_bdcen,'pelouse_seche',geom_col='geom',if_exists='replace' ) # [CROSS_MNT] ddf = df[['idfinal','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() # zonal_tmp = zonal_stats(ddf.geom,path2_mnt+file_mnt25,stats="mean max",nodata=0,all_touched=True) # import rasterio # with rasterio.open(path2_mnt+file_mnt5) as src: # affine = src.transform # array = src.read(1) # zonal_tmp3 = zonal_stats(ddf, array, affine=affine,stats="mean max",nodata=0,all_touched=True) # zonal_tmp = pd.DataFrame(zonal_tmp) # zonal_tmp3 = pd.DataFrame(zonal_tmp3) # zonal_tmp.columns = ['max_alti','mean_alti'] # zonal_tmp3.columns = ['max_alti','mean_alti'] # tmp = pd.merge(df,zonal_tmp2,left_index=True,right_index=True) from pycen import con as con_aza sql = 'SELECT * FROM ps.v_pelouseseches_noalti' ddf = gpd.read_postgis(sql,con_aza) zonal_tmp2 = zonal_stats(ddf,path2_mnt+file_mnt5,stats="max",nodata=0,all_touched=True) zonal_tmp2 = pd.DataFrame(zonal_tmp2) zonal_tmp2.columns = ['max_alti'] tmp = pd.concat([ddf,zonal_tmp2], axis=1) tmp = tmp.set_geometry('geom', crs=2154) tmp['infeq_1200'] = tmp.max_alti <= 1200 tmp2 = tmp[['site_code','infeq_1200']].copy() tmp2.rename(columns={'site_code':'id_site'},inplace=True) tmp2.to_sql('r_infeq_1200m',con_aza,'ps',if_exists='append',index=False) df.to_postgis( 'ps'+table[3:-1],pycen.con_bdcen,'pelouse_seche',geom_col='geom',if_exists='replace' ) PSinfeq1200 = tmp[tmp.max_alti <= 1200].copy() PSsup1200 = tmp[tmp.max_alti > 1200].copy() # 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.dropna(axis=1, how='all',inplace=True) if table in [ '"cr_VERCORS_habitats_CBNA_1999-2007"','"cr_TRIEVES+VERCORS_habitats_CBNA_2014"']: PSinfeq1200.loc[PSinfeq1200.statut=="Habitat d'intérêt communautaire",'statut'] = \ "Communautaire" PSinfeq1200.loc[PSinfeq1200.statut=="Habitat d'intérêt communautaire retenu prioritaire",'statut'] = \ "Prioritaire" PSinfeq1200.loc[PSinfeq1200.statut=="Habitat d'intérêt communautaire, retenu prioritaire",'statut'] = \ "Prioritaire" PSinfeq1200.loc[ PSinfeq1200.statut=="""Habitat d'intérêt communautaire, retenu prioritaire pour les sites riches en orchidées""", 'statut'] = "Prioritaire" PSinfeq1200.loc[ PSinfeq1200.statut=="""Habitat d'intérêt communautaire retenu prioritaire pour les sites riches en orchidées""", 'statut'] = "Prioritaire" PSinfeq1200.loc[ PSinfeq1200.statut=="""Habitat communautaire, retenu prioritaire pour les sites riches en orchidées""", 'statut'] = "Prioritaire" PSinfeq1200.loc[PSinfeq1200.statut=="Habitat d'intérêt communautaire retenu prioritaire",'statut'] = \ "Prioritaire" if table in [ '"cr_CHARTREUSE_habitats_CBNA_2000-2012"', '"cr_ECRIN_habitats_CBNA_2014"']: dic = {'PR':'Prioritaire','IC':'Communautaire','NC':'Non communautaire'} PSinfeq1200.statut.replace(dic,inplace=True) DF_INI = PSinfeq1200.copy() df = DF_INI.copy() # df['structure'] = 'CBNA' # sql = """ # DELETE FROM sites.sites # WHERE sites.id in (SELECT id_site from sites.r_sites_geom where id_lot = 1); # """ # with pycen.con.begin() as cnx: # cnx.execute(sql) path0 = '/home/colas/Documents/9_PROJETS/2_PS/TO IMPORT/' cols_date = PSsup1200.columns[PSsup1200.columns.str.contains('date')] PSsup1200[cols_date] = PSsup1200[cols_date].astype(str) PSsup1200.to_file(path0 + 'PS_ECRIN_CBNA_2014_SUP1200.shp') cols_date = DF_INI.columns[DF_INI.columns.str.contains('date')] DF_INI[cols_date] = DF_INI[cols_date].astype(str) DF_INI.to_file(path0 + 'PS_ECRIN_CBNA_2014_INF1200.shp') op = ''' gdalwarp -overwrite -s_srs EPSG:2154 -t_srs EPSG:2154 -co FORMAT=GPKG -of GTiff -tr 25.0 -25.0 -tap -cutline \ "PG:dbname='bd-cen-38' host=192.168.0.189 port=5432 sslmode=disable user='cgeier' password='adm1n*bdCen'" \ -csql "SELECT * FROM habitat."cr_TRIEVES_VERCORS_habitats_CBNA_2014" WHERE idfinal in (24060)" \ "{mnt}" {out} '''.format( mnt=path2_mnt+'out.tif', out=path2_mnt+'24060.tif') system(op) mnt = rio.open(Path_tmp+mnt_out) sql = 'SELECT * FROM pelouse_seche."PB_codehab_nonPresent_dans_corineBiotope"' df = gpd.read_postgis(sql,con) lst = ['cr_TRIEVES+VERCORS_habitats_CBNA_2014','cr_VERCORS_habitats_CBNA_1999-2007','cr_ECRIN_habitats_CBNA_2014','cr_CHARTREUSE_habitats_CBNA_2000-2012'] df = df[df.table_org.isin(lst)] df.loc[df.auteur=='PnE','structure'] = 'PnE' sql = 'SELECT * FROM pelouse_seche."PB_codehabCBNA_nonPresent_dans_corineBiotope"' df2 = gpd.read_postgis(sql,con) df3 = pd.concat([df,df2]) df3.set_geometry('geom',inplace=True,crs=2154) df3.to_postgis( name='PB_codehabCBNA_nonPresent_dans_corineBiotope', con=con, schema='pelouse_seche', if_exists='replace', geom_col='geom' )