Commit e2cf2e80 authored by Seneque Colin's avatar Seneque Colin Committed by Dumoulin Nicolas
Browse files

Critère de productivité et correction du MCMC

parent c031074c
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
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]]
#!/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 rtree
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) :
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
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'
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("../Etape1_shp_propre/MCMC/testSimulationSurf_MCMC.shp", "./productivite_shape/Parcelle_PAT_valCad.shp")
*** 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'
#!/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)
#!/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:
f
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)
......@@ -9,6 +9,7 @@ import random
from random import randint
from proximite import proximite
# 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
......@@ -31,7 +32,7 @@ def return_culture_code(culture):
if culture == "Protéagineux".decode('utf-8') :
return 6
# Fonction de création d'un scénario aléatoire modifiant un nombre donné en paramètre 'nb_cult_modif' de parcelles
def MCMC_creation(shape, nb_cult_modif, cereales, culturesI, fourrages, fruit_legumes, oleagineux, prairies, proteagineux):
with fiona.open(shape, 'r') as input:
schema = input.schema.copy()
......@@ -40,8 +41,7 @@ def MCMC_creation(shape, nb_cult_modif, cereales, culturesI, fourrages, fruit_le
resultat_final = []
#Correspond au nombre de parcelle de culture que le programme modifiera
#Parametre
#Correspond au nombre de parcelles de culture que le programme modifiera
nombre_culture_a_modifie = nb_cult_modif
#Surfaces souhaitées pour chacun des types de culture pour la simulation
......@@ -59,7 +59,7 @@ def MCMC_creation(shape, nb_cult_modif, cereales, culturesI, fourrages, fruit_le
#Initialisation des surfaces de chaque culture à 0, ceci incrémentera à la création de chaque parcelle lors de la simulation
cultures = {}
for i in range (7) :
for i in range (8) :
cultures[i] = { 'surf': 0.0, 'parcelles': []}
......@@ -68,15 +68,19 @@ def MCMC_creation(shape, nb_cult_modif, cereales, culturesI, fourrages, fruit_le
i = 0
for f in input :
total_culture += f["properties"]["SURF_PARC"]
total_culture += float(f["properties"]["SURF_PARC"])
#Initialisation de la surface des cultures de la couche d'entrée
# Initialisation de la surface des cultures de la couche d'entrée
culture = f["properties"]["rpg_lfgc_s"]
if culture != None :
#total_culture += float(f["properties"]["SURF_PARC"])
cultures[return_culture_code(culture)]['surf'] += float(f["properties"]["SURF_PARC"])
# Le FID de la parcelle est ajouté à la liste 'cultures[typeCulture][parcelles]
cultures[return_culture_code(culture)]['parcelles'].append(i)
else :
cultures[7]['parcelles'].append(i)
cultures[7]['surf'] += float(f["properties"]["SURF_PARC"])
#Incrémentation du fid des parcelles
i += 1
......@@ -88,20 +92,31 @@ def MCMC_creation(shape, nb_cult_modif, cereales, culturesI, fourrages, fruit_le
liste_cultures_a_augmenter.append(i)
else :
liste_cultures_a_reduire.append(i)
#print cultures[i]['surf']*100/total_culture
culture_tab = ["Céréales", "Cultures industrielles", "Fourrages", "Fruits et légumes", "Oléagineux", "Prairies", "Protéagineux" ]
with fiona.open('MCMC/testSimulationSurf_MCMC.shp','w',driver=driver, crs=crs, schema=schema) as parcelles :
cpt = 0
while cpt < nombre_culture_a_modifie and len(liste_cultures_a_augmenter) > 0 :
while cpt < nombre_culture_a_modifie and len(liste_cultures_a_augmenter) + len(liste_cultures_a_reduire) > 1 :
#Récupération aléatoire d'une parcelle possédant une culture devant diminuer (stocké dans le dictionnaire 'cultures[i]["parcelles"]'
rand_cul = random.choice(liste_cultures_a_reduire)
if len(liste_cultures_a_reduire) != 0 :
rand_cul = random.choice(liste_cultures_a_reduire)
else :
# Cas où les proportions des parcelles sont déjà optimisé mais où il reste des parcelles dans plusieurs types
# On ajoute un type aléatoire des types à augmenter (qui sera extrèmement proche de la valeur souhaité) dans celle des valeurs à réduire afin de continuer l'algorithme jusqu'à réduire au maximum le nombre de parcelle restante
# Permet d'optimiser au maximum la randomisation du résultat
rand_cul = random.choice(liste_cultures_a_augmenter)
liste_cultures_a_augmenter.remove(rand_cul)
liste_cultures_a_reduire.append(rand_cul)
random_parcelles= random.choice(cultures[rand_cul]['parcelles'])
if len(liste_cultures_a_augmenter) == 0 :
rand_augmenter = random.choice(liste_cultures_a_reduire)
liste_cultures_a_augmenter.append(rand_augmenter)
liste_cultures_a_reduire.remove(rand_augmenter)
random_parcelles = random.choice(cultures[rand_cul]['parcelles'])
code_culture_rand = input[random_parcelles]["properties"]["rpg_lfgc_s"]
......@@ -118,27 +133,34 @@ def MCMC_creation(shape, nb_cult_modif, cereales, culturesI, fourrages, fruit_le
while succes == 0 :
rand_cult = random.choice(liste_cultures_a_augmenter)
cultures[rand_cult]['surf'] += input[random_parcelles]["properties"]["SURF_PARC"]
recW["properties"]['rpg_lfgc_s'] = culture_tab[rand_cult]
cultures[rand_cult]['surf'] += float(input[random_parcelles]["properties"]["SURF_PARC"] )
recW["properties"]['rpg_lfgc_s'] = culture_tab[rand_cult]
#Retirer la parcelle 'random_parcelles' de 'liste_nombre'
cultures[return_culture_code(input[random_parcelles]["properties"]["rpg_lfgc_s"])]['parcelles'].remove(random_parcelles)
#Retirer la surface de l'ancienne valeur de la parcelle contenu dans 'culture[type_culture]'
cultures[return_culture_code(input[random_parcelles]["properties"]["rpg_lfgc_s"])]['surf'] -= float(input[random_parcelles]["properties"]["SURF_PARC"])
parcelles.write(recW)
resultat_final.append(recW)
# Si la proportion de surface d'une culture du territoire dépasse la proportion souhaitée
# On l'a retire de la liste des cultures à augmenter pour la mettre dans celle à réduire
if cultures[rand_cult]['surf']*100/total_culture >= pourcentage_surface[rand_cult] :
if float(cultures[rand_cult]['surf']*100/total_culture) >= pourcentage_surface[rand_cult] :
liste_cultures_a_augmenter.remove(rand_cult)
liste_cultures_a_reduire.append(rand_cult)
# Pareil mais l'inverse pour les cultures qu'il ne faut plus réduire mais qu'il faut augmenter
if cultures[return_culture_code(input[random_parcelles]["properties"]["rpg_lfgc_s"])]['surf'] <= pourcentage_surface[rand_cult] :
if float(cultures[return_culture_code(input[random_parcelles]["properties"]["rpg_lfgc_s"])]['surf']*100/total_culture) <= pourcentage_surface[return_culture_code(input[random_parcelles]["properties"]["rpg_lfgc_s"])] :
liste_cultures_a_reduire.remove(return_culture_code(input[random_parcelles]["properties"]["rpg_lfgc_s"]))
liste_cultures_a_augmenter.append(return_culture_code(input[random_parcelles]["properties"]["rpg_lfgc_s"]))
succes = 1
cpt += 1
# Décommenter pour voir en 'live' les changements des proportions des types dans la console
#for i in range(7) :
#print cultures[i]['surf']*100/total_culture
# Si toutes les parcelles d'un type de culture ont déjà été modifiées, on retire ce type de culture des listes
for i in range(7):
if len(cultures[i]["parcelles"]) == 0 :
......@@ -146,8 +168,9 @@ def MCMC_creation(shape, nb_cult_modif, cereales, culturesI, fourrages, fruit_le
liste_cultures_a_reduire.remove(i)
elif i in liste_cultures_a_augmenter :
liste_cultures_a_augmenter.remove(i)
for i in range(7) :
# Toutes les parcelles restantes sont ajoutées dans la nouvelle couche, les parcelles avec un type vide non traitées précédement sont ajoutées dans cette phase
for i in range(8) :
for e in cultures[i]["parcelles"] :
recW = {}
recW["geometry"]=input[e]["geometry"]
......@@ -158,27 +181,27 @@ def MCMC_creation(shape, nb_cult_modif, cereales, culturesI, fourrages, fruit_le
recW["properties"]['Bdv']=input[e]["properties"]["Bdv"]
recW["properties"]['INSEE_COM']=input[e]["properties"]["INSEE_COM"]
recW["properties"]['rpg_lfgc_s']=input[e]["properties"]["rpg_lfgc_s"]
#parcelles.write(recW)
parcelles.write(recW)
resultat_final.append(recW)
print
print "Pourcentage culture apres : "
print ("Pourcentage culture apres MCMC : ")
for i in range(7) :
print cultures[i]['surf']*100/total_culture
print (cultures[i]['surf']*100/total_culture)
return resultat_final
#MCMC_creation(shape, nb_cult_modif, cereales, culturesI, fourrages, fruit_legumes, oleagineux, prairies, proteagineux)
cereales = 20
cereales = 10
culturesI = 5
fourrages = 10
fruits_legumes = 15
fruits_legumes = 10
oleagineux = 10
prairies = 15
proteagineux = 25
prairies = 45
proteagineux = 10
nombre_parcelle_a_modifier = 78000
liste_resultat = MCMC_creation('./MCMC/Parcelle_PAT_MCMC.shp', nombre_parcelle_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux)
proximite(liste_resultat, fruits_legumes)
#proximite(liste_resultat, fruits_legumes)
*** MCMC ***
L'algorithme MCMC a pour objectif de modifier un scénario initial pour le transformer de façon partiellement aléatoire.
En paramètre, il est possible de manipuler la proportion de céréales, cultures industrielles, fourrages, fruits et légumes, oléagineux, prairies et protéagineux que l'ont souhaite avoir dans notre nouveau scénario.
Il y a également la possibilité de définir le nombre de parcelles que l'ont souhaite modifier.
==============================================
Mode de fonctionnement :
- Stocker l'état initial du scénario dans un dictionnaire.
- Définir quels sont les types de culture qui doivent gagner en surfaces et ceux qui doivent en perdre.
- Remplacer les parcelles à réduire par celles à augmenter afin de se rapprocher au maximum des proportions définies en paramètre.
- Remplacer jusqu'à ce qu'on arrive au nombre de parcelle à modifier défini en paramètre
Résultat :
Une nouvelle couche contenant un scénario se rapprochant au maximum des proportions définies en paramètre
Algorithme : 'MCMC.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 rtree
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') :