Fix in proximity: score computed related to the target demand

optimizations on area computation

scenariosAleatoires/MCMC_class.py

@@ -94,7 +94,7 @@ class Scenario :
 			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)
+			self.proximite = Proximite(self.initial_patches, self.mcmc_config['indicators_config']['proximite']['surf_totale_cible'])
 
 			# If ScenarioInitial is a filename -> transform it in a list
 			if liste == False :
@@ -321,7 +321,8 @@ class Scenario :
 	These functions returns values of indices of this scenario
 	'''
 	def indice_proximite(self) :
-		indice_prox = self.proximite.compute_indicator(self.resultat_final, self.fruits_legumes, False)
+		patches = gpd.GeoDataFrame.from_features(self.resultat_final)
+		indice_prox = self.proximite.compute_indicator(patches, self.fruits_legumes, False)
 		#print "Résultats de l'indice de proximité : " + str(indice_prox)
 		return indice_prox
 
@@ -343,7 +344,7 @@ class Scenario :
 		#Matrix containing all the data about biodiversity. Created with Numpy and the formula of Equivalent Connected Area (Saura et al., 2011)
 		biodiv_config = self.mcmc_config['indicators_config']['biodiversity']
 		matrixfilename = biodiv_config['matrixfilename']
-		biodiversity = Biodiversity(biodiv_config['dmax'], biodiv_config['epsilon'])
+		biodiversity = Biodiversity(self.initial_patches, biodiv_config['dmax'], biodiv_config['epsilon'])
 		biodiversity.load_probabilities(matrixfilename)
 		patches = gpd.GeoDataFrame.from_features(self.resultat_final)
 		indice_bio = biodiversity.compute_indicator(patches)

scenariosAleatoires/indice_biodiversite_2.py

@@ -3,7 +3,7 @@
 
 import numpy as np
 from scipy.sparse import csr_matrix, lil_matrix
-from math import exp,log
+from math import exp,log,pow
 from tqdm import tqdm
 import geopandas as gpd
 
@@ -22,7 +22,7 @@ class Biodiversity:
 	function P(i,j) = exp(-l*d). The coefficient l (lambda) is determined for a given
 	minimal probability at the maximal distance d.
 	'''
-	def __init__(self, dmax, epsilon=0.001):
+	def __init__(self, initial_patches, dmax, epsilon=0.001):
 		'''
 		Creates an indicator with initial parameters
 
@@ -32,6 +32,7 @@ class Biodiversity:
 		'''
 		self.dmax = dmax
 		self.epsilon = epsilon
+		self.final_result_ratio = pow(initial_patches['SURF_PARC'].sum(), 2) / 100000
 
 	def compute_dispersal_probabilities(self, patches, progress=False):
 		'''
@@ -74,11 +75,10 @@ class Biodiversity:
 
 		:param patches: List of PAT's patches and their data
 		'''
-		surface = patches.area.sum()
 		meadowmask = patches["cultgeopat"]=='Prairies'
 		biodiversity = ((self.probabilities*meadowmask).T*meadowmask).sum()
 		# We multiply by 100000 to get an indices with a value easy to understand for the users
-		return biodiversity/(surface * surface) *100000
+		return biodiversity/self.final_result_ratio
 
 	def save_probabilities(self, filename):
 		'''
@@ -105,7 +105,7 @@ if __name__ == '__main__':
 	patches = gpd.GeoDataFrame.from_file("../output/PAT_patches/PAT_patches.shp", encoding='utf-8')
 	matrix_filename='../output/matrix_biodiversity.npz'
 	import os
-	biodiv = Biodiversity(1000, 0.001)
+	biodiv = Biodiversity(patches, 1000, 0.001)
 	if not os.path.isfile(matrix_filename):
 		biodiv.compute_dispersal_probabilities(patches)
 		biodiv.save_probabilities(matrix_filename)

scenariosAleatoires/productivite.py

@@ -15,7 +15,7 @@ class Productivity:
 
 	def compute_indicator(self, layer_scenario):
 		valeurs_cad = layer_scenario[layer_scenario['cultgeopat']=='Fruits et légumes']
-		valeurs_cad = valeurs_cad['VALEUR_CAD'] / valeurs_cad.area
+		valeurs_cad = valeurs_cad['VALEUR_CAD'] / valeurs_cad['SURF_PARC']
 		return valeurs_cad.sum() / len(valeurs_cad)
 
 if __name__ == '__main__':

scenariosAleatoires/proximite.py

@@ -2,73 +2,35 @@
 # -*- coding: utf-8 -*-
 
 class Proximite:
-
-	def __init__(self, initial_patches) :
-		self.total_culture = initial_patches.area.sum() #153256.0
-		self.Population_PAT = initial_patches.groupby('INSEE_COM')['POPULATION'].first().sum() #498873
-		# Create a dictionary of living area (bassin de vie) bv[code_Living_Area] = Population_Living_Area
-		self.bv = {
-			bdv: {
-				0: initial_patches[initial_patches['Bdv']==bdv].groupby('INSEE_COM')['POPULATION'].first().sum(),
-				1: {3: 0.0}
-			}
-			for bdv in initial_patches.groupby('Bdv')['Bdv'].first().values
-		}
-
-	def compute_indicator(self, shape, fruit_legumes, affichage):
-		# 3 is the code for Fruits and vegetables
-		self.pourcentage_surface = {}
-		self.pourcentage_surface[3] = fruit_legumes
-		# Incrementing the fruits and vegetables area for each living area (bassin de vie)
-		f = 0
-		cpt = 0
-		while f < len(shape) :
-			culture = shape[f]["properties"]["cultgeopat"]
-			code_culture = 3
-			if culture == "Fruits et legumes" :
-				self.bv[shape[f]["properties"]["Bdv"]][1][code_culture] += shape[f]["properties"]["SURF_PARC"]
-			f += 1
-		total_diff = 0.0
-		nombre_bv = 0
-		# Calculating the indice for each living area
-		for x in self.bv :
-			popBv = self.bv[x][0]
-			nombre_bv += 1
-			if affichage == True:
-				print()
-				print(x)
-			Sbv = (self.total_culture * float(self.pourcentage_surface[3])/100 * popBv / self.Population_PAT)
-
-			surf_bv = self.bv[x][1][3]
-
-			if Sbv - surf_bv < 0 :
-				indice_proximite = 0
-				if affichage == True:
-					print("Il y a trop de fruits/legumes dans ce bassin de vie")
-				indice_proximite = indice_proximite/100
-			else :
-				indice_proximite = (Sbv - surf_bv ) / Sbv*100
-				if affichage == True:
-					print("Il n'y a pas assez de fruits/légumes dans ce bassin de vie")
-
-			total_diff +=  indice_proximite * popBv
-			if affichage == True:
-				print("Indice de proximité : " + str(indice_proximite ))
-
-		# Average of indices of each living area
-		diff_moyenne = total_diff / (self.Population_PAT)
-
-		# Near 0 -> good value, the nearest from 100, the baddest the value and baddest the repartition of fruits and vegetables is.
-		if affichage == True:
-			print()
-			print("Indice de proximité du scénario : " + str(diff_moyenne) + "%")
-
-		return diff_moyenne
+	'''
+	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.
+	'''
+	def __init__(self, initial_patches, surf_totale_cible) :
+		Population_PAT = initial_patches.groupby('INSEE_COM')['POPULATION'].first().sum() #498873
+		popByBdv = patches.groupby('Bdv').apply(lambda x:x.groupby('INSEE_COM')['POPULATION'].first().sum())
+		# OK but slower
+		# popByBdv = patches.groupby(['Bdv','INSEE_COM']).first().groupby('Bdv')['POPULATION'].sum()
+		self.targetSurfByBdv = surf_totale_cible * popByBdv/Population_PAT
+
+	def compute_indicator(self, patches, fruit_legumes, affichage):
+		'''
+		: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
+		'''
+		# get area of "fruits et légumes" in the scenario
+		flSurfByBdv = patches[patches['cultgeopat']=='Fruits et légumes'].groupby('Bdv').agg({'SURF_PARC':sum})
+		result = flSurfByBdv.div(self.targetSurfByBdv,axis=0)
+		# put 0 where target is reached
+		result = result.where(result<1,0).sum() / result.where(result<1).count()
+		return result.values[0]
 
 if __name__ == '__main__':
 	import geopandas as gpd
 	patches = gpd.GeoDataFrame.from_file("../output/PAT_patches/PAT_patches.shp", encoding='utf-8')
-	prox = Proximite(patches)
-	del patches['geometry']
-	shape2 = [{'properties':patches.ix[i].to_dict()} for i in range(patches.shape[0])]
-	print(prox.compute_indicator(shape2, 3, True))
+	prox = Proximite(patches, 20165939.605135553)
+	print(prox.compute_indicator(patches, 3, True))