Commit 96195249 authored by Dumoulin Nicolas's avatar Dumoulin Nicolas
Browse files

fixes for data produced by PATBuilder

optimization fixes
parent 9fd92c31
Céréales Cultures industrielles Fourrages Fruits et légumes Oléagineux Prairies Protéagineux
Céréales 0.0 0.5 0.5 1 0.5 0 0.5
Cultures industrielles 0.5 0.0 0.5 1 0.5 0 0.5
Fourrages 0.5 0.5 0.0 1 0.5 0 0.5
Fruits et légumes 0.5 0.5 0.5 0 0.5 0 0.5
Oléagineux 0.5 0.5 0.5 1 0.0 0 0.5
Prairies 0.5 0.5 0.5 1 0.5 0 0.5
Protéagineux 0.5 0.5 0.5 1 0.5 0 0.0
...@@ -16,12 +16,12 @@ import sqlite3 ...@@ -16,12 +16,12 @@ import sqlite3
import sys, os import sys, os
import yaml import yaml
from proximite import proximite from proximite import Proximite
from resilience_list import resilience from resilience_list import Resilience
from productivite import indice_productivite from productivite import Productivity
#from indice_biodiversite import biodiversite #from indice_biodiversite import biodiversite
from indice_biodiversite_2 import biodiversity from indice_biodiversite_2 import Biodiversity
from social import social from social import Social
from EPCI_culture_repartition import EPCI_culture_repartition from EPCI_culture_repartition import EPCI_culture_repartition
from radar import scenarios_graph from radar import scenarios_graph
...@@ -42,7 +42,7 @@ class Scenario : ...@@ -42,7 +42,7 @@ class Scenario :
:param culture: Take the value of a kind of cultural type (Cereales, Culture industrielles, Fruits et legumes :param culture: Take the value of a kind of cultural type (Cereales, Culture industrielles, Fruits et legumes
Oleagineux, Prairies, Proteagineux) Oleagineux, Prairies, Proteagineux)
''' '''
if culture == "Cereales": if culture == "Céréales":
#print culture #print culture
return 0 return 0
...@@ -52,18 +52,18 @@ class Scenario : ...@@ -52,18 +52,18 @@ class Scenario :
if culture == "Fourrages" : if culture == "Fourrages" :
return 2 return 2
if culture == "Fruits et legumes" : if culture == "Fruits et légumes" :
#print culture #print culture
return 3 return 3
if culture == "Oleagineux" : if culture == "Oléagineux" :
#print culture #print culture
return 4 return 4
if culture == "Prairies" : if culture == "Prairies" :
return 5 return 5
if culture == "Proteagineux" : if culture == "Protéagineux" :
return 6 return 6
def __init__(self, scenarioInitial, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste, rng) : def __init__(self, scenarioInitial, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, liste, rng) :
...@@ -93,10 +93,13 @@ class Scenario : ...@@ -93,10 +93,13 @@ class Scenario :
sys.exit(1) sys.exit(1)
self.mcmc_config = yaml.load(open(mcmc_config_filename,'r')) self.mcmc_config = yaml.load(open(mcmc_config_filename,'r'))
self.initial_patches = gpd.GeoDataFrame.from_file(scenarioInitial, encoding='utf-8')
self.proximite = Proximite(self.initial_patches)
# If ScenarioInitial is a filename -> transform it in a list # If ScenarioInitial is a filename -> transform it in a list
if liste == False : if liste == False :
scenario = [] scenario = []
with fiona.open(self.shape, 'r') as input: with fiona.open(self.shape, 'r', encoding='utf-8') as input:
for f in input : for f in input :
scenario.append(f) scenario.append(f)
self.shape = scenario self.shape = scenario
...@@ -137,7 +140,7 @@ class Scenario : ...@@ -137,7 +140,7 @@ class Scenario :
total_culture += float(f["properties"]["SURF_PARC"]) 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 # 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"]["rpg_lfgc_s"] culture = f["properties"]["cultgeopat"]
if culture != None and len(culture)>0: 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)]['surf'] += float(f["properties"]["SURF_PARC"])
cultures[self.return_culture_code(culture)]['parcelles'].append(i) cultures[self.return_culture_code(culture)]['parcelles'].append(i)
...@@ -156,9 +159,7 @@ class Scenario : ...@@ -156,9 +159,7 @@ class Scenario :
else : else :
liste_cultures_a_reduire.append(i) liste_cultures_a_reduire.append(i)
culture_tab = ["Cereales", "Cultures industrielles", "Fourrages", "Fruits et legumes", "Oleagineux", "Prairies", "Proteagineux", None] culture_tab = ["Céréales", "Cultures industrielles", "Fourrages", "Fruits et légumes", "Oléagineux", "Prairies", "Protéagineux", None]
cpt = 0 cpt = 0
while cpt < nombre_culture_a_modifie and len(liste_cultures_a_augmenter) + len(liste_cultures_a_reduire) > 1 : while cpt < nombre_culture_a_modifie and len(liste_cultures_a_augmenter) + len(liste_cultures_a_reduire) > 1 :
...@@ -180,7 +181,7 @@ class Scenario : ...@@ -180,7 +181,7 @@ class Scenario :
random_parcelles = rng.choice(cultures[rand_cul]['parcelles']) random_parcelles = rng.choice(cultures[rand_cul]['parcelles'])
code_culture_rand = input[random_parcelles]["properties"]["rpg_lfgc_s"] code_culture_rand = input[random_parcelles]["properties"]["cultgeopat"]
recW = {} recW = {}
recW["geometry"]=input[random_parcelles]["geometry"] recW["geometry"]=input[random_parcelles]["geometry"]
...@@ -190,19 +191,20 @@ class Scenario : ...@@ -190,19 +191,20 @@ class Scenario :
recW["properties"]['SURF_PARC']=input[random_parcelles]["properties"]["SURF_PARC"] recW["properties"]['SURF_PARC']=input[random_parcelles]["properties"]["SURF_PARC"]
recW["properties"]['Bdv']=input[random_parcelles]["properties"]["Bdv"] recW["properties"]['Bdv']=input[random_parcelles]["properties"]["Bdv"]
recW["properties"]['INSEE_COM']=input[random_parcelles]["properties"]["INSEE_COM"] recW["properties"]['INSEE_COM']=input[random_parcelles]["properties"]["INSEE_COM"]
recW["properties"]['VALEUR_CAD']=input[random_parcelles]["properties"]["VALEUR_CAD"]
succes = 0 succes = 0
while succes == 0 : while succes == 0 :
rand_cult = rng.choice(liste_cultures_a_augmenter) rand_cult = rng.choice(liste_cultures_a_augmenter)
cultures[rand_cult]['surf'] += float(input[random_parcelles]["properties"]["SURF_PARC"] ) cultures[rand_cult]['surf'] += float(input[random_parcelles]["properties"]["SURF_PARC"] )
recW["properties"]['rpg_lfgc_s'] = culture_tab[rand_cult] recW["properties"]['cultgeopat'] = culture_tab[rand_cult]
#Remove the patch 'random_parcelles' of 'liste_nombre' #Remove the patch 'random_parcelles' of 'liste_nombre'
cultures[self.return_culture_code(input[random_parcelles]["properties"]["rpg_lfgc_s"])]['parcelles'].remove(random_parcelles) 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] #Remove the area of the old value contained in 'culture[type_culture]
cultures[self.return_culture_code(input[random_parcelles]["properties"]["rpg_lfgc_s"])]['surf'] -= float(input[random_parcelles]["properties"]["SURF_PARC"]) cultures[self.return_culture_code(input[random_parcelles]["properties"]["cultgeopat"])]['surf'] -= float(input[random_parcelles]["properties"]["SURF_PARC"])
resultat_final.append(recW) resultat_final.append(recW)
# If the proportion of the area of one kind of culture exceeds the proportion wanted # If the proportion of the area of one kind of culture exceeds the proportion wanted
...@@ -212,9 +214,9 @@ class Scenario : ...@@ -212,9 +214,9 @@ class Scenario :
liste_cultures_a_reduire.append(rand_cult) liste_cultures_a_reduire.append(rand_cult)
# Same than last block but in reverse, reduce_list to increase_list # Same than last block but in reverse, reduce_list to increase_list
if float(cultures[self.return_culture_code(input[random_parcelles]["properties"]["rpg_lfgc_s"])]['surf']*100/total_culture) <= pourcentage_surface[self.return_culture_code(input[random_parcelles]["properties"]["rpg_lfgc_s"])] : 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"]["rpg_lfgc_s"])) 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"]["rpg_lfgc_s"])) liste_cultures_a_augmenter.append(self.return_culture_code(input[random_parcelles]["properties"]["cultgeopat"]))
succes = 1 succes = 1
cpt += 1 cpt += 1
...@@ -241,8 +243,10 @@ class Scenario : ...@@ -241,8 +243,10 @@ class Scenario :
recW["properties"]['CODE_CULTU'] = input[e]["properties"]["CODE_CULTU"] recW["properties"]['CODE_CULTU'] = input[e]["properties"]["CODE_CULTU"]
recW["properties"]['SURF_PARC'] = input[e]["properties"]["SURF_PARC"] recW["properties"]['SURF_PARC'] = input[e]["properties"]["SURF_PARC"]
recW["properties"]['Bdv'] = input[e]["properties"]["Bdv"] 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"]['INSEE_COM'] = input[e]["properties"]["INSEE_COM"]
recW["properties"]['rpg_lfgc_s'] = culture_tab[i] recW["properties"]['VALEUR_CAD'] = input[e]["properties"]["VALEUR_CAD"]
recW["properties"]['cultgeopat'] = culture_tab[i]
resultat_final.append(recW) resultat_final.append(recW)
...@@ -276,11 +280,11 @@ class Scenario : ...@@ -276,11 +280,11 @@ class Scenario :
Write the scenario in a shapefile, permits to visualize the map Write the scenario in a shapefile, permits to visualize the map
:param shape: the name of the file which will contains the scenario :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', 'rpg_lfgc_s' : 'str'}} 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} 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" driver = "ESRI Shapefile"
with fiona.open(shape,'w',driver=driver, crs=crs, schema=schema) as parcelles : with fiona.open(shape,'w',driver=driver, crs=crs, schema=schema, encoding='utf-8') as parcelles :
for f in self.resultat_final : for f in self.resultat_final :
parcelles.write(f) parcelles.write(f)
...@@ -317,37 +321,42 @@ class Scenario : ...@@ -317,37 +321,42 @@ class Scenario :
These functions returns values of indices of this scenario These functions returns values of indices of this scenario
''' '''
def indice_proximite(self) : def indice_proximite(self) :
indice_prox = proximite(self.resultat_final, self.fruits_legumes, False) indice_prox = self.proximite.compute_indicator(self.resultat_final, self.fruits_legumes, False)
#print "Résultats de l'indice de proximité : " + str(indice_prox) #print "Résultats de l'indice de proximité : " + str(indice_prox)
return indice_prox return indice_prox
def indice_resilience(self): def indice_resilience(self):
indice_resi = resilience(self.mcmc_config['indicators_config']['resilience'], self.resultat_final, False) resilience = Resilience(self.mcmc_config['indicators_config']['resilience'])
patches = gpd.GeoDataFrame.from_features(self.resultat_final)
indice_resi = resilience.compute_indicator(patches, False)
#print "Résultats de l'indice de résilience : " + str(indice_resi) #print "Résultats de l'indice de résilience : " + str(indice_resi)
return indice_resi return indice_resi
def indice_productivitee(self): def indice_productivitee(self):
indice_prod = indice_productivite(self.resultat_final, self.mcmc_config['indicators_config']['productivity'], False) productivity = Productivity()
patches = gpd.GeoDataFrame.from_features(self.resultat_final)
indice_prod = productivity.compute_indicator(patches)
#print "Résultats de l'indice de productivite : " + str(indice_prod) #print "Résultats de l'indice de productivite : " + str(indice_prod)
return indice_prod return indice_prod
def indice_biodiversite(self): def indice_biodiversite(self):
#Matrix containing all the data about biodiversity. Created with Numpy and the formula of Equivalent Connected Area (Saura et al., 2011) #Matrix containing all the data about biodiversity. Created with Numpy and the formula of Equivalent Connected Area (Saura et al., 2011)
matrice = self.mcmc_config['indicators_config']['biodiversity'] biodiv_config = self.mcmc_config['indicators_config']['biodiversity']
indice_bio = biodiversity(self.resultat_final, matrice) matrixfilename = biodiv_config['matrixfilename']
biodiversity = Biodiversity(biodiv_config['dmax'], biodiv_config['epsilon'])
biodiversity.load_probabilities(matrixfilename)
patches = gpd.GeoDataFrame.from_features(self.resultat_final)
indice_bio = biodiversity.compute_indicator(patches)
#print "Résultats de l'indice de biodiversite : " + str(indice_bio) #print "Résultats de l'indice de biodiversite : " + str(indice_bio)
return indice_bio return indice_bio
def indice_social(self): def indice_social(self):
bdv = self.mcmc_config['indicators_config']['social']['bdv'] social_config = self.mcmc_config['indicators_config']['social']
shape_exploitant = self.mcmc_config['indicators_config']['social']['exploitant'] social = Social(social_config['cost_matrix_filename'], social_config['bdv_threshold'])
indice_social = social(shape_exploitant, self.resultat_final, bdv) indice_social = social.compute_indicator(gpd.GeoDataFrame.from_features(self.shape), gpd.GeoDataFrame.from_features(self.resultat_final))
#print "Résultats de l'indice de sociabilité : " + str(indice_social) #print "Résultats de l'indice de sociabilité : " + str(indice_social)
return indice_social return indice_social
dict_scenario = {} # Contains all scenarios created by the function "générer_x_scenarios" dict_scenario = {} # Contains all scenarios created by the function "générer_x_scenarios"
''' '''
...@@ -575,7 +584,7 @@ def EPCI_repartition(best_scenario): ...@@ -575,7 +584,7 @@ def EPCI_repartition(best_scenario):
''' '''
list_EPCI_scenario = [] list_EPCI_scenario = []
EPCI_layer = [] EPCI_layer = []
with fiona.open("./EPCI_clean.shp", 'r') as epci : with fiona.open("./EPCI_clean.shp", 'r', encoding='utf-8') as epci :
for f in epci : for f in epci :
EPCI_layer.append(f) EPCI_layer.append(f)
...@@ -687,7 +696,7 @@ def MCMC(scenarInitial, nbr_loop, nbr_scenario, nbr_scenario_to_stock, nbr_parce ...@@ -687,7 +696,7 @@ def MCMC(scenarInitial, nbr_loop, nbr_scenario, nbr_scenario_to_stock, nbr_parce
scenario = [] scenario = []
with fiona.open(scenarInitial, 'r') as input: with fiona.open(scenarInitial, 'r', encoding='utf-8') as input:
for f in input : for f in input :
scenario.append(f) scenario.append(f)
...@@ -827,47 +836,50 @@ def MCMC(scenarInitial, nbr_loop, nbr_scenario, nbr_scenario_to_stock, nbr_parce ...@@ -827,47 +836,50 @@ def MCMC(scenarInitial, nbr_loop, nbr_scenario, nbr_scenario_to_stock, nbr_parce
json_str = json.dumps(scenario_json) json_str = json.dumps(scenario_json)
fid_gz.write(json_str) 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 # 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 cereales = 22
culturesI = 5 culturesI = 5
fourrages = 5 fourrages = 5
fruits_legumes = 3 fruits_legumes = 3
oleagineux = 1 oleagineux = 1
prairies = 59 prairies = 59
proteagineux = 2 proteagineux = 2
# Change to select how many patches you will change from the initial scenario # Change to select how many patches you will change from the initial scenario
nbr_parcelles_a_modifier = 20000 nbr_parcelles_a_modifier = 20000
# Random seed # Random seed
seed = random.randrange(sys.maxsize) seed = random.randrange(sys.maxsize)
rng = random.Random(seed) rng = random.Random(seed)
print(("Seed was:", seed)) print(("Seed was:", seed))
fichier_seed = open("./lastest_seed.txt", "w") fichier_seed = open("./lastest_seed.txt", "w")
fichier_seed.write(str(seed)) 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 = Scenario("./Parcelle_PAT_MCMC.shp",nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, False)
#s1.indice_proximite() #s1.indice_proximite()
#s1.indice_resilience('../Parcelle_PAT/Grille_resilience.shp') #s1.indice_resilience('../Parcelle_PAT/Grille_resilience.shp')
#s1.indice_productivitee("../Parcelle_PAT/Parcelle_PAT_valCad.shp") #s1.indice_productivitee("../Parcelle_PAT/Parcelle_PAT_valCad.shp")
#s1.indice_biodiversite('./compute_dispersal_probabilities_1000.npz') #s1.indice_biodiversite('./compute_dispersal_probabilities_1000.npz')
#s1.indice_biodiversite('./matrix_biodiversity.npz') #s1.indice_biodiversite('./matrix_biodiversity.npz')
#s1.indice_social('./Bassin_de_vie_Pat.shp', './Parcelle_PAT_exploitant.shp') #s1.indice_social('./Bassin_de_vie_Pat.shp', './Parcelle_PAT_exploitant.shp')
#s1.write_in_json() #s1.write_in_json()
#s1.write_in_shape('./testSimulationSurf_MCMC.shp') #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) #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_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_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_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_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) #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 # 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) #scenario_initial = Scenario("./Parcelle_PAT_MCMC.shp",0, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, False, rng)
scenario_initial.write_indices_in_json() scenario_initial = Scenario("../output/PAT_patches/PAT_patches.shp",0, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, False, rng)
MCMC("./Parcelle_PAT_MCMC.shp", 2, 11, 2, nbr_parcelles_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux, 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()
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)
...@@ -68,14 +68,14 @@ class Biodiversity: ...@@ -68,14 +68,14 @@ class Biodiversity:
self.probabilities[fid,i] = 0 self.probabilities[fid,i] = 0
self.probabilities = self.probabilities.tocsr() self.probabilities = self.probabilities.tocsr()
def biodiversity(self, patches): def compute_indicator(self, patches):
''' '''
Return the indice of biodiversity Return the indice of biodiversity
:param patches: List of PAT's patches and their data :param patches: List of PAT's patches and their data
''' '''
surface = 1537520000 # 1 537 520 000 surface = patches.area.sum()
meadowmask = patches["properties"]["rpg_lfgc_s"]=='Prairies' meadowmask = patches["cultgeopat"]=='Prairies'
biodiversity = ((self.probabilities*meadowmask).T*meadowmask).sum() biodiversity = ((self.probabilities*meadowmask).T*meadowmask).sum()
# We multiply by 100000 to get an indices with a value easy to understand for the users # We multiply by 100000 to get an indices with a value easy to understand for the users
return biodiversity/(surface * surface) *100000 return biodiversity/(surface * surface) *100000
...@@ -100,19 +100,14 @@ class Biodiversity: ...@@ -100,19 +100,14 @@ class Biodiversity:
shape=loader['shape'] shape=loader['shape']
) )
def biodiversity( patches, filename): if __name__ == '__main__':
''' import geopandas as gpd
Load the matrix of probabilities and initial parameters from a file (npz format). patches = gpd.GeoDataFrame.from_file("../output/PAT_patches/PAT_patches.shp", encoding='utf-8')
''' matrix_filename='../output/matrix_biodiversity.npz'
patches = gpd.GeoDataFrame.from_features(patches) import os
loader = np.load(filename) biodiv = Biodiversity(1000, 0.001)
dmax = loader['dmax'] if not os.path.isfile(matrix_filename):
epsilon = loader['epsilon'] biodiv.compute_dispersal_probabilities(patches)
probabilities = csr_matrix( biodiv.save_probabilities(matrix_filename)
(loader['data'], loader['indices'], loader['indptr']), biodiv.load_probabilities(matrix_filename)
shape=loader['shape'] print(biodiv.compute_indicator(patches))
)
surface = 1537520000 # 1 537 520 000
meadowmask = patches.rpg_lfgc_s=='Prairies'
biodiversity = ((probabilities*meadowmask).T*meadowmask).sum()
return biodiversity/(surface * surface) * 100000
#!/usr/bin/python3 #!/usr/bin/python3
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
from shapely.geometry import mapping, shape import geopandas as gpd
from shapely.ops import cascaded_union import pandas as pd
import fiona
from fiona import collection
import rtree
import math
# Transform the name of the culture in a code used the algorithm class Productivity:
def return_culture_code(culture):
if culture == "Cereales" :
return 0
if culture == "Cultures industrielles" :
return 1
if culture == "Fourrages" :
return 2
if culture == "Fruits et legumes" :
return 3
if culture == "Oleagineux" :
return 4
if culture == "Prairies" :
return 5
if culture == "Proteagineux" :
return 6
# Give the average value of the cadastral value by squareroot of the scenario for each culture's type
def indice_productivite(layer_scenario, shape_valeur_cadastrale, affichage) :
''' '''
This function calculate the productivity indice. This indicator calculate the productivity indice.
The cadastral's values of the patches of the PAT's territory are used to get an approximation idea of the value of each patch. The cadastral's values of the patches of the PAT's territory are used to get an approximation idea of the value of each patch.
The indice reflect the value of the patch used to cultivate the fruits and vegetables. The indice reflect the value of the patch used to cultivate the fruits and vegetables.
:param layer_scenario: The scenario to analyse as a list
:param shape_valeur_cadastrale: ShapeFile of the PAT's territory with the data about the casdastral's value of each patch as an attribute.
:param affichage: True if we want the display in the console, False is the other case
''' '''
liste_i = {} def __init__(self):
culture_tab = ["Cereales", "Cultures industrielles", "Fourrages", "Fruits et legumes", "Oleagineux", "Prairies", "Proteagineux" ] pass
var_indice_productivite = 0
for i in range(7): # Initialize the list of patches for each type def compute_indicator(self, layer_scenario):
liste_i[i] = [] valeurs_cad = layer_scenario[layer_scenario['cultgeopat']=='Fruits et légumes']
valeurs_cad = valeurs_cad['VALEUR_CAD'] / valeurs_cad.area
for i in range(len(layer_scenario)): # For each patch of the scenario return valeurs_cad.sum() / len(valeurs_cad)
if layer_scenario[i]["properties"]["rpg_lfgc_s"] != None : # If the culture type is not 'None' if __name__ == '__main__':
liste_i[return_culture_code(layer_scenario[i]["properties"]["rpg_lfgc_s"])].append(i) # Add the FID in the list according to his type import geopandas as gpd
patches = gpd.GeoDataFrame.from_file("../output/PAT_patches/PAT_patches.shp", encoding='utf-8')
with fiona.open(shape_valeur_cadastrale, 'r') as layer_valeur_cadastrale: # Open the layer of cadastral values prod = Productivity()
for y in range(7) : # for each type of culture print(prod.compute_indicator(patches))
somme_valeur_cadastrale = 0.0 # Initialize the cadastral value
for i in liste_i[y]: # For each patch of the 'y' type (see return_culture_code function for the signification of y)
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 # Increment the cadastral values(m2)
if len(liste_i[y]) != 0 : # If there is less than one patch possesing the type 'y'
if affichage == True :
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é")
if y == 3 :
var_indice_productivite = somme_valeur_cadastrale/len(liste_i[y])
return var_indice_productivite
#indice_productivite("../Etape1_shp_propre/MCMC/testSimulationSurf_MCMC.shp", "./productivite_shape/Parcelle_PAT_valCad.shp")
#!/usr/bin/python #!/usr/bin/python
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
from shapely.geometry import mapping, shape
from shapely.ops import cascaded_union
import fiona
from fiona import collection
import random
from random import randint
import re
def proximite(shape, fruit_legumes, affichage) : class Proximite:
'''
Function calculating the proximity indice.
The purpose of this indice is to reflect the proximity between productions site and consommation site.
We are calculating it with a comparison of the area of patches cultivating fruits and vegetables and the population of living area (bassin de vie) in the PAT's territory
A living area with a low population will require a lower cultivating area to get a good value.
:param shape: The scenario to analyse as a list
:param fruit_legume: The proportion of area we want in the PAT for the "fruits and vegetables" type
:param affichage: True if we want the display in the console, False is the other case
'''