Merge branch 'parcelles_aleatoires' into 'master'

Parcelles aleatoires

See merge request !1

.gitignore

+resources/dl
+resources/private
+output
+*.pyc
+lastest_seed.txt
+MCMC_config.yml

PATBuilder.py

+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import pandas as pd
+import geopandas as gpd
+import os, glob
+import rasterio
+import json
+
+def load_MAJIC(majic_dir):
+    patches = gpd.GeoDataFrame()
+    for f in glob.glob(majic_dir+'/*.shp'):
+        patches.append(gpd.GeoDataFrame.from_file(f))
+    return patches
+
+def build_initial_PAT(rpg_parcelles_filename, rpg_ilots_filename, pat_cultural_classes_filename, majic_dir,
+    bdalti_dir, bdalti_tiles, pat_municipalities_dir, bdv_filename, bdv_fixes_filename,
+    adminexpress_com_filename, adminexpress_epci_filename):
+    '''
+    This method build the initial patches set by consolidating the data sources
+    used by the project:
+     - RPG: Registre Parcellaire Graphique
+     - PAT municipalities: list of csv files containing for each administrative sets (EPCI, PNR), the list of municipalities considered
+     - PAT cultural classes: give the classification of cultural code of RPG considered for the PAT
+     - AdminExpress: datasets containing the shapes and population of municipalities
+     - BDV: list of "Bassins de vie" that aggregates municipalities in functional sets
+     - MAJIC: landed property data, used for extracting "revenu cadastral" for each patch
+
+     Returns an object
+    '''
+    municipalities = build_PAT_municipalities(pat_municipalities_dir, bdv_filename, bdv_fixes_filename,
+        adminexpress_com_filename, adminexpress_epci_filename)
+    municipalities = municipalities[['INSEE_COM','Nom','CODE_EPCI','NOM_EPCI', 'geometry','POPULATION','BV2012','LIBBV2012']]
+    rpg = gpd.GeoDataFrame.from_file(rpg_parcelles_filename, encoding='utf-8')[['ID_PARCEL','CODE_CULTU','geometry']]
+    # we will use the centroid for doing the sjoin with municipalities
+    rpg['patches'] = rpg.geometry
+    rpg.geometry = rpg.geometry.centroid
+    patches = gpd.sjoin(rpg, municipalities, op='intersects') # 'op' is useless because sjoin between points and polygons
+    # sjoin with rpg_ilots for retrieving id of "exploitant"
+    rpg_ilots = gpd.GeoDataFrame.from_file(rpg_ilots_filename, encoding='utf-8')[['id_ilot','id_expl','surf_ilot','geometry']]
+    # restoring geometry instead of centroid
+    patches.geometry = patches['patches']
+    del patches['patches']
+    del patches['index_right']
+    patches_expl = gpd.sjoin(patches, rpg_ilots, how='left', op='intersects')
+    # removing orphan patches (without exploitant)
+    orphan_patches = patches_expl[patches_expl['id_ilot'].isnull()]
+    print('{} patches without "id_expl" ({:.4f} % of total area)'.format(
+        len(orphan_patches),
+        orphan_patches.geometry.area.sum() / patches_expl.geometry.area.sum()
+        ))
+    patches_expl = patches_expl[~patches_expl['id_ilot'].isnull()]
+    patches_expl.geometry = patches_expl.geometry.buffer(0)
+    rpg_ilots.geometry = rpg_ilots.geometry.buffer(0)
+    # selecting the greatest intersection between duplicates patches from sjoin
+    patches_expl['intersection_surf'] = patches_expl.apply(
+       lambda row: rpg_ilots[rpg_ilots['id_ilot']==row['id_ilot']]['geometry'].intersection(row['geometry']).area.sum(),
+       axis=1)
+    patches_expl['max_intersection_surf'] = patches_expl.groupby('ID_PARCEL')['intersection_surf'].transform(max)
+    patches_expl = patches_expl[patches_expl['max_intersection_surf'] == patches_expl['intersection_surf']]
+    # cleaning
+    del patches_expl['max_intersection_surf']
+    del patches_expl['intersection_surf']
+    del patches_expl['index_right']
+    patches_expl['id_ilot'] = patches_expl['id_ilot'].astype(int)
+    patches_expl['id_expl'] = patches_expl['id_expl'].astype(int)
+    # for the last remaining duplicates, we select the "ilot" from the "exploitant" that have the greatest surface
+    patches_expl['expl_surf'] = patches_expl.apply(
+        lambda row: patches_expl[patches_expl['id_expl']==row['id_expl']]['surf_ilot'].sum(),
+        axis=1)
+    patches_expl['max_expl_surf'] = patches_expl.groupby('ID_PARCEL')['expl_surf'].transform(max)
+    patches_expl = patches_expl[patches_expl['max_expl_surf'] == patches_expl['expl_surf']]
+    del patches_expl['max_expl_surf']
+    del patches_expl['expl_surf']
+    # cultural code joining
+    codes = pd.read_csv(pat_cultural_classes_filename)
+    patches_expl = patches_expl.merge(codes[['category','Classe GéoPAT']],how='left',left_on='CODE_CULTU',right_on='category')
+    patches_expl = patches_expl.rename(columns={
+        'Classe GéoPAT':'cultgeopat',
+        'BV2012':'Bdv'
+        })
+    # Adding elevation from BD Alti
+    rasters = [] # reading tiles
+    for tile in bdalti_tiles:
+        dataset = rasterio.open(bdalti_dir+tile+'.asc')
+        rasters.append([dataset, dataset.read(1)])
+    def getElevation(point): # Function for getting the elevation from a point
+        # Looping on every loaded rasters
+        for raster in rasters:
+            coords = raster[0].index(point.x, point.y)
+            # If coords are in bounding box, returning the elevation
+            if coords[0] >= 0 and coords[1] >= 0 and coords[0] < 1000 and coords[1] < 1000:
+                return raster[1][coords]
+    patches_expl['elevation'] = patches_expl.geometry.centroid.apply(getElevation)
+    # TODO join with MAJIC
+    patches_expl['SURF_PARC'] = patches_expl.area
+    return [patches_expl, municipalities]
+
+def build_PAT_municipalities(pat_municipalities_dir, bdv_filename, bdv_fixes_filename,
+        adminexpress_com_filename, adminexpress_epci_filename):
+    # Read list of considered municipalities
+    municipalities = pd.DataFrame()
+    for f in glob.glob(pat_municipalities_dir+'/*.csv'):
+        df=pd.read_csv(f)
+        df[os.path.splitext(os.path.basename(f))[0]] = 'True'
+        if municipalities.empty:
+            municipalities = df
+        else:
+            municipalities = pd.merge(municipalities,df,on=['Code Insee','Nom'],how='outer')
+    municipalities = municipalities.fillna('False')
+    municipalities['Code Insee']=municipalities['Code Insee'].astype(str)
+    # Load "Bassins de vie"
+    bvxls=pd.read_excel(bdv_filename,sheetname='Composition_communale', header=5)
+    municipalities = pd.merge(municipalities,bvxls[['CODGEO','BV2012','LIBBV2012']],left_on='Code Insee',right_on='CODGEO')
+    del municipalities['CODGEO']
+    # fixes for "Bassins de vie"
+    bdvfixes=pd.read_csv(bdv_fixes_filename, delimiter=';').astype(str)[['Code Insee','BV2012','LIBBV2012']]
+    municipalities = pd.merge(municipalities, bdvfixes, how='left', on='Code Insee', suffixes=('_old',''))
+    municipalities['BV2012'] = municipalities['BV2012'].fillna(municipalities['BV2012_old'])
+    municipalities['LIBBV2012'] = municipalities['LIBBV2012'].fillna(municipalities['LIBBV2012_old'])
+    # Load shapes and population
+    admin_express_com = gpd.GeoDataFrame.from_file(
+            adminexpress_com_filename, encoding='utf-8'
+        )[['INSEE_COM','geometry','POPULATION','CODE_EPCI']]
+    municipalities = pd.merge(admin_express_com,municipalities,how='inner',left_on='INSEE_COM',right_on='Code Insee')
+    del municipalities['Code Insee']
+    # load EPCI
+    admin_express_epci = gpd.GeoDataFrame.from_file(adminexpress_epci_filename)[['CODE_EPCI','NOM_EPCI']]
+    municipalities = municipalities.merge(admin_express_epci, how='inner', on='CODE_EPCI')
+    return municipalities
+
+def buildHexagons(patches, muns, hexagonsFilename):
+    muns.crs={'init':'epsg:2154'}
+    hexagons = gpd.GeoDataFrame.from_file(hexagonsFilename)
+    hexagonsPAT=gpd.sjoin(hexagons, muns[['geometry']], op='intersects')
+    # drop duplicates on index
+    hexagonsPAT.reset_index(inplace=True)
+    hexagonsPAT.drop_duplicates(subset='index', inplace=True)
+    # keep only geometry
+    hexagonsPAT.drop(['index','left','bottom','right','top','index_right'], axis=1, inplace=True)
+    hexagonsPAT.reset_index(drop=True, inplace=True)
+    # Get indexes as a columns for keeping indexes after 'overlay' call
+    hexagonsPAT.rename_axis('FID', inplace=True)
+    hexagonsPAT.reset_index(inplace=True)
+    # Intersection between hexagons and patches
+    patches_inter_hexagons = gpd.overlay(patches[['ID_PARCEL','geometry','id_expl','elevation', 'cultgeopat','INSEE_COM']], hexagonsPAT.reset_index(), how="intersection")
+    patches_inter_hexagons.to_file('output/hexagons/patches_inter_hexagons.shp')
+    # computes area of intersection
+    patches_inter_hexagons['area_ha'] = patches_inter_hexagons.area / 10000
+    # Adding elevation mean
+    patches_inter_hexagons['tmp'] = patches_inter_hexagons['elevation']*patches_inter_hexagons['area_ha']
+    grouped = patches_inter_hexagons.groupby('FID')[['tmp','area_ha']].sum()
+    hexagonsPAT['elevation_mean'] = grouped['tmp']/grouped['area_ha']
+    # remove hexagons without intersection with patches
+    hexagonsPAT = hexagonsPAT[~hexagonsPAT['elevation_mean'].isnull()]
+    # Adding number of exploitants
+    hexagonsPAT['nb_expl'] = patches_inter_hexagons.groupby('FID')['id_expl'].nunique()
+    # Adding initial cultural surfaces
+    surfaces = patches_inter_hexagons.groupby('FID').apply(lambda x:x.groupby('cultgeopat')['area_ha'].sum())
+    surfaces = surfaces.reset_index().pivot(index='FID', columns='cultgeopat', values='area_ha')
+    surfaces.columns = ['init_surf_cult' + str(col) for col in surfaces.columns]
+    surfaces.fillna(0, inplace=True)
+    hexagonsPAT = hexagonsPAT.join(surfaces)
+    hexagonsPAT.to_file('output/hexagons/hexagons1000PAT.shp')
+    # Extracting list of municipalities included in each hexagon and the cumulated area of patches
+    areasByHexByMunDataFrame = patches_inter_hexagons.groupby('FID').apply(lambda h: h.groupby('INSEE_COM').apply(lambda p:p.area.sum() / 10000)).to_frame('area_ha')
+    # Converting as a dict
+    areasDict = areasByHexByMunDataFrame.groupby(level=0).apply(lambda dff:dff.xs(dff.name)['area_ha'].to_dict()).to_dict()
+    areasDict['__README']='This dict gives area of patches (in ha) for each municipality in each hexagon. The first key is the ID of the hexagon.'
+    # Writing in file
+    with open("output/hexagons/hexagons_muns.json", "w") as output:
+        json.dump(areasDict, output)
+    return hexagonsPAT
+
+if __name__ == '__main__':
+    import yaml
+    try:
+        from yaml import CLoader as Loader
+    except ImportError:
+        from yaml import Loader
+    resources={}
+    config = yaml.load(open('resources/INDEX.yml','r'), Loader=Loader)
+    for k,v in config.items():
+        resources[v['variable']] = 'resources/'+v['file'] if 'file' in v else k
+    [patches, municipalities] = build_initial_PAT(
+        resources['rpg_parcelles_filename'],
+        resources['rpg_ilots_filename'],
+        resources['pat_cultural_classes_filename'],
+        resources['majic_filename'],
+        resources['bdalti_dir'], config['BDALTI']['tiles'],
+        resources['pat_municipalities_dir'],
+        resources['bdv_filename'], resources['bdv_fixes'], resources['adminexpress_com_filename'],
+        resources['adminexpress_epci_filename'])
+    # TODO merge with MAJIC data not automatic
+    # FIXME use raw data instead these preprocessed data
+    majic = gpd.GeoDataFrame.from_file('Productivite/productivite_shape/Parcelle_PAT_valCad_test.shp', encoding='utf-8')[['ID_PARCEL','VALEUR_CAD']]
+    patches = patches.merge(majic,how='left',on='ID_PARCEL')
+    # some missing data due to the preprocessed data
+    patches = patches[~patches['VALEUR_CAD'].isnull()]
+
+    patches.to_file('output/PAT_patches', encoding='utf-8')
+    gpd.GeoDataFrame(municipalities).to_file('output/municipalities', encoding='utf-8')
+    buildHexagons(patches, municipalities, resources['hexagons_filename'])

Productivite/Parcelle_PAT/Parcelle_PAT.prj

+PROJCS["Lambert_Conformal_Conic",GEOGCS["GCS_GRS 1980(IUGG, 1980)",DATUM["D_unknown",SPHEROID["GRS80",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["standard_parallel_1",44],PARAMETER["standard_parallel_2",49],PARAMETER["latitude_of_origin",46.5],PARAMETER["central_meridian",3],PARAMETER["false_easting",700000],PARAMETER["false_northing",6600000],UNIT["Meter",1]]
\ No newline at end of file

Productivite/Parcelle_PAT/Parcelle_PAT.qpj

+PROJCS["unnamed",GEOGCS["GRS 1980(IUGG, 1980)",DATUM["unknown",SPHEROID["GRS80",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",44],PARAMETER["standard_parallel_2",49],PARAMETER["latitude_of_origin",46.5],PARAMETER["central_meridian",3],PARAMETER["false_easting",700000],PARAMETER["false_northing",6600000],UNIT["Meter",1]]

Productivite/productivite.py

+#!/usr/bin/python3
+# -*- coding: utf-8 -*-
+from shapely.geometry import mapping, shape
+from shapely.ops import cascaded_union
+import fiona
+from fiona import collection
+import math
+
+# Transforme le nom du type de culture en un code utilisé dans l'algorithme
+def return_culture_code(culture):
+	if culture == "Céréales".decode('utf-8') :
+		return 0
+
+	if culture == "Cultures industrielles" :
+		return 1
+
+	if culture == "Fourrages" :
+		return 2
+
+	if culture == "Fruits et légumes".decode('utf-8') :
+		return 3
+
+	if culture == "Oléagineux".decode('utf-8') :
+		return 4
+
+	if culture == "Prairies" :
+		return 5
+
+	if culture == "Protéagineux".decode('utf-8') :
+		return 6
+# Donne la moyenne de la valeur cadastrale au mètre carré du scénario donné pour chacun des types de parcelle
+def indice_productivite(shape_scenario, shape_valeur_cadastrale, statsOutputFilename) :
+	liste_i = {}
+	culture_tab = ["Céréales", "Cultures industrielles", "Fourrages", "Fruits et légumes", "Oléagineux", "Prairies", "Protéagineux" ]
+
+	for i in range(7): # Initialise la liste des parcelles pour chacun des types
+		liste_i[i] = []
+
+	with fiona.open(shape_scenario, 'r') as layer_scenario : # Ouvre la couche scenario
+		for i in range(len(layer_scenario)): # Pour chacune des parcelles du scenario
+
+			if layer_scenario[i]["properties"]["rpg_lfgc_s"] != None : # Si le type de culture est non nul
+				liste_i[return_culture_code(layer_scenario[i]["properties"]["rpg_lfgc_s"])].append(i) # Ajouter son FID dans la liste en fonction de son type
+
+	with fiona.open(shape_valeur_cadastrale, 'r') as layer_valeur_cadastrale: # Ouvre la couche des valeurs cadastrales
+		with open(statsOutputFilename,'w') as statsFile:
+			statsFile.write('type culture;valeur cadastrale moyenne par m²')
+			for y in range(7) : # Pour chacun des types de culture
+				somme_valeur_cadastrale = 0.0 # Initialisation de sa valeur cadastrale
+				for i in liste_i[y]: # Pour chacune des parcelles du type numéro 'y' (voir fonction return_culture_code)
+					valeur_cadastrale =  layer_valeur_cadastrale[i]["properties"]["VALEUR_CAD"]
+					surf_parcelle = shape(layer_valeur_cadastrale[i]["geometry"]).area
+
+					somme_valeur_cadastrale += valeur_cadastrale / surf_parcelle # Incrémentation de la valeur cadastrale en m2
+				if len(liste_i[y]) != 0 : # Si il y a au moins une parcelle possédant le type numéro 'y'
+					statsFile.write("{};{}".format(culture_tab[y],somme_valeur_cadastrale/len(liste_i[y])))
+					print("La moyenne des valeurs cadastrales pour les parcelles de type " + str(culture_tab[y]) + " est de " + str(somme_valeur_cadastrale/len(liste_i[y])) + " par mètre carré")
+
+indice_productivite("./Parcelle_PAT/Parcelle_PAT.shp", "./productivite_shape/Parcelle_PAT_valCad_test.shp", "productivite.csv")

Productivite/productivite.txt

+
+						***  PRODUCTIVITE  *** 
+
+Le critère de productivité est là pour donner l'information relative à la valeur cadastrale des parcelles de chacun des types de culture du PAT. Il doit nous permettre de faire la distinction entre les scénarios où les fruits et légumes sont implantés sur des parcelles 'riches' et ceux où ils sont sur des parcelles moins 'riches'.
+
+Afin d'arriver à ce résultat final, plusieurs étapes sont nécessaires. Afin de traiter les données, nous avons eu besoin de les modifier afin qu'elles reflètent au maximum la réalité des choses. Ce document présentera les différentes étapes que nous avons réalisés et se terminera par la dernière qui correspondra au résultat final que le critère a pour objectif de transmettre.
+
+====================================================================
+
+Etape 1 : Fusion des différentes couches afin d'obtenir une couche recouvrant l'ensemble des parcelles du PAT avec leur valeur cadastrale
+
+Etat initial : 
+Nous avons a notre disposition 3 couches correspondant aux parcelles du territoire PAT pour chacun des départements présents à l'intérieur du territoire.
+42 : Loire, 43 : HAUTE-LOIRE, 63 : PUY-DE-DOME	
+
+Déroulement :
+Un script fût créé afin de regrouper nos 3 couches tout en conservant les données correspondantes aux valeurs cadastrales.
+
+Problèmes : 
+Suite à l'analyse de notre nouvelle couche, nous avons remarqué de multiples imperfections aux niveaux des géométries des entités. Grâce à quelques tests, nous avons pu confirmer la superposition de parcelles cadastrales aux alentours de la séparation des départements. Nous avons donc déduit que la réalisation des couches cadastrales de département n'ont pas été réalisées en se concertant, et qu'il y a donc des parcelles appartenant à 2 départements qui sont présentes dans ces couches.
+
+Résultat : Une couche de la fusion des 3 couches correspondantes aux parcelles PAT de chacun des départements.
+
+Algorithme : 'union_cadastrale.py'
+
+======================================================================
+
+Etape 2 : Rendre la couche précise et nettoyer les problèmes de géométries tout en tenant en compte des variations potentielles de la valeur cadastrale
+
+Etat initial :
+Nous utiliserons la couches des parcelles fusionnées de l'étape 1 et la couche des parcelles du PAT.
+
+Déroulement :
+
+Le concept de cette phase est de nettoyer la couche des parcelles cadastrales de ses superpositions tout en modifiant la valeur cadastrale des entités proportionnellement à la surface perdue. Pour cela, différents choix ont dû être fait lors de la conception de l'algorithme.
+
+La question principale fût la façon de déterminer quelle parcelle couper lors de la superposition de 2 parcelles. Nous avons décidés que le choix de donner une priorité aux parcelles de plus grandes surfaces ou possédant une valeur cadastrale supérieure n'était pas vraiment le meilleur que nous pouvions réaliser et peut s'avérer compliqué dû à la réduction de ces valeurs à chaque tour de boucle. 
+Nous avons donc choisi de découper en fonction de l'ordre d'apparition des entités. Dans le cas où la parcelle 1 intersecte la parcelle 2, la parcelle 1 sera conservée entièrement et la parcelle 2 sera découpée de sa partie empiétant sur la parcelle 1, sa valeur cadastrale sera modifiée en conséquence.
+
+Le principe de l'algorithme tourne autour de ce mécanisme et permet d'obtenir une couche reflétant des résultats plus proches de la réalité.
+
+Résultat :
+Une couche des parcelles cadastrales sans superposition d'entité et avec des valeurs cadastrales modifiées en conséquence.
+
+Algorithme : valeur_cadastrale.py
+
+=====================================================================
+
+Etape 3 : Obtenir une couche des parcelles du PAT avec un champ contenant la vraie valeur cadastrale de chacune des parcelles
+
+Etat initial :
+La couche des parcelles du PAT et la couche résultante de l'étape 2, c'est-à-dire celle des parcelles cadastrales sans les superpositions.
+
+Déroulement :
+Dans cette dernière étape de préparation des données, nous faisons en sorte d'obtenir les valeurs cadastrales des parcelles PAT à l'aide de la couche des parcelles cadastrales nettoyée.
+Nous vérifions donc chacune des parcelles cadastrales intersectant une parcelle PAT. Nous ajoutons ensuite à la parcelle PAT la valeur cadastrale correspondant en fonction du % de la surface de la parcelle cadastrale sur la parcelle PAT.
+Par exemple, si la parcelle cadastrale 'P1' est entièrement incluse dans celle du PAT, celle du PAT se verra ajouter 100% de la valeur cadastrale de la parcelle 'P1'.
+Cependant, si la parcelle 'P2' n'intersecte qu'à 50% la parcelle du PAT, alors il n'y aura qu'un ajout de 50% de la valeur cadastrale de 'P2' dans la parcelle du PAT. 
+
+Résultat : 
+La couche des parcelles du PAT avec 2 nouveaux champs. Le champ 'VALEUR_CAD' contenant la valeur cadastrale de la parcelle et le champ 'SURF_NON_UTILISE' contenant la surface de la parcelle PAT ne possédant pas d'intersection avec des parcelles cadastrales.
+ 
+Algorithme : valeur_cadastrale.py
+
+======================================================================
+
+Dernière partie : Etude d'un scénario pour définir ses valeurs de productivités pour chacun des types de culture
+
+Etat initial : 
+Une couche correspondant à un scénario de répartition des parcelles possibles et la couche résultat de l'étape précédente.
+
+Déroulement :
+Pour commencer, on stocke la valeur cadastrale au mètre carré pour chacun des types de parcelle du PAT. Ensuite, on divise cette somme par le nombre d'occurrence de chacun des types afin d'obtenir la moyenne de la valeur cadastrale pour chaque type.
+
+Résultat :
+Un affichage de la moyenne de la valeur cadastrale par mètre carré pour chacun des types de culture.
+
+Algorithme : 'productivite.py'
+
+
+
+

Productivite/productivite_shape/union_cadastre.py

+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+from shapely.geometry import mapping, shape, Polygon
+from shapely.ops import cascaded_union
+from shapely.wkt import loads
+import fiona
+from fiona import collection
+import random
+from random import randint
+import rtree
+from collections import deque
+
+with fiona.open('./42_jointure_valCad.shp', 'r') as input42:
+	with fiona.open('./43_jointure_valCad.shp', 'r') as input43:
+		with fiona.open('./63_jointure_valCad.shp', 'r') as input63:
+			schema = {'geometry': 'Polygon', 'properties': {'idpar' : 'str', 'idcom' : 'str', 'valeur_cadastrale' : 'int'}}
+			crs = input42.crs
+			driver = input42.driver
+			
+			with fiona.open('./valeurs_cadastrales.shp','w',driver=driver, crs=crs, schema=schema) as file_valeurs_cadastrales :
+
+				for f in input42 :
+					recW = {}
+					recW["geometry"]=f["geometry"]
+					recW["properties"] = {}
+					
+					recW["properties"]["idpar"] = f["properties"]["idpar"]
+					recW["properties"]["idcom"] = f["properties"]["idcom"]
+					recW["properties"]["valeur_cadastrale"] = f["properties"]["inst42_201"]
+					
+					file_valeurs_cadastrales.write(recW)
+					
+				for f in input43 :
+					recW = {}
+					recW["geometry"]=f["geometry"]
+					recW["properties"] = {}
+					
+					recW["properties"]["idpar"] = f["properties"]["idpar"]
+					recW["properties"]["idcom"] = f["properties"]["idcom"]
+					recW["properties"]["valeur_cadastrale"] = f["properties"]["inst43_201"]
+					
+					file_valeurs_cadastrales.write(recW)
+					
+				for f in input63 :
+					recW = {}
+					recW["geometry"]=f["geometry"]
+					recW["properties"] = {}
+					
+					recW["properties"]["idpar"] = f["properties"]["idpar"]
+					recW["properties"]["idcom"] = f["properties"]["idcom"]
+					recW["properties"]["valeur_cadastrale"] = f["properties"]["inst63_201"]
+					
+					file_valeurs_cadastrales.write(recW)
+			
+		
+

Productivite/valeur_cadastrale.py

+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+from shapely.geometry import mapping, shape, Polygon
+from shapely.ops import cascaded_union
+from shapely.wkt import loads
+from shapely.validation import explain_validity
+import fiona
+from fiona import collection
+import random
+from random import randint
+import rtree
+import re
+
+
+
+with fiona.open('./Parcelle_PAT/Parcelle_PAT.shp', 'r') as inputShape:
+
+	schema = {'geometry': 'Polygon', 'properties': {'ID_PARCEL' : 'str', 'INSEE_COM' : 'str', 'CODE_CULTU' : 'str', 'Bdv' : 'str', 'rpg_lfgc_s' : 'str' , 'SURF_PARC' : 'float', 'VALEUR_CAD' : 'float', 'SURF_NON_U' : 'float'}}
+	schemaCadastre = {'geometry' : 'Polygon', 'properties': {'idpar' : 'str', 'idcom' : 'str', 'valeur_cad' : 'float'}}
+	crs = inputShape.crs
+	driver = inputShape.driver
+
+	#m = re.search('(.*)(.shp)', './Parcelle_PAT/Parcelle_PAT.shp')
+
+
+	# PARTIE 1 : Transforme la couche de valeur cadastrale contenant des erreurs de géométries et des superpositions en couche de meilleure qualité
+
+	# Ouvre la couche de la fusion des 3 départements existants dans le PAT
+	with fiona.open('./productivite_shape/valeur_cadastrale_propre.shp' ,'w',driver=driver, crs=crs, schema=schemaCadastre) as shape_cadPropre:
+		with fiona.open('./Parcelle_PAT/Parcelle_PAT.shp', 'r') as input:
+			with fiona.open('./productivite_shape/valeurs_cadastrales.shp', 'r') as inputCadastre : # Ouvre la couche des parcelles cadastrales contenant leurs valeurs
+				# Création de l'index spatial sur les parcelles cadastrales
+				index = rtree.index.Index()
+				i = 0
+				for feat1 in inputCadastre : # Pour chacune des géométries
+					if feat1['geometry'] != None : # Si la géométrie est non nulle
+						geom1 = shape(feat1['geometry'])
+						index.insert(i, geom1.bounds)
+					i += 1
+
+				liste_parcelles_traitees = [] # Liste qui contiendra le FID des parcelles déjà traitées
+				for i in range(len(input)): # Nouvelle boucle traversant toutes les parcelles du PAT
+							geom2 = shape(input[i]["geometry"]) # Stocke la géométrie de la parcelle traversée
+
+
+							for fid in list(index.intersection(geom2.bounds)) : # Pour chacune des parcelles cadastrales ayant une intersection avec la parcelle du PAT
+								if inputCadastre[fid]["geometry"] != None and fid not in liste_parcelles_traitees : # Si la parcelle n'a pas déjà été traitée et que sa geométrie est non nulle
+									feat1 = inputCadastre[fid]
+									geom1 = shape(feat1["geometry"])
+
+									liste_parcelle_intersects = [] # Liste qui contiendra toutes les parcelles intersectant la parcelle de fid 'fid2'
+									for fid2 in list(index.intersection(geom1.bounds)) : # Pour les parcelles intersectant les limites de cette parcelle
+										if inputCadastre[fid2]["geometry"] != None  and fid != fid2: # S'il possède une géométrie et qu'il n'est pas lui-même
+											feat12 = inputCadastre[fid2]
+											geom12 = shape(feat12["geometry"])
+											if geom1.intersection(geom12).area > 0 : # Si la surface d'intersection est supérieure à 0, donc si il y a superposition
+												if fid not in liste_parcelles_traitees : # Si la parcelle 'fid' n'a pas déjà été traitée
+													liste_parcelles_traitees.append(fid) # L'ajouter à la liste des parcelles traitées
+												liste_parcelle_intersects.append(fid2) # Ajouter 'fid' à la liste des parcelles intersectant 'fid'
+
+									if len(liste_parcelle_intersects) > 0 : # Si 'fid' intersecte 1 ou plusieurs parcelles
+										geom1Bis = feat1["geometry"]
+										surface = geom1.area
+										valeur_cadastraleBis = inputCadastre[fid]["properties"]["valeur_cad"]
+										newGeom = geom1Bis
+										for y in range(len(liste_parcelle_intersects)) : # Pour chacune des parcelles intersectant 'fid'
+											if liste_parcelle_intersects[y] in liste_parcelles_traitees :  # Si la parcelle a déjà été traitée
+												geomY = shape(inputCadastre[liste_parcelle_intersects[y]]["geometry"]) # Stocker sa géométrie
+												surfaceIntersect = shape(newGeom).intersection(geomY).area # Stocker la surface d'intersection
+
+												# La nouvelle géometrie correspondra à l'ancienne sans la partie intersectant l'autre parcelle déjà traitée
+												newGeom = mapping(shape(geom1Bis).difference(geomY))
+
+												if newGeom["type"] != 'GeometryCollection' : # Si la nouvelle géométrie est vide (superposition complète de la parcelle 'fid2' sur 'fid')
+													geom1Bis = newGeom
+													valeur_cadastraleBis = valeur_cadastraleBis * (surface - surfaceIntersect) / surface
+
+										# Recuperer le fid des parcelles dans liste_parcelles_intersects pour recuperer la valeur des properties
+										recW = {}
+										recW["geometry"] = geom1Bis
+
+										recW["properties"] = {}
+										recW["properties"]["idpar"] = inputCadastre[fid]["properties"]["idpar"]
+										recW["properties"]["idcom"] = inputCadastre[fid]["properties"]["idcom"]
+										recW["properties"]["valeur_cad"] = valeur_cadastraleBis
+
+										shape_cadPropre.write(recW) # Ajoute cette entitée dans la nouvelle couche
+
+				# Ajouter toutes les parcelles non traitées à la couche. Ces parcelles ne possèdent aucun problème d'intersection pouvant poser des problèmes potentiels
+				cpt = 0
+				for i in range(len(inputCadastre)) :
+					if i not in liste_parcelles_traitees :
+						#print cpt
+						cpt += 1
+						recW = {}
+						recW["geometry"] = inputCadastre[i]["geometry"]
+
+						recW["properties"] = {}
+						recW["properties"]["idpar"] = inputCadastre[i]["properties"]["idpar"]
+						recW["properties"]["idcom"] = inputCadastre[i]["properties"]["idcom"]
+						recW["properties"]["valeur_cad"] = inputCadastre[i]["properties"]["valeur_cad"]
+
+						shape_cadPropre.write(recW)
+
+	# PARTIE 2 : Utilise la couche de valeur cadastrale de bonne qualité afin de l'associer avec la couche des parcelles du PAT
+	# Cela résulte en une couche des parcelles du PAT avec un champ contenant leurs valeurs cadastrales tenant compte des superpositions des parcelles cadastrales avec celles du PAT
+
+	# Création d'index spatial sur les parcelles du PAT
+	index = rtree.index.Index()
+	liste_valeur_cadastrale = {}
+	liste_surface_non_utilise_PAT = {}
+	i = 0
+	for feat1 in inputShape :
+		geom1 = shape(feat1['geometry'])
+		index.insert(i, geom1.bounds)
+		# Champ de la valeur cadastrale des parcelles PAT
+		liste_valeur_cadastrale[i] = 0
+		# Champ de la surface de la parcelle ne possédant pas de valeur cadastrale connu avec nos données MAJIC
+		liste_surface_non_utilise_PAT[i] = geom1.area
+
+		i += 1
+
+	with fiona.open('./productivite_shape/Parcelle_PAT_valCad_test.shp' ,'w',driver=driver, crs=crs, schema=schema) as shape_cad:
+		with fiona.open('./productivite_shape/valeur_cadastrale_propre.shp', 'r') as inputCadastre : # shape des parcelles cadastrales et de leurs valeurs
+
+			i = 0
+			for i in range(len(inputCadastre)): #Pour chacune des parcelles
+				print(i)
+				if inputCadastre[i]["geometry"] != None : # Si il y a une géométrie
+					geom2 = shape(inputCadastre[i]["geometry"]).buffer(0) # Faire un buffer 0 pour corriger des erreurs de géométries potentielles
+
+					for fid in list(index.intersection(geom2.bounds)): # Pour chaque parcelle de l'index i qui intersecte les limites d'une parcelle cadastrale
+						feat1 = inputShape[fid]
+						geom1 = shape(feat1["geometry"])
+
+						if (geom2).intersects(geom1) : # Si les 2 géométries s'intersectent
+
+							surface_intersection = geom2.intersection(geom1).area # Surface d'intersection
+
+							# Incrémentation de la valeur cadastrale de la parcelle avec fonction du % de superposition entre la parcelle du PAT et celle du cadastre
+							liste_valeur_cadastrale[fid] += inputCadastre[i]["properties"]["valeur_cad"] * surface_intersection / geom2.area
+
+							# Décrémentation permettant de connaître la surface d'une parcelle PAT ne contenant pas de valeur cadastrale
+							liste_surface_non_utilise_PAT[fid] -= float(surface_intersection)
+				i += 1
+
+		# Ajout des parcelles PAT avec les nouveaux champs contenant la valeur cadastrale et la surface de la parcelle n'ayant pas eu d'intersection avec une parcelle cadastrale
+		for i in range(len(liste_valeur_cadastrale)) :
+			recW = {}
+			recW["geometry"]=inputShape[i]["geometry"]
+			recW["properties"] = {}
+			recW["properties"]["ID_PARCEL"] = inputShape[i]["properties"]["ID_PARCEL"]
+			recW["properties"]["INSEE_COM"] = inputShape[i]["properties"]["INSEE_COM"]
+			recW["properties"]["CODE_CULTU"] = inputShape[i]["properties"]["CODE_CULTU"]
+			recW["properties"]["Bdv"] = inputShape[i]["properties"]["Bdv"]
+			recW["properties"]["rpg_lfgc_s"] = inputShape[i]["properties"]["rpg_lfgc_s"]
+			recW["properties"]["SURF_PARC"] = inputShape[i]["properties"]["SURF_PARC"]
+			recW["properties"]["VALEUR_CAD"] = liste_valeur_cadastrale[i]
+			recW["properties"]["SURF_NON_U"] = liste_surface_non_utilise_PAT[i]
+			shape_cad.write(recW)

download_data.py

+import yaml
+import urllib.request