#!/usr/bin/python # -*- coding: utf-8 -*- import fiona import numpy as np import re from shapely.geometry import mapping, shape from shapely.ops import cascaded_union from fiona import collection import geopandas as gpd import random from random import randint import json import gzip import sqlite3 import sys, os import yaml from proximite import Proximite from resilience_list import Resilience from productivite import Productivity #from indice_biodiversite import biodiversite from indice_biodiversite_2 import Biodiversity from social import Social from EPCI_culture_repartition import EPCI_culture_repartition from radar import scenarios_graph class Indicators: def __init__(self, mcmc_config_filename, initial_patches): if not os.path.isfile(mcmc_config_filename): print('Error: file not found "{}"'.format(mcmc_config_filename)) print('Please copy the template file "MCMC_config.sample.yml" and adjust to your settings and run again this program') sys.exit(1) self.mcmc_config = yaml.load(open(mcmc_config_filename,'r')) self._proximity = Proximite(initial_patches, self.mcmc_config['indicators_config']['proximite']['surf_totale_cible']) self._resilience = Resilience(self.mcmc_config['indicators_config']['resilience']) self._productivity = Productivity() biodiv_config = self.mcmc_config['indicators_config']['biodiversity'] matrixfilename = biodiv_config['matrixfilename'] self._biodiversity = Biodiversity(initial_patches, biodiv_config['dmax'], biodiv_config['epsilon']) self._biodiversity.load_probabilities(matrixfilename) social_config = self.mcmc_config['indicators_config']['social'] self._social = Social(initial_patches, social_config['cost_matrix_filename'], social_config['bdv_threshold']) def proximity(self, patches): return self._proximity.compute_indicator(patches, False) def resilience(self, patches): return self._resilience.compute_indicator(patches, False) def productivity(self, patches): return self._productivity.compute_indicator(patches) def biodiversity(self, patches): return self._biodiversity.compute_indicator(patches) def social(self, patches): return self._social.compute_indicator(patches) class Scenario : ''' This class is used to create a new randomized scenario from an initial scenario. A few functions are included to permit the manipulation of this object. The particularity of this class is the way he permit to chose a direction for the random. You have to chose the proportion of each kind of cultural patches and the algorithme will generate a new scenario who will try to approach this proportion. ''' def return_culture_code(self, culture): ''' Return the code associated to the kind of cultural type put as a parameter of the function :param culture: Take the value of a kind of cultural type (Cereales, Culture industrielles, Fruits et legumes Oleagineux, Prairies, Proteagineux) ''' if culture == "Céréales": #print culture return 0 if culture == "Cultures industrielles" : return 1 if culture == "Fourrages" : return 2 if culture == "Fruits et légumes" : #print culture return 3 if culture == "Oléagineux" : #print culture return 4 if culture == "Prairies" : return 5 if culture == "Protéagineux" : return 6 def __init__(self, indicators, scenarioInitial, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste, rng) : ''' Create a new randomized scenario according to the parameters :param scenarioInitial: the new scenario will originate from this one, he will take his values then change them :param nbr_parcelles_a_modifier: the number of patches who will be changed from the initial scenario :param cereales, culturesI, fruits_legumes, oleagineux, prairies, proteagineux: the proportion of area of each kind of cultural patches we want to approach in the new scenario :param liste: if 'ScenarioInitial' is a filename, chose False and if it's a list of patches, chose True. It's helpful to manage both kind of initial scenario. ''' self.shape = scenarioInitial self.nbr_parcelles_a_modifier = nbr_parcelles_a_modifier self.cereales = cereales self.culturesI = culturesI self.fourrages = fourrages self.fruits_legumes = fruits_legumes self.oleagineux = oleagineux self.prairies = prairies self.proteagineux = proteagineux self.indicators = indicators self.initial_patches = gpd.GeoDataFrame.from_file(scenarioInitial, encoding='utf-8') # If ScenarioInitial is a filename -> transform it in a list if liste == False : scenario = [] with fiona.open(self.shape, 'r', encoding='utf-8') as input: for f in input : scenario.append(f) self.shape = scenario else: scenario = scenarioInitial # Sort the initial list because the work done later will change the order. This way, we will be able to order the final list too et not be confused when we use the FID. newList = sorted(scenario, key=lambda k: k['properties']['ID_PARCEL']) self.shape = newList input = self.shape resultat_final = [] nombre_culture_a_modifie = self.nbr_parcelles_a_modifier #Area wanted for each kind of cultural patches for the simulation # 0 : Céréales, 1 : culture industrielle, 2 : fourrages, 3 : fruits et légumes, 4 : oléagineux, 5 : prairies, 6 : protéagineux pourcentage_surface = {} pourcentage_surface[0] = self.cereales #Céréales pourcentage_surface[1] = self.culturesI # Cultures industrielles pourcentage_surface[2] = self.fourrages # Fourrages pourcentage_surface[3] = self.fruits_legumes # Fruits et légumes pourcentage_surface[4] = self.oleagineux # Oléagineux pourcentage_surface[5] = self.prairies # Prairies pourcentage_surface[6] = self.proteagineux # Protéagineux #Initialize the area of each kind of cultural patches to 0. They will be incremented accordingly of the state of the scenario cultures = {} for i in range (8) : cultures[i] = { 'surf': 0.0, 'parcelles': []} #Total area of PAT patches total_culture = 0.0 #153256 i = 0 for f in input : total_culture += float(f["properties"]["SURF_PARC"]) # Initialize the area of the initial scenario for each kind of culture and store the FID of patches accordingly of the culture of the patch culture = f["properties"]["cultgeopat"] if culture != None and len(culture)>0: cultures[self.return_culture_code(culture)]['surf'] += float(f["properties"]["SURF_PARC"]) cultures[self.return_culture_code(culture)]['parcelles'].append(i) else : cultures[7]['parcelles'].append(i) cultures[7]['surf'] += float(f["properties"]["SURF_PARC"]) #Increment the FID of patches i += 1 liste_cultures_a_reduire = [] liste_cultures_a_augmenter = [] for i in range(7) : if cultures[i]['surf']*100/total_culture <= pourcentage_surface[i] : liste_cultures_a_augmenter.append(i) else : liste_cultures_a_reduire.append(i) culture_tab = ["Céréales", "Cultures industrielles", "Fourrages", "Fruits et légumes", "Oléagineux", "Prairies", "Protéagineux", None] cpt = 0 while cpt < nombre_culture_a_modifie and len(liste_cultures_a_augmenter) + len(liste_cultures_a_reduire) > 1 : #Get a random patches who need to reduce in area if len(liste_cultures_a_reduire) != 0 : rand_cul = rng.choice(liste_cultures_a_reduire) else : # Case when the desired area is reached in all kind of cultural patches but there is still patches to modify # We add a random culture to increase in the list of culture we need to reduce, so the algorithme go on until he reach "nombre_culture_a_modifie" rand_cul = rng.choice(liste_cultures_a_augmenter) liste_cultures_a_augmenter.remove(rand_cul) liste_cultures_a_reduire.append(rand_cul) if len(liste_cultures_a_augmenter) == 0 : rand_augmenter = rng.choice(liste_cultures_a_reduire) liste_cultures_a_augmenter.append(rand_augmenter) liste_cultures_a_reduire.remove(rand_augmenter) random_parcelles = rng.choice(cultures[rand_cul]['parcelles']) code_culture_rand = input[random_parcelles]["properties"]["cultgeopat"] recW = {} recW["geometry"]=input[random_parcelles]["geometry"] recW["properties"] = {} recW["properties"]['ID_PARCEL']=input[random_parcelles]["properties"]["ID_PARCEL"] recW["properties"]['CODE_CULTU']=input[random_parcelles]["properties"]["CODE_CULTU"] recW["properties"]['SURF_PARC']=input[random_parcelles]["properties"]["SURF_PARC"] recW["properties"]['Bdv']=input[random_parcelles]["properties"]["Bdv"] recW["properties"]['INSEE_COM']=input[random_parcelles]["properties"]["INSEE_COM"] recW["properties"]['VALEUR_CAD']=input[random_parcelles]["properties"]["VALEUR_CAD"] succes = 0 while succes == 0 : rand_cult = rng.choice(liste_cultures_a_augmenter) cultures[rand_cult]['surf'] += float(input[random_parcelles]["properties"]["SURF_PARC"] ) recW["properties"]['cultgeopat'] = culture_tab[rand_cult] #Remove the patch 'random_parcelles' of 'liste_nombre' cultures[self.return_culture_code(input[random_parcelles]["properties"]["cultgeopat"])]['parcelles'].remove(random_parcelles) #Remove the area of the old value contained in 'culture[type_culture] cultures[self.return_culture_code(input[random_parcelles]["properties"]["cultgeopat"])]['surf'] -= float(input[random_parcelles]["properties"]["SURF_PARC"]) resultat_final.append(recW) # If the proportion of the area of one kind of culture exceeds the proportion wanted # We remove it from the list of culture to increase and we add in the one list of the one to reduce 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) # Same than last block but in reverse, reduce_list to increase_list if float(cultures[self.return_culture_code(input[random_parcelles]["properties"]["cultgeopat"])]['surf']*100/total_culture) <= pourcentage_surface[self.return_culture_code(input[random_parcelles]["properties"]["cultgeopat"])] : liste_cultures_a_reduire.remove(self.return_culture_code(input[random_parcelles]["properties"]["cultgeopat"])) liste_cultures_a_augmenter.append(self.return_culture_code(input[random_parcelles]["properties"]["cultgeopat"])) succes = 1 cpt += 1 # Remove the comment to see the proportion change during the algorithm #for i in range(7) : #print cultures[i]['surf']*100/total_culture # If all patches of a kind of culture as already been modified, we remove this culture type from the list (increase_list or reduce_list) for i in range(7): if len(cultures[i]["parcelles"]) == 0 : if i in liste_cultures_a_reduire : liste_cultures_a_reduire.remove(i) elif i in liste_cultures_a_augmenter : liste_cultures_a_augmenter.remove(i) # Every patch remaining will be added in the new layer. Patches with an empty type will be added now for i in range(8) : for e in cultures[i]["parcelles"] : recW = {} recW["geometry"] = input[e]["geometry"] recW["properties"] = {} recW["properties"]['ID_PARCEL'] = input[e]["properties"]["ID_PARCEL"] recW["properties"]['CODE_CULTU'] = input[e]["properties"]["CODE_CULTU"] recW["properties"]['SURF_PARC'] = input[e]["properties"]["SURF_PARC"] recW["properties"]['Bdv'] = input[e]["properties"]["Bdv"] recW["properties"]['id_expl'] = input[e]["properties"]["id_expl"] recW["properties"]['INSEE_COM'] = input[e]["properties"]["INSEE_COM"] recW["properties"]['VALEUR_CAD'] = input[e]["properties"]["VALEUR_CAD"] recW["properties"]['cultgeopat'] = culture_tab[i] resultat_final.append(recW) #print #print ("Pourcentage culture apres MCMC : ") #for i in range(7) : #print (cultures[i]['surf']*100/total_culture) # Sort the final list to make it ordered as the same way of the initial list resultat_final = sorted(resultat_final, key=lambda k: k['properties']['ID_PARCEL']) self.resultat_final = resultat_final self.resultat_final_gpd = gpd.GeoDataFrame.from_features(resultat_final) def get_scenario(self): ''' Return the object himself ''' return self def get_resultat_final(self): ''' Return all entities of the scenario in the list format (faster to work on than a shapefile) ''' return self.resultat_final def write_in_shape(self, shape) : ''' Write the scenario in a shapefile, permits to visualize the map :param shape: the name of the file which will contains the scenario ''' schema = {'geometry': 'Polygon', 'properties': {'CODE_CULTU' : 'str', 'Bdv' : 'str', 'ID_PARCEL' : 'str', 'INSEE_COM' : 'str', 'SURF_PARC' : 'float', 'cultgeopat' : 'str'}} crs = {'lon_0': 3, 'ellps': 'GRS80', 'y_0': 6600000, 'no_defs': True, 'proj': 'lcc', 'x_0': 700000, 'units': 'm', 'lat_2': 44, 'lat_1': 49, 'lat_0': 46.5} driver = "ESRI Shapefile" with fiona.open(shape,'w',driver=driver, crs=crs, schema=schema, encoding='utf-8') as parcelles : for f in self.resultat_final : parcelles.write(f) def write_in_json(self): ''' *** not used, can be removed/modified *** Write the scenario in a Json file ''' data = self.resultat_final with gzip.GzipFile('./best_scenario.gz','w') as fid_gz: json_str = json.dumps(self.resultat_final) fid_gz.write(json_str) def write_indices_in_json(self): ''' Write the values of indices of the initial scenario in a JSON file This function is used for display and statistics in the application ''' dict_indices = {} dict_indices["proximite"] = self.indice_proximite() dict_indices["resilience"] = self.indice_resilience() dict_indices["productivite"] = self.indice_productivitee() dict_indices["biodiversite"] = self.indice_biodiversite() dict_indices["social"] = self.indice_social() with gzip.GzipFile('./initial_scenario_indices.gz','w') as fid_gz: json_str = json.dumps(dict_indices, ensure_ascii=False) fid_gz.write(json_str.encode('utf-8')) ''' These functions returns values of indices of this scenario ''' def indice_proximite(self) : return self.indicators.proximity(self.resultat_final_gpd) def indice_resilience(self): return self.indicators.resilience(self.resultat_final_gpd) def indice_productivitee(self): return self.indicators.productivity(self.resultat_final_gpd) def indice_biodiversite(self): return self.indicators.biodiversity(self.resultat_final_gpd) def indice_social(self): return self.indicators.social(self.resultat_final_gpd) dict_scenario = {} # Contains all scenarios created by the function "générer_x_scenarios" ''' These few function are used to get some data and make stats. They are not called in the main program but as they can be useful, we keep them. ''' def generer_x_scenarios(x, scenario_initial, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste ): #Génère x scénarios et les ranges dans "dict_scenario" ''' Generate X scenarios and stocks them in a dict. It will be used to compare all scenario in the dict and keep the betters one ''' dict_scenario = {} for i in range(x): s = Scenario(scenario_initial, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste) s_r = s.get_scenario() dict_scenario[i] = s_r #comparer_2_scenarios(dict_scenario[0], dict_scenario[1]) return dict_scenario def x_all_indices(scenarInitial, x, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux): ''' Create X numbers of scenario, each new scenario is created from an initial one. ''' scenar = scenarInitial liste = False for i in range(x) : s = Scenario(scenar, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste) s.indice_proximite() s.indice_resilience() s.indice_productivitee() s.indice_biodiversite() scenar = s.get_resultat_final() liste = True def x_proximite(scenarInitial, x, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux): ''' Calculate X number of proximity indices Create an X numbers of scenario, each new scenario is created from an initial one. It will create 4 set of data and put them into a file readable as a csv. Each set of data change a different number of patches. ''' scenar = scenarInitial liste = False fichier = open("./indice_proximite_300.txt", "w") for i in range(x) : print(i) list = [300, 1000, 5000, 10000] for nbr in list : nbr_parcelles_a_modifier = nbr s = Scenario(scenar, nbr, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste) ind = s.indice_proximite() if nbr != 10000: fichier.write(str(ind) + ";") else : fichier.write(str(ind) + "\n") def x_resilience(scenarInitial, x, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux): ''' Calculate X number of resiliency indices Create an X numbers of scenario, each new scenario is created from an initial one. It will create 4 set of data and put them into a file readable as a csv. Each set of data change a different number of patches. ''' scenar = scenarInitial liste = False fichier = open("./indice_resilience.txt", "w") for i in range(x) : print(i) list = [300, 1000, 5000, 10000] for nbr in list : nbr_parcelles_a_modifier = nbr s = Scenario(scenar, nbr, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste) #print "Scenario " + str(i) + " : " ind = s.indice_resilience() #scenar = s.get_resultat_final() #liste = True if nbr != 10000: fichier.write(str(ind) + ";") else : fichier.write(str(ind) + "\n") def x_productivite(scenarInitial, x, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux): ''' Calculate X number of productivity indices Create an X numbers of scenario, each new scenario is created from an initial one. It will create 4 set of data and put them into a file readable as a csv. Each set of data change a different number of patches. ''' scenar = scenarInitial liste = False fichier = open("./indice_productivite.txt", "w") for i in range(x) : print(i) list = [300, 1000, 5000, 10000] for nbr in list : s = Scenario(scenar, nbr, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste) #print "Scenario " + str(i) + " : " ind = s.indice_productivitee() #scenar = s.get_resultat_final() #liste = True if nbr != 10000: fichier.write(str(ind) + ";") else : fichier.write(str(ind) + "\n") def x_biodiversite(scenarInitial, x, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux): ''' Calculate X number of biodiversity indices Create an X numbers of scenario, each new scenario is created from an initial one. It will create 4 set of data and put them into a file readable as a csv. Each set of data change a different number of patches. ''' scenar = scenarInitial liste = False fichier = open("./indice_biodiversite.txt", "w") for i in range(x) : print(i) list = [300, 1000, 5000, 10000] for nbr in list : s = Scenario(scenar, nbr, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste) #print "Scenario " + str(i) + " : " ind = s.indice_biodiversite() #scenar = s.get_resultat_final() #liste = True if nbr != 10000: fichier.write(str(ind) + ";") else : fichier.write(str(ind) + "\n") def x_social(scenarInitial, x, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux): ''' Calculate X number of social indices Create an X numbers of scenario, each new scenario is created from an initial one. It will create 4 set of data and put them into a file readable as a csv. Each set of data change a different number of patches. ''' scenar = scenarInitial liste = False fichier = open("./indice_social.txt", "w") for i in range(x) : print(i) list = [300, 1000, 5000, 10000] for nbr in list : s = Scenario(scenar, nbr, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste) #print "Scenario " + str(i) + " : " ind = s1.indice_social() #scenar = s.get_resultat_final() #liste = True if nbr != 10000: fichier.write(str(ind) + ";") else : fichier.write(str(ind) + "\n") def comparer_2_scenarios(scenar1, scenar2) : ''' *** Not really useful anymore *** Compare 2 scenarios who are in the global variable 'dict_scenario' ''' dico = {} dico[0] = {} dico[1] = {} s1_proximite = scenar1.indice_proximite() dico[0]["proximite"] = s1_proximite s2_proximite = scenar2.indice_proximite() dico[1]["proximite"] = s2_proximite s1_resilience = scenar1.indice_resilience() dico[0]["resilience"] = s1_resilience s2_resilience = scenar2.indice_resilience() dico[1]["resilience"] = s2_resilience s1_productivite = scenar1.indice_productivitee() dico[0]["productivite"] = s1_productivite s2_productivite = scenar2.indice_productivitee() dico[1]["productivite"] = s2_productivite s1_biodiversite = scenar1.indice_biodiversite() dico[0]["biodiversite"] = s1_biodiversite s2_biodiversite = scenar2.indice_biodiversite() dico[1]["biodiversite"] = s2_biodiversite dico[0]["social"] = 0 dico[1]["social"] = 0 print() if s1_proximite > s2_proximite : print("Le scénario 1 possède un meilleur indice de proximité que le scenario 2 de : " + str((s1_proximite - s2_proximite) / s1_proximite * 100) + "%") else : print("Le scénario 2 possède un meilleur indice de proximité que le scénario 1 de : " + str((s2_proximite - s1_proximite) / s2_proximite * 100) + "%") print() if s1_resilience > s2_resilience : print("Le scénario 1 possède un meilleur indice de résilience que le scenario 2 de : " + str((s1_resilience - s2_resilience) / s1_resilience * 100) + "%") else : print("Le scénario 2 possède un meilleur indice de résilience que le scénario 1 de : " + str((s2_resilience - s1_resilience) / s2_resilience * 100) + "%") print() if s1_productivite > s2_productivite : print("Le scénario 1 possède un meilleur indice de productivité que le scenario 2 de : " + str((s1_productivite - s2_productivite) / s1_productivite * 100) + "%") else : print("Le scénario 2 possède un meilleur indice de productivité que le scénario 1 de : " + str((s2_productivite - s1_productivite) / s2_productivite * 100) + "%") if s1_biodiversite > s2_biodiversite : print("Le scénario 1 possède un meilleur indice de biodiversite que le scenario 2 de : " + str((s1_biodiversite - s2_biodiversite) / s1_biodiversite * 100) + "%") else : print("Le scénario 2 possède un meilleur indice de biodiversite que le scénario 1 de : " + str((s2_biodiversite - s1_biodiversite) / s2_biodiversite * 100) + "%") scenarios_graph(2, dico) #Use to stock in JSON file the data about the area of each EPCI for each best scenario def EPCI_repartition(best_scenario): ''' Function used to get the area of each kind of culture for each EPCI in the PAT territory :param best_scenario: dictionnary of the final best scenario produced by the MCMC function ''' list_EPCI_scenario = [] EPCI_layer = [] with fiona.open("./EPCI_clean.shp", 'r', encoding='utf-8') as epci : for f in epci : EPCI_layer.append(f) for i in range(len(best_scenario)) : print(best_scenario[i]['proximite']) list_EPCI_scenario.append(EPCI_culture_repartition(EPCI_layer, best_scenario[i]['scenario'])) with gzip.GzipFile('./EPCI_repartition.gz','w') as fid_gz: json_str = json.dumps(list_EPCI_scenario) fid_gz.write(json_str) def get_best_scenario(dict_scenar, nbr_scenar_to_stock): ''' Function used is the MCMC function Get the 'nbr_scenar_to_stock' best scenario for each kind of cultural patches from a dictionary 'dict_scenar' :param dict_scenar: dictionary containing a definite by the MCMC function numbers of scenarios :param nbr_scenar_to_stock: number of scenario we want to keep for each indices ''' dict_resultat = {} list_indices = ['proximite', 'resilience', 'productivite', 'biodiversite', 'social'] for i in range(5): dict_resultat[i] = {} dict_resultat[i]['values'] = [] dict_resultat[i]['scenario'] = [] # Add nbr_scenar_to_stock number of scenarios in the dictionnary to start the sort of the best scenarios for y in range(nbr_scenar_to_stock) : for indices in range(5): dict_resultat[indices]['scenario'].append(dict_scenar[y]) dict_resultat[indices]['values'].append(dict_scenar[y][list_indices[indices]]) # Loop who select the best scenarios for i in range(len(dict_scenar)): for y in range(nbr_scenar_to_stock) : for indices in range(5): # Low = best : proximite high = best : resilience, productivite, biodiversite, social # For the proximity indice, the lowest value is the best, a special loop is created for this indice if i == 0 : if dict_scenar[i][list_indices[indices]] < dict_resultat[indices]['values'][y] and not(dict_scenar[i][list_indices[indices]] in dict_resultat[indices]['values']) : del dict_resultat[indices]['values'][y] dict_resultat[indices]['values'].append(dict_scenar[i][list_indices[indices]]) del dict_resultat[indices]['scenario'][y] dict_resultat[indices]['scenario'].append(dict_scenar[i]) elif dict_scenar[i][list_indices[indices]] > dict_resultat[indices]['values'][y] and not(dict_scenar[i][list_indices[indices]] in dict_resultat[indices]['values']) : del dict_resultat[indices]['values'][y] dict_resultat[indices]['values'].append(dict_scenar[i][list_indices[indices]]) del dict_resultat[indices]['scenario'][y] dict_resultat[indices]['scenario'].append(dict_scenar[i]) return dict_resultat dict_scenario = {} dump_scenario = {} dump_cpt = 0 def MCMC(scenarInitial, nbr_loop, nbr_scenario, nbr_scenario_to_stock, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, rng): ''' Final function of the project. His purpose is to find X best scenario from a Y number of generated scenario multiplied by the number of loops we want. The first loop will create scenarios from the initial scenario put in the parameters. The others loops will used the already chosen best scenario to create new ones. This method will permit to analyse a high variety of different scenario and enhance the size of the spectrum of data created. The function save in an Sqlite database all the scenarios created to keep trace of them. Database name : scenario_dump Table name : scenario :param scenarInitial: Initial scenario who will be used to create the scenarios in the first loop :param nbr_loop: The number of loop who will run in this algorithme. The first loop is counted as one, so the value need to be at least 1 :param nbr_scenario: The number of scenario created in each loop :param nbr_scenario_to_stock: The number of best scenario saved in a dictionnary at the end of each loop for each indices. The final number of best scenario saved is 'nbr_scenario_to_stock' * number of indices (5) :param nbr_parcelles_a_modifier: the number of patches who will be changed from the initial scenario :param cereales, culturesI, fruits_legumes, oleagineux, prairies, proteagineux: the proportion of area of each kind of cultural patches we want to approach in the new scenario :param liste: if 'ScenarioInitial' is a filename, chose False and if it's a list of patches, chose True. It's helpful to manage both kind of initial scenario. ''' liste = True global dump_cpt # Put in a sqlite database all the scenario created during the algorithme conn = sqlite3.connect('./scenario_dump.db') cursor = conn.cursor() cursor.execute(""" DROP TABLE IF EXISTS scenario """) conn.commit() cursor.execute(""" CREATE TABLE scenario ( ID integer primary key, scenario VARCHAR(10000000), indice_proximite double, indice_resilience double, indice_productivite double, indice_biodiversite double, indice_social double ) """) conn.commit() scenario = [] with fiona.open(scenarInitial, 'r', encoding='utf-8') as input: for f in input : scenario.append(f) # Get the values of indices of the initial scenario si = Scenario(scenario, 0, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, True, rng) si_indices = {} si_indices['proximite'] = si.indice_proximite() si_indices['resilience'] = si.indice_resilience() si_indices['productivite'] = si.indice_productivitee() si_indices['biodiversite'] = si.indice_biodiversite() si_indices['social'] = si.indice_social() s_new_indices_avg = {} s_new_indices_avg['proximite'] = 0 s_new_indices_avg['resilience'] = 0 s_new_indices_avg['productivite'] = 0 s_new_indices_avg['biodiversite'] = 0 s_new_indices_avg['social'] = 0 for nbr_scenar in range(nbr_scenario): s = Scenario(scenario, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste, rng) dict_scenario[nbr_scenar] = {} dict_scenario[nbr_scenar]['scenario'] = s.get_scenario() dict_scenario[nbr_scenar]['proximite'] = s.indice_proximite() dict_scenario[nbr_scenar]['resilience'] = s.indice_resilience() dict_scenario[nbr_scenar]['productivite'] = s.indice_productivitee() dict_scenario[nbr_scenar]['biodiversite'] = s.indice_biodiversite() dict_scenario[nbr_scenar]['social'] = s.indice_social() s_new_indices_avg['proximite'] += dict_scenario[nbr_scenar]['proximite'] s_new_indices_avg['resilience'] += dict_scenario[nbr_scenar]['resilience'] s_new_indices_avg['productivite'] += dict_scenario[nbr_scenar]['productivite'] s_new_indices_avg['biodiversite'] += dict_scenario[nbr_scenar]['biodiversite'] s_new_indices_avg['social'] += dict_scenario[nbr_scenar]['social'] #dump_scenario[dump_cpt] = s.get_resultat_final() cursor.execute(""" INSERT INTO scenario(id, scenario, indice_proximite, indice_resilience, indice_productivite, indice_biodiversite, indice_social) VALUES(?, ?, ?, ?, ?, ?, ?) """, (dump_cpt ,str(s.get_resultat_final()), dict_scenario[nbr_scenar]['proximite'], dict_scenario[nbr_scenar]['resilience'], dict_scenario[nbr_scenar]['productivite'], dict_scenario[nbr_scenar]['biodiversite'], dict_scenario[nbr_scenar]['social']) ) dump_cpt += 1 conn.commit() list_indices = ['proximite', 'resilience', 'productivite', 'biodiversite', 'social'] ''' for i in list_indices : print i +" initiale : " + str(si_indices[i]) print i + " average of new scenarios : " + str(s_new_indices_avg[i] / nbr_scenario) ''' value_end_loop = {} for i in list_indices: value_end_loop[i] = abs(si_indices[i] - (s_new_indices_avg[i] / nbr_scenario)) / 100 best_scenario = get_best_scenario(dict_scenario, nbr_scenario_to_stock) #Triggers of end condition trigger = {} trigger['proximite'] = 0 trigger['resilience'] = 0 trigger['productivite'] = 0 trigger['biodiversite'] = 0 trigger['social'] = 0 cpt_loop = 0 fichier = open("./difference_values.txt", "w") list_indices = ['proximite', 'resilience', 'productivite', 'biodiversite', 'social'] while cpt_loop < nbr_loop-1 and (trigger['proximite'] + trigger['resilience'] + trigger['productivite'] + trigger['biodiversite'] + trigger['social']) < 5 : trigger['proximite'] = 0 trigger['resilience'] = 0 trigger['productivite'] = 0 trigger['biodiversite'] = 0 trigger['social'] = 0 cpt = 0 for i in best_scenario: for val in best_scenario[i]['scenario']: dict_scenario[cpt] = val for nbr_scenar in range(nbr_scenario - nbr_scenario_to_stock*5, nbr_scenario) : # Create new scenario using the best scenario chosen in the previous step s = Scenario(dict_scenario[nbr_scenar%(nbr_scenario_to_stock * 5)]['scenario'].get_resultat_final(), nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste, rng) dict_scenario[nbr_scenar]['scenario'] = s.get_scenario() dict_scenario[nbr_scenar]['proximite'] = s.indice_proximite() dict_scenario[nbr_scenar]['resilience'] = s.indice_resilience() dict_scenario[nbr_scenar]['productivite'] = s.indice_productivitee() dict_scenario[nbr_scenar]['biodiversite'] = s.indice_biodiversite() dict_scenario[nbr_scenar]['social'] = s.indice_social() old_best_scenario = best_scenario best_scenario = get_best_scenario(dict_scenario, nbr_scenario_to_stock) cpt_loop += 1 for ind in range(5): sum_old = 0 sum_new = 0 for y in range(nbr_scenario_to_stock) : sum_old += old_best_scenario[ind]['values'][y] sum_new += best_scenario[ind]['values'][y] difference = abs(((sum_old/nbr_scenario_to_stock) - (sum_new/nbr_scenario_to_stock)) / (sum_new/nbr_scenario_to_stock) * 100) fichier.write(str(difference) + ";") #print "Difference for old and new " + list_indices[ind] + " is : " + str(difference) #print si_indices[list_indices[ind]] #print value_end_loop[list_indices[ind]] #print if difference <= value_end_loop[list_indices[ind]] : trigger[list_indices[ind]] = 1 scenario_json = {} for i in range(nbr_scenario_to_stock*5): scenario_json[i] = {} cpt_json = 0 for i in range(len(best_scenario)): for y in range(nbr_scenario_to_stock): scenario_json[cpt_json]['scenario'] = best_scenario[i]['scenario'][y]['scenario'].get_resultat_final() scenario_json[cpt_json]['proximite'] = best_scenario[i]['scenario'][y]['proximite'] scenario_json[cpt_json]['resilience'] = best_scenario[i]['scenario'][y]['resilience'] scenario_json[cpt_json]['productivite'] = best_scenario[i]['scenario'][y]['productivite'] scenario_json[cpt_json]['biodiversite'] = best_scenario[i]['scenario'][y]['biodiversite'] scenario_json[cpt_json]['social'] = best_scenario[i]['scenario'][y]['social'] cpt_json += 1 EPCI_repartition(scenario_json) with gzip.GzipFile('./best_scenario.gz','w') as fid_gz: json_str = json.dumps(scenario_json) fid_gz.write(json_str) if __name__ == '__main__': # Modifie these values to define how many of each cultures you want in your scenario, the modifications will tends to approach these values cereales = 22 culturesI = 5 fourrages = 5 fruits_legumes = 3 oleagineux = 1 prairies = 59 proteagineux = 2 # Change to select how many patches you will change from the initial scenario nbr_parcelles_a_modifier = 20000 # Random seed seed = random.randrange(sys.maxsize) rng = random.Random(seed) print(("Seed was:", seed)) fichier_seed = open("./lastest_seed.txt", "w") fichier_seed.write(str(seed)) #s1 = Scenario("./Parcelle_PAT_MCMC.shp",nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, False) #s1.indice_proximite() #s1.indice_resilience('../Parcelle_PAT/Grille_resilience.shp') #s1.indice_productivitee("../Parcelle_PAT/Parcelle_PAT_valCad.shp") #s1.indice_biodiversite('./compute_dispersal_probabilities_1000.npz') #s1.indice_biodiversite('./matrix_biodiversity.npz') #s1.indice_social('./Bassin_de_vie_Pat.shp', './Parcelle_PAT_exploitant.shp') #s1.write_in_json() #s1.write_in_shape('./testSimulationSurf_MCMC.shp') #dict_scenario = generer_x_scenarios(2, "./Parcelle_PAT_MCMC.shp", nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, False) #x_proximite("./Parcelle_PAT_MCMC.shp", 100 , nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux) #x_resilience("./Parcelle_PAT_MCMC.shp", 100 , nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux) #x_productivite("./Parcelle_PAT_MCMC.shp", 100 , nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux) #x_biodiversite("./Parcelle_PAT_MCMC.shp", 1 , nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux) #x_social("./Parcelle_PAT_MCMC.shp", 1 , nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux) # scenarInitial, nbr_loop, nbr_scenario, nbr_scenario_to_stock, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux #scenario_initial = Scenario("./Parcelle_PAT_MCMC.shp",0, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, False, rng) patches = gpd.GeoDataFrame.from_file("../output/PAT_patches/PAT_patches.shp", encoding='utf-8') indicators = Indicators('MCMC_config.yml', patches) scenario_initial = Scenario(indicators,"../output/PAT_patches/PAT_patches.shp",0, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, False, rng) #scenario_initial = Scenario("../output/patches_RiomLimagneVolcans.shp",0, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, False, rng) scenario_initial.write_indices_in_json() #MCMC("./Parcelle_PAT_MCMC.shp", 2, 11, 2, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, rng)