social.py

#!/usr/bin/python3
# -*- coding: utf-8 -*-
import pandas as pd
import numpy as np
import os

class Social:
	'''
	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
	 the conversion cost of the total conversion between the initial scenario and the old one.

	Each transition has a value between 0 and 1.
	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):
		'''
		: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.
		'''
		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

	def compute_patches_transition_cost(self, initial_patches):
		costs=[]
		for bdv,grouped in initial_patches.groupby('Bdv'):
			grouped=grouped.copy()
			surfaces_initial = grouped.groupby('cultgeopat')['SURF_PARC'].sum()
			cult_threshold = surfaces_initial > surfaces_initial.sum() * self.bdv_threshold
		for culture in self.matrice_transition.columns.get_values():
			for row in initial_patches.itertuples():
				cost = 0.0
				# patches that have different farming that others in same Bdv
				if culture not in cult_threshold or not cult_threshold[culture]:
					# patches that have different farming that others of the same farmer
					if culture not in initial_patches[initial_patches['id_expl']==row.id_expl]['cultgeopat']:
						cost = self.matrice_transition.loc[row.cultgeopat,culture]
				costs.append([row.ID_PARCEL, row.cultgeopat, culture, cost])
		self.patches_transition_cost = pd.DataFrame(costs, columns=['ID_PARCEL','initialcult','cultgeopat','cost'])

	def save_patches_transition_cost(self, filename):
		'''
		Save the matrix of costs computed for each patch and initial parameters in a file (npz format).
		'''
		m = self.patches_transition_cost.values
		np.savez(filename, bdv_threshold=self.bdv_threshold, patches_md5sum=self.patches_md5sum, cost_matrix_md5sum=self.cost_matrix_md5sum,
			data=m.data)

	def load_patches_transition_cost(self, filename):
		'''
		Load the matrix of costs and initial parameters from a file (npz format).

		The correspondance will be checked between
		 class attributes and metadata stored in the file
		'''
		loader = np.load(filename)
		error_message = ('Matrix metadata {0} does\'nt correspond to metadata given at indicator construction'
			'\n\t{0} in metadata: {1} while asked: {2}'
			'\n\tmatrix filename: ') + filename
		for k1,k2 in [[self.bdv_threshold,'bdv_threshold'],[self.patches_md5sum,'patches_md5sum'],[self.cost_matrix_md5sum,'cost_matrix_md5sum']]:
			if k2 not in loader or k1!=loader[k2]:
				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'])

	def compute_indicator(self, scenario_patches):
		'''
		:param scenario: The scenario to analyse
		'''
		cost = 0
		patches_altered = scenario_patches[scenario_patches['cultgeopat'] != scenario_patches['init_cult']]
		cost = pd.merge(patches_altered, self.patches_transition_cost,on=['ID_PARCEL','cultgeopat'], how='left')['cost'].sum()
		return 1/(1+cost)

if __name__ == '__main__':
	import geopandas as gpd
	import time
	import os
	import hashlib
	from MCMC import md5sum
	# 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']
	# 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)
	# computing init score
	scenario_patches = patches.copy()
	print(social.compute_indicator(scenario_patches))
	# small reallocation and computing score
	scenario_patches.loc[1,'cultgeopat']='Fruits et légumes'
	# 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 i in range(N):
			scenario_patches.ix[i,'cultgeopat']='Protéagineux'
		print('{} {}'.format(N, social.compute_indicator(scenario_patches)))