Commit 8a8b4d7b authored by Seneque Colin's avatar Seneque Colin
Browse files

Add of indices,MCMC and some samples files with english commentaries

parent 92cbc481
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
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"]]
UTF-8
\ No newline at end of file
PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]
\ No newline at end of file
PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]
#!/usr/bin/python
# -*- coding: utf-8 -*-
import numpy as np
from scipy.sparse import csr_matrix, lil_matrix
from math import exp,log
from tqdm import tqdm
import geopandas as gpd
def EPCI_culture_repartition(EPCI, scenario):
'''
This function calculate the area for each type of culture for each EPCI and return it as a dictionary.
It will be mostly used in the web application.
:param EPCI: ShapeFile of EPCI converts to a list
:param scenario:
'''
EPCI_layer = gpd.GeoDataFrame.from_features(EPCI)
scenario_layer = gpd.GeoDataFrame.from_features(scenario)
index = scenario_layer.sindex
dict_EPCI = {}
for i in range(len(EPCI_layer)) :
dict_EPCI[i] = {}
dict_EPCI[i]["properties"] = {}
#dict_EPCI[i]["geometry"] = EPCI_layer.iloc.geometry
dict_EPCI[i]["properties"]['CODE_EPCI'] = EPCI_layer.iloc[i].CODE_EPCI
dict_EPCI[i]["properties"]['NOM_EPCI'] = EPCI_layer.iloc[i].NOM_EPCI
dict_EPCI[i]["properties"]['POPULATION'] = EPCI_layer.iloc[i].POPULATION
dict_EPCI[i]["properties"]['Cereales'] = 0
dict_EPCI[i]["properties"]['Cultures industrielles'] = 0
dict_EPCI[i]["properties"]['Fourrages'] = 0
dict_EPCI[i]["properties"]['Fruits et legumes'] = 0
dict_EPCI[i]["properties"]['Oleagineux'] = 0
dict_EPCI[i]["properties"]['Prairies'] = 0
dict_EPCI[i]["properties"]['Proteagineux'] = 0
dict_EPCI[i]["properties"]['None'] = 0
for i in range(len(EPCI_layer)):
geom_i = EPCI_layer.iloc[i].geometry
for fid in list(index.intersection(geom_i.bounds)) :
if scenario_layer.iloc[fid].rpg_lfgc_s == None :
dict_EPCI[i]["properties"]["None"] += scenario_layer.iloc[fid].SURF_PARC
else :
dict_EPCI[i]["properties"][scenaro_layer.iloc[fid].rpg_lfgc_s] += scenario_layer.iloc[fid].SURF_PARC
return dict_EPCI
#!/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
# Transforme le nom du type de culture en un code utilisé dans l'algorithme
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
# Fonction de création d'un scénario aléatoire modifiant un nombre donné en paramètre 'nb_cult_modif' de parcelles
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 parcelles de culture que le programme modifiera
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 (8) :
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 += float(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 :
#total_culture += float(f["properties"]["SURF_PARC"])
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)
else :
cultures[7]['parcelles'].append(i)
cultures[7]['surf'] += float(f["properties"]["SURF_PARC"])
#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)
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) + len(liste_cultures_a_reduire) > 1 :
#Récupération aléatoire d'une parcelle possédant une culture devant diminuer (stocké dans le dictionnaire 'cultures[i]["parcelles"]'
if len(liste_cultures_a_reduire) != 0 :
rand_cul = random.choice(liste_cultures_a_reduire)
else :
# Cas où les proportions des parcelles sont déjà optimisé mais où il reste des parcelles dans plusieurs types
# On ajoute un type aléatoire des types à augmenter (qui sera extrèmement proche de la valeur souhaité) dans celle des valeurs à réduire afin de continuer l'algorithme jusqu'à réduire au maximum le nombre de parcelle restante
# Permet d'optimiser au maximum la randomisation du résultat
rand_cul = random.choice(liste_cultures_a_augmenter)
liste_cultures_a_augmenter.remove(rand_cul)
liste_cultures_a_reduire.append(rand_cul)
if len(liste_cultures_a_augmenter) == 0 :
rand_augmenter = random.choice(liste_cultures_a_reduire)
liste_cultures_a_augmenter.append(rand_augmenter)
liste_cultures_a_reduire.remove(rand_augmenter)
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'] += float(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 float(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 float(cultures[return_culture_code(input[random_parcelles]["properties"]["rpg_lfgc_s"])]['surf']*100/total_culture) <= pourcentage_surface[return_culture_code(input[random_parcelles]["properties"]["rpg_lfgc_s"])] :
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
# Décommenter pour voir en 'live' les changements des proportions des types dans la console
#for i in range(7) :
#print cultures[i]['surf']*100/total_culture
# 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)
# Toutes les parcelles restantes sont ajoutées dans la nouvelle couche, les parcelles avec un type vide non traitées précédement sont ajoutées dans cette phase
for i in range(8) :
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 MCMC : ")
for i in range(7) :
print (cultures[i]['surf']*100/total_culture)
return resultat_final
cereales = 10
culturesI = 5
fourrages = 10
fruits_legumes = 10
oleagineux = 10
prairies = 45
proteagineux = 10
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)
This diff is collapsed.
UTF-8
\ No newline at end of file
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
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"]]
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