MCMC -> création d'un shape + calcul indice de proximité du scenario créé

scenariosAleatoires/MCMC.py

+#!/usr/bin/python
+# -*- 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
+from proximite import proximite 
+
+def return_culture_code(culture):
+	if culture == "Céréales".decode('utf-8') :
+		return 0
+		
+	if culture == "Cultures industrielles" :
+		return 1
+		
+	if culture == "Fourrages" :
+		return 2
+		
+	if culture == "Fruits et légumes".decode('utf-8') :
+		return 3
+		
+	if culture == "Oléagineux".decode('utf-8') :
+		return 4
+		
+	if culture == "Prairies" :
+		return 5
+		
+	if culture == "Protéagineux".decode('utf-8') :
+		return 6
+	
+
+def MCMC_creation(shape, nb_cult_modif, cereales, culturesI, fourrages, fruit_legumes, oleagineux, prairies, proteagineux):
+	with fiona.open(shape, 'r') as input:
+		schema = input.schema.copy()
+		crs = input.crs
+		driver = input.driver
+		
+		resultat_final = []
+		
+		#Correspond au nombre de parcelle de culture que le programme modifiera
+		#Parametre
+		nombre_culture_a_modifie = nb_cult_modif
+		
+		#Surfaces souhaitées pour chacun des types de culture pour la simulation
+		# 0 : Céréales, 1 : culture industrielle, 2 : fourrages, 3 : fruits et légumes, 4 : oléagineux, 5 : prairies, 6 : protéagineux
+		#Valeurs choisi au hasard
+		pourcentage_surface = {}
+		pourcentage_surface[0] = cereales #Céréales
+		pourcentage_surface[1] = culturesI # Cultures industrielles
+		pourcentage_surface[2] = fourrages # Fourrages
+		pourcentage_surface[3] = fruit_legumes # Fruits et légumes
+		pourcentage_surface[4] = oleagineux # Oléagineux
+		pourcentage_surface[5] = prairies # Prairies
+		pourcentage_surface[6] = proteagineux # Protéagineux
+		
+		#Initialisation des surfaces de chaque culture à 0, ceci incrémentera à la création de chaque parcelle lors de la simulation
+		cultures = {}
+	
+		for i in range (7) :
+			cultures[i] = { 'surf': 0.0, 'parcelles': []}
+
+		
+		#Surface totale des parcelles du PAT
+		total_culture = 0.0 #153256
+					
+		i = 0
+		for f in input : 
+			total_culture += f["properties"]["SURF_PARC"]
+			
+			#Initialisation de la surface des cultures de la couche d'entrée
+			culture = f["properties"]["rpg_lfgc_s"]
+			
+			if culture != None :
+				cultures[return_culture_code(culture)]['surf'] += float(f["properties"]["SURF_PARC"])
+				# Le FID de la parcelle est ajouté à la liste 'cultures[typeCulture][parcelles]
+				cultures[return_culture_code(culture)]['parcelles'].append(i)
+			#Incrémentation du fid des parcelles
+			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)
+			#print cultures[i]['surf']*100/total_culture
+					
+		culture_tab = ["Céréales", "Cultures industrielles", "Fourrages", "Fruits et légumes", "Oléagineux", "Prairies", "Protéagineux" ]
+		
+		
+		with fiona.open('MCMC/testSimulationSurf_MCMC.shp','w',driver=driver, crs=crs, schema=schema) as parcelles :
+		  cpt = 0
+		  while cpt < nombre_culture_a_modifie and len(liste_cultures_a_augmenter) > 0 :
+		
+							
+					#Récupération aléatoire d'une parcelle possédant une culture devant diminuer (stocké dans le dictionnaire 'cultures[i]["parcelles"]'
+					rand_cul = random.choice(liste_cultures_a_reduire)
+					
+					random_parcelles= random.choice(cultures[rand_cul]['parcelles'])
+					
+					code_culture_rand = input[random_parcelles]["properties"]["rpg_lfgc_s"]
+				 
+					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"]
+								
+					succes = 0
+					while succes == 0 :		
+						rand_cult = random.choice(liste_cultures_a_augmenter)
+						
+						cultures[rand_cult]['surf'] += input[random_parcelles]["properties"]["SURF_PARC"] 
+						recW["properties"]['rpg_lfgc_s'] = culture_tab[rand_cult]	
+						#Retirer la parcelle 'random_parcelles' de 'liste_nombre'			
+						cultures[return_culture_code(input[random_parcelles]["properties"]["rpg_lfgc_s"])]['parcelles'].remove(random_parcelles)
+						#Retirer la surface de l'ancienne valeur de la parcelle contenu dans 'culture[type_culture]'
+						cultures[return_culture_code(input[random_parcelles]["properties"]["rpg_lfgc_s"])]['surf'] -= float(input[random_parcelles]["properties"]["SURF_PARC"])
+						parcelles.write(recW)
+						resultat_final.append(recW)
+						# Si la proportion de surface d'une culture du territoire dépasse la proportion souhaitée
+						# On l'a retire de la liste des cultures à augmenter pour la mettre dans celle à réduire
+						if 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)
+							
+						# Pareil mais l'inverse pour les cultures qu'il ne faut plus réduire mais qu'il faut augmenter
+						if cultures[return_culture_code(input[random_parcelles]["properties"]["rpg_lfgc_s"])]['surf'] <= pourcentage_surface[rand_cult] :
+							liste_cultures_a_reduire.remove(return_culture_code(input[random_parcelles]["properties"]["rpg_lfgc_s"]))
+							liste_cultures_a_augmenter.append(return_culture_code(input[random_parcelles]["properties"]["rpg_lfgc_s"]))
+							
+						succes = 1
+						cpt += 1
+					# Si toutes les parcelles d'un type de culture ont déjà été modifiées, on retire ce type de culture des listes 
+					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)
+					
+		  for i in range(7) :
+			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"]['INSEE_COM']=input[e]["properties"]["INSEE_COM"]
+				recW["properties"]['rpg_lfgc_s']=input[e]["properties"]["rpg_lfgc_s"]
+				#parcelles.write(recW)
+				resultat_final.append(recW)
+		
+		print
+		print "Pourcentage culture apres : "	
+		for i in range(7) :
+			print cultures[i]['surf']*100/total_culture
+		return resultat_final
+
+		
+#MCMC_creation(shape, nb_cult_modif, cereales, culturesI, fourrages, fruit_legumes, oleagineux, prairies, proteagineux)
+cereales = 20
+culturesI = 5
+fourrages = 10
+fruits_legumes = 15
+oleagineux = 10
+prairies = 15
+proteagineux = 25
+
+nombre_parcelle_a_modifier = 78000
+
+liste_resultat = MCMC_creation('./MCMC/Parcelle_PAT_MCMC.shp', nombre_parcelle_a_modifier, cereales, culturesI, fourrages, fruits_legumes, oleagineux, prairies, proteagineux)
+proximite(liste_resultat, fruits_legumes)
+	

scenariosAleatoires/MCMC/Parcelle_PAT_MCMC.prj

+PROJCS["RGF93_Lambert_93",GEOGCS["GCS_RGF93",DATUM["D_RGF_1993",SPHEROID["GRS_1980",6378137,298.257222101]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["standard_parallel_1",49],PARAMETER["standard_parallel_2",44],PARAMETER["latitude_of_origin",46.5],PARAMETER["central_meridian",3],PARAMETER["false_easting",700000],PARAMETER["false_northing",6600000],UNIT["Meter",1]]
\ No newline at end of file

scenariosAleatoires/MCMC/Parcelle_PAT_MCMC.qpj

+PROJCS["RGF93 / Lambert-93",GEOGCS["RGF93",DATUM["Reseau_Geodesique_Francais_1993",SPHEROID["GRS 1980",6378137,298.257222101,AUTHORITY["EPSG","7019"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6171"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4171"]],PROJECTION["Lambert_Conformal_Conic_2SP"],PARAMETER["standard_parallel_1",49],PARAMETER["standard_parallel_2",44],PARAMETER["latitude_of_origin",46.5],PARAMETER["central_meridian",3],PARAMETER["false_easting",700000],PARAMETER["false_northing",6600000],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],AUTHORITY["EPSG","2154"]]

scenariosAleatoires/proximite.py

+#!/usr/bin/python
+# -*- 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
+
+def proximite(shape, fruit_legumes) :
+		#with fiona.open(shape, 'r') as input:
+		#Surfaces souhaitées pour chacun des types de culture pour la simulation
+		# 0 : Céréales, 1 : culture industrielle, 2 : fourrages, 3 : fruits et légumes, 4 : oléagineux, 5 : prairies, 6 : protéagineux
+		#Valeurs choisi au hasard
+		pourcentage_surface = {}
+
+		pourcentage_surface[3] = fruit_legumes # Fruits et légumes
+
+		
+		#Surface totale des parcelles du PAT
+		total_culture = 153256.0 #
+		
+		#Calculé en fonction de la couche 'ADMIN_EXPRESS-GOG_1-1__SHP__FRA_2018-04-03
+		Population_PAT = 498873
+		
+		#Création du dictionnaire des bassins de vie, code_BV/Population_BV et stockage de la surface actuelle de chaque culture différencié grâce à des emplacement de 0 à 7
+		bv = {}
+		
+		bv["63003"] = {}
+		bv["63010"] = {}
+		bv["63040"] = {}
+		bv["42019"] = {}
+		bv["63050"] = {}
+		bv["43040"] = {} 
+		bv["63113"] = {}
+		bv["63125"] = {}
+		bv["43080"] = {}
+		bv["63178"] = {}
+		bv["63231"] = {}
+		bv["43112"] = {}
+		bv["43157"] = {}
+		bv["63195"] = {} 
+		bv["63210"] = {}
+		bv["42147"] = {}
+		bv["42159"] = {}
+		bv["63284"] = {}
+		bv["63285"] = {}
+		bv["63291"] = {}
+		bv["63300"] = {}
+		bv["42204"] = {}
+		bv["63430"] = {}
+		bv["63214"] = {}
+		bv["63457"] = {}
+		
+
+			
+		bv["63003"][0] = 15130
+		bv["63010"][0] = 6796
+		bv["63040"][0] = 16922
+		bv["42019"][0] = 833 
+		bv["63050"][0] = 3151
+		bv["43040"][0] = 874 
+		bv["63113"][0] = 306616
+		bv["63125"][0] = 12787
+		bv["43080"][0] = 4127
+		bv["63178"][0] = 3834
+		bv["63231"][0] = 7953
+		bv["43112"][0] = 1032
+		bv["43157"][0] = 1398
+		bv["63195"][0] = 630 
+		bv["63210"][0] = 1481
+		bv["42147"][0] = 861
+		bv["42159"][0] = 2397
+		bv["63284"][0] = 16302
+		bv["63285"][0] = 2066
+		bv["63291"][0] = 6319
+		bv["63300"][0] = 50952
+		bv["42204"][0] = 121
+		bv["63430"][0] = 16996
+		bv["63214"][0] = 11121
+		bv["63457"][0] = 8174
+		#Initialisation du stockage de la surface culture actuelle à 0
+		# bv["63003"][15130][0] = 0 ---> La surface de céréale allouée au bassin de vie '63003' en fonction des parcelles simulées est actuellement de 0
+		# Ce compteur incrémentera l'endroit adapté en fonction de la parcelle créée
+		for x in bv :
+			bv[x][1] = {}
+			bv[x][1][3] = 0.0
+
+		f = 0
+		while f < len(shape) :
+			#print shape
+			#print shape[f]["properties"]["rpg_lfgc_s"]
+			culture = shape[f]["properties"]["rpg_lfgc_s"]
+			code_culture = 3
+			if culture == "Fruits et légumes" :
+				print culture
+				code_culture = 3
+				bv[shape[f]["properties"]["Bdv"]][1][code_culture] += shape[f]["properties"]["SURF_PARC"]
+			f += 1
+			
+		
+		total_diff = 0.0
+		nombre_bv = 0
+		
+		for x in bv :
+			popBv = bv[x][0]
+			nombre_bv += 1
+
+			print
+			print x
+			#print
+			
+			Sbv = (total_culture * float(pourcentage_surface[3])/100 * popBv / Population_PAT)
+
+			surf_bv = bv[x][1][3]
+
+			
+		
+			indice_proximite = abs(Sbv - surf_bv ) / Sbv*100
+			if Sbv - surf_bv < 0 :
+				print "Il y a trop de fruits/legumes dans ce bassin de vie"
+				indice_proximite = indice_proximite/100
+			else :
+				print "Il n'y a pas assez de fruits/légumes dans ce bassin de vie"
+				
+			total_diff +=  indice_proximite 
+			print  "Indice de proximité : " + str(indice_proximite )
+				
+				
+		#print total_diff
+		diff_moyenne = total_diff / (nombre_bv) 
+		
+		# Plus on est proche de 0, plus l'indice est bon, si on est proche de 100%, il y a une très mauvaise répartition des fruits et légumes
+		print 
+		print
+		print "Indice de proximité du scénario : " + str(diff_moyenne) + "%"
+		
+#proximite(shape, fruit_legumes)
+#proximite('../Etape1_shp_propre/MCMC/Parcelle_PAT_MCMC.shp', 10 )
+
+