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

new abstraction class for indicators

parent 7cad49cc
......@@ -3,3 +3,4 @@ geopandas
rtree
pandas
tqdm
overrides
#!/usr/bin/python
# -*- coding: utf-8 -*-
class Indicator:
def __init__(self, config, initial_patches, patches_md5sum, targetPAT):
pass
def compute_indicator(self, patches):
pass
......@@ -11,7 +11,7 @@ from productivite import Productivity
from indice_biodiversite_2 import Biodiversity
from social import Social
from tqdm import tqdm
import hashlib
from patutils import md5sum
class Scenario:
def __init__(self, patches):
......@@ -44,24 +44,11 @@ class Scenario:
class Indicators:
def __init__(self, config, initial_patches, patches_md5sum, targetPAT):
self._proximity = Proximite(initial_patches, targetPAT['Fruits et légumes'])
self._proximity = Proximite(None, initial_patches, None, targetPAT)
self._resilience = Resilience(config['resilience'], initial_patches)
self._productivity = Productivity()
biodiv_config = config['biodiversity']
matrixfilename = biodiv_config['matrixfilename']
self._biodiversity = Biodiversity(initial_patches, patches_md5sum, biodiv_config['dmax'], biodiv_config['epsilon'])
if not os.path.isfile(matrixfilename):
self._biodiversity.compute_dispersal_probabilities(initial_patches)
self._biodiversity.save_probabilities(matrixfilename)
self._biodiversity.load_probabilities(matrixfilename)
social_config = config['social']
self._social = Social(initial_patches, patches_md5sum, social_config['cost_matrix_filename'],
md5sum(social_config['cost_matrix_filename']), social_config['bdv_threshold'])
patches_costs_filename = social_config['patches_costs_filename']
if not os.path.isfile(patches_costs_filename):
self._social.compute_patches_transition_cost(initial_patches)
self._social.save_patches_transition_cost(patches_costs_filename)
self._social.load_patches_transition_cost(patches_costs_filename)
self._biodiversity = Biodiversity(config['biodiversity'], initial_patches, patches_md5sum)
self._social = Social(config['social'], initial_patches, patches_md5sum)
self.indicators_names = ['proximity', 'resilience', 'productivity', 'biodiversity', 'social']
def compute_indicators(self, patches):
......@@ -75,10 +62,10 @@ class Indicators:
return pd.DataFrame(rows, columns=self.indicators_names)
def proximity(self, patches):
return self._proximity.compute_indicator(patches, False)
return self._proximity.compute_indicator(patches)
def resilience(self, patches):
return self._resilience.compute_indicator(patches, False)
return self._resilience.compute_indicator(patches)
def productivity(self, patches):
return self._productivity.compute_indicator(patches)
......@@ -89,13 +76,6 @@ class Indicators:
def social(self, patches):
return self._social.compute_indicator(patches)
def md5sum(filename):
hash_md5 = hashlib.md5()
with open(filename, "rb") as f:
for chunk in iter(lambda: f.read(4096), b""):
hash_md5.update(chunk)
return hash_md5.hexdigest()
class MCMC:
def __init__(self, mcmc_config_filename):
if not os.path.isfile(mcmc_config_filename):
......
......@@ -4,10 +4,12 @@
import numpy as np
from scipy.sparse import csr_matrix, lil_matrix
from math import exp,log,pow
import os
from tqdm import tqdm
import geopandas as gpd
from Indicator import Indicator
from overrides import overrides
class Biodiversity:
class Biodiversity(Indicator):
'''
This class is used to compute the biodiversity indicator for a set of cultural
patches.
......@@ -22,21 +24,25 @@ 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, initial_patches, patches_md5sum, dmax, epsilon=0.001):
@overrides
def __init__(self, biodiversity_config, initial_patches, patches_md5sum, targetPAT=None):
'''
Creates an indicator with initial parameters
The original file giving the patches is needed for checking the correspondance
between this file and the matrix cached in a file.
:param patches_filename file will be used for computing checksum
:param dmax: the maximum distance threshold for considering a pair of patches
:param epsilon: the minimum probability threshold that should resulted of
the exponential function at distance dmax
'''
self.md5sum = patches_md5sum
self.dmax = dmax
self.epsilon = epsilon
# the maximum distance threshold for considering a pair of patches
self.dmax = biodiversity_config['dmax']
# the minimum probability threshold that should resulted of the exponential function at distance dmax
self.epsilon = biodiversity_config['epsilon']
self.final_result_ratio = pow(initial_patches['SURF_PARC'].sum(), 2) / 100000
matrix_filename = biodiversity_config['matrixfilename']
if not os.path.isfile(matrix_filename):
self.compute_dispersal_probabilities(initial_patches)
self.save_probabilities(matrix_filename)
else:
self.load_probabilities(matrix_filename)
def compute_dispersal_probabilities(self, patches, progress=False):
'''
......@@ -73,6 +79,7 @@ class Biodiversity:
self.probabilities[fid,i] = 0
self.probabilities = self.probabilities.tocsr()
@overrides
def compute_indicator(self, patches):
'''
Return the indice of biodiversity
......@@ -112,20 +119,9 @@ class Biodiversity:
)
if __name__ == '__main__':
import geopandas as gpd
import os
import hashlib
filename = "../output/PAT_patches/PAT_patches.shp"
patches = gpd.GeoDataFrame.from_file(filename, encoding='utf-8')
patches = patches[patches['cultgeopat']!='Non Considérée']
from patutils import load_pat_patches, md5sum
patches_filename = "../output/PAT_patches/PAT_patches.shp"
patches = load_pat_patches(patches_filename)
matrix_filename='../output/matrix_biodiversity_test.npz'
hash_md5 = hashlib.md5()
with open(filename, "rb") as f:
for chunk in iter(lambda: f.read(4096), b""):
hash_md5.update(chunk)
biodiv = Biodiversity(patches, hash_md5.hexdigest(), 1000, 0.001)
if not os.path.isfile(matrix_filename):
biodiv.compute_dispersal_probabilities(patches)
biodiv.save_probabilities(matrix_filename)
biodiv.load_probabilities(matrix_filename)
biodiv = Biodiversity({'dmax': 1000, 'epsilon': 0.001, 'matrixfilename': matrix_filename}, patches, md5sum(patches_filename))
print(biodiv.compute_indicator(patches))
......@@ -2,17 +2,21 @@
# -*- coding: utf-8 -*-
import geopandas as gpd
import pandas as pd
from Indicator import Indicator
from overrides import overrides
class Productivity:
class Productivity(Indicator):
'''
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 indice reflect the value of the patch used to cultivate the fruits and vegetables.
'''
def __init__(self):
@overrides
def __init__(self, social_config=None, initial_patches=None, patches_md5sum=None, targetPAT=None):
pass
@overrides
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['SURF_PARC']
......@@ -20,6 +24,7 @@ class Productivity:
if __name__ == '__main__':
import geopandas as gpd
patches = gpd.GeoDataFrame.from_file("../output/PAT_patches/PAT_patches.shp", encoding='utf-8')
from patutils import load_pat_patches
patches = load_pat_patches("../output/PAT_patches/PAT_patches.shp")
prod = Productivity()
print(prod.compute_indicator(patches))
#!/usr/bin/python
# -*- coding: utf-8 -*-
from Indicator import Indicator
from overrides import overrides
class Proximite:
class Proximite(Indicator):
'''
Function calculating the proximity indice.
The purpose of this indice is to reflect the proximity between productions site and consommation site.
......@@ -10,14 +12,17 @@ class Proximite:
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) :
@overrides
def __init__(self, config, initial_patches, patches_md5sum=None, targetPAT=None):
# def __init__(self, config, initial_patches, patches_md5sum, surf_totale_cible) :
Population_PAT = initial_patches.groupby('INSEE_COM')['POPULATION'].first().sum() #498873
popByBdv = initial_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
self.targetSurfByBdv = targetPAT['Fruits et légumes'] * popByBdv/Population_PAT
def compute_indicator(self, patches, affichage):
@overrides
def compute_indicator(self, patches):
'''
: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
......@@ -36,11 +41,12 @@ class Proximite:
return result
if __name__ == '__main__':
import geopandas as gpd
import pandas as pd
patches = gpd.GeoDataFrame.from_file("../output/PAT_patches/PAT_patches.shp", encoding='utf-8')
from patutils import load_pat_patches
patches_filename = "../output/PAT_patches/PAT_patches.shp"
patches = load_pat_patches(patches_filename)
target = pd.read_csv('../resources/targetPAT.csv', sep=';',index_col=0)
targetRatio = (target['2050']-target['2016'])/target['2016']
targetPAT = patches.groupby('cultgeopat')['SURF_PARC'].sum()*(1+targetRatio)
prox = Proximite(patches, targetPAT['Fruits et légumes'])
print(prox.compute_indicator(patches, True))
prox = Proximite(None, patches, None, targetPAT)
print(prox.compute_indicator(patches))
......@@ -2,15 +2,18 @@
# -*- coding: utf-8 -*-
import geopandas as gpd
import math
from Indicator import Indicator
from overrides import overrides
class Resilience:
def __init__(self, grille, initial_patches):
class Resilience(Indicator):
@overrides
def __init__(self, resilience_config, initial_patches, patches_md5sum=None, targetPAT=None):
'''
Initialize the indicator.
:param grille: Grid divising the PAT in squares
'''
self.grid = gpd.GeoDataFrame.from_file(grille)
self.grid = gpd.GeoDataFrame.from_file(resilience_config)
intersection = gpd.sjoin(
self.grid[['id','geometry']],
initial_patches[['geometry','SURF_PARC']],
......@@ -20,7 +23,8 @@ class Resilience:
self.surfaces = intersection[['id','SURF_PARC']].groupby('id').sum()
# The parameters are : the grid layer, a list corresponding to a scenario and a boolean value True or False signifying if we want to see the printing lines or not
def compute_indicator(self, patches, affichage=False):
@overrides
def compute_indicator(self, patches):
'''
This function calculate the resiliency indice
......@@ -37,18 +41,19 @@ class Resilience:
'''
resilience = 0
for cell,group_by_cell in self.grouped_by_cell:
surf_tot = self.surfaces.loc[cell].values[0]
intersection = patches[['ID_PARCEL','cultgeopat','SURF_PARC']].loc[group_by_cell['index_right']]
prop_cult = intersection.groupby('cultgeopat')['SURF_PARC'].sum() / surf_tot
res_cell = - (prop_cult * prop_cult.apply(math.log2)).sum()
surf_tot = self.surfaces.loc[cell].values[0] # 0.085ms
intersection = patches[['ID_PARCEL','cultgeopat','SURF_PARC']].loc[group_by_cell['index_right']] # 4.26ms
prop_cult = intersection.groupby('cultgeopat')['SURF_PARC'].sum() / surf_tot # 0.5ms
res_cell = - (prop_cult * prop_cult.apply(math.log2)).sum() # 0.22ms
resilience += res_cell
return resilience / self.size
if __name__ == '__main__':
import geopandas as gpd
patches = gpd.GeoDataFrame.from_file("../output/PAT_patches/PAT_patches.shp", encoding='utf-8')
from patutils import load_pat_patches
patches_filename = "../output/PAT_patches/PAT_patches.shp"
patches = load_pat_patches(patches_filename)
import time
res = Resilience("Parcelle_PAT/Grille_resilience.shp")
res = Resilience("Parcelle_PAT/Grille_resilience.shp", patches)
start = time.time()
print(res.compute_indicator(patches))
elapsed = time.time()-start
......
......@@ -3,8 +3,11 @@
import pandas as pd
import numpy as np
import os
from patutils import md5sum
from Indicator import Indicator
from overrides import overrides
class Social:
class Social(Indicator):
'''
This function is used to calculate the social indice.
The idea is to propose a transition matrix for each type of culture and use it to visualize
......@@ -14,18 +17,26 @@ class Social:
These values were defined during reunions and with the review of our experts.
'''
def __init__(self, initial_patches, patches_md5sum, cost_matrix_filename, cost_matrix_md5sum, bdv_threshold):
@overrides
def __init__(self, social_config, initial_patches, patches_md5sum, targetPAT=None):
'''
:param initial_patches Initial state of the patches (attributes cultugeopat and SURF_PARC required)
:param patches_md5sum md5 checksum of the input file of patches
:param cost_matrix_filename input file containing the matrix in CSV format (separator: \t)
:param bdv_threshold proportion of a cultural type needed for considering that this cultural type
is suffisantly present in the locality for ignoring the conversion cost.
:param social_config dict giving parameters
:param targetPAT is ignored in this implementation
'''
self.patches_md5sum = patches_md5sum
self.matrice_transition = pd.read_csv(cost_matrix_filename, sep='\t',index_col=0)
self.cost_matrix_md5sum = cost_matrix_md5sum
self.bdv_threshold = bdv_threshold
self.matrice_transition = pd.read_csv(social_config['cost_matrix_filename'], sep='\t',index_col=0)
self.cost_matrix_md5sum = md5sum(social_config['cost_matrix_filename'])
# proportion of a cultural type needed for considering that this cultural type
# is suffisantly present in the locality for ignoring the conversion cost.
self.bdv_threshold = social_config['bdv_threshold']
costs_filename = social_config['patches_costs_filename']
if not os.path.isfile(costs_filename):
self.compute_patches_transition_cost(initial_patches)
self.save_patches_transition_cost(costs_filename)
else:
self.load_patches_transition_cost(costs_filename)
def compute_patches_transition_cost(self, initial_patches):
costs=[]
......@@ -68,6 +79,7 @@ class Social:
raise ValueError(error_message.format(k2, dict(loader.items()).get(k2, 'NA'), k1))
self.patches_transition_cost = pd.DataFrame(loader['data'], columns=['ID_PARCEL','initialcult','cultgeopat','cost'])
@overrides
def compute_indicator(self, scenario_patches):
'''
:param scenario: The scenario to analyse
......@@ -78,24 +90,17 @@ class Social:
return 1/(1+cost)
if __name__ == '__main__':
import geopandas as gpd
import time
import os
import hashlib
from MCMC import md5sum
from patutils import load_pat_patches
# loading and initializing patches data
patches_filename = "../output/PAT_patches/PAT_patches.shp"
matrix_filename = '../resources/matrice_transition.csv'
patches = gpd.GeoDataFrame.from_file(patches_filename, encoding='utf-8')
patches = patches[patches['cultgeopat']!='Non Considérée']
patches['init_cult'] = patches['cultgeopat']
patches = load_pat_patches(patches_filename)
# init indicator
social = Social(patches, md5sum(patches_filename), matrix_filename, md5sum(matrix_filename), 0.1)
costs_filename = '../output/patches_transition_costs.npz'
if not os.path.isfile(costs_filename):
social.compute_patches_transition_cost(patches)
social.save_patches_transition_cost(costs_filename)
social.load_patches_transition_cost(costs_filename)
social = Social({
'cost_matrix_filename': matrix_filename, 'bdv_threshold': 0.1,
'patches_costs_filename': '../output/patches_transition_costs.npz'
}, patches, md5sum(patches_filename))
# computing init score
scenario_patches = patches.copy()
print(social.compute_indicator(scenario_patches))
......@@ -104,7 +109,7 @@ if __name__ == '__main__':
# should give 1/(1+1) the cost for one transition to Fruits et légumes
print(social.compute_indicator(scenario_patches))
# computing with more reallocation
for N in [2000,5000,10000]:
for N in [200,500,2000]:#,10000]:
for i in range(N):
scenario_patches.ix[i,'cultgeopat']='Protéagineux'
print('{} {}'.format(N, social.compute_indicator(scenario_patches)))
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment