Commit 82bb61e8 authored by Nicolas Dumoulin's avatar Nicolas Dumoulin
Browse files

Post-process script

parent 697b4667
import pandas as pd
import geopandas as gpd
import patutils
import itertools
import os
from tqdm import tqdm
from sklearn import preprocessing
import multiprocessing as mp
from functools import partial
import glob
def cultModificationForScenario(coeffs, cult, cultToDecrease, patches, scenarios, scores):
scen = scenarios.copy()
# Keeping only modified patches for each scenario
scen[scen == patches.loc[scen.columns,'init_cult']] = -1
if cult in cultToDecrease:
# For cultural types that decrease, we only keep modified patches which initial cultural type equals to current cultural type
scen.loc[:,patches.loc[scen.columns,'init_cult'] != cult] = -1
# For cultural types that increase, we only study patches with given cultural
scen[scen != cult] = -1
# Putting 0 where patches is concerned by cultural modification
scen.mask(scen > -1, 0, inplace=True)
# Replacing -1 by 0 and 0 by 1
scen = pd.eval('scen + 1')
return scen
def buildHexagons(combi, cult, cultToDecrease, patches, scenarios, scores, patches_inter_hexagons, onlyFirstDecile = False):
coeffs = pd.Series(combi, index=patutils.indicators)
scores = scores[coeffs.index].copy()
scen = cultModificationForScenario(coeffs, cult, cultToDecrease, patches, scenarios, scores)
# Now we have "1" for patches modified and "0" for unmodified
# counting the number of reallocation per patch
patches_reallocation_counts = scen.sum(axis=0)
# Reverting so that best values are the highest one (results are optimized for minimum)
scores = 1 / scores
# reshaping scores and applying coeffs
#scores = scores / scores.max()
scores[coeffs.index] = preprocessing.MinMaxScaler().fit_transform(scores)
scores = scores.mul(coeffs, axis=1).sum(axis=1)
if onlyFirstDecile:
deciles = pd.cut(scores.rank(method='first'), 10, labels=False)
scores = scores[deciles == deciles.max()]
scen = scen.loc[scores.index,:]
# Now we have the sum of scores for each patch
# Reshaping again the aggregated scores, so the best scenario score to 1
scores_by_scen = (scores-scores.min())/(scores.max()-scores.min())
# Computing the mean of scores for scenarios where a reallocation has occured (>0)
scores_by_patch = (scen.mul(scores_by_scen, axis=0).sum(axis=0) / patches_reallocation_counts).fillna(0)
# merging with geometry of intersection patches/hexagons
# Assigning patches scores to the intersections patches/hexagons
# hex_scores = gpd.GeoDataFrame(pd.merge(scores.to_frame(key), patches_inter_hexagons, on='ID_PARCEL'))
hex_scores = patches_inter_hexagons.merge(scores_by_patch.to_frame(key), on='ID_PARCEL')
# Multiplying by the area of the intersection patch/hexagon
hex_scores[key] = hex_scores[key] * hex_scores.area
hex_scores = hex_scores.groupby('FID')[key].sum()
# Rescaling the surfaces of hexagons with the mean of the surfaces of scenarios
modifiedAreaFromScen = scen.mul(patches.geometry.area, axis=1).sum(axis=1).mean()
modifiedAreaFromHexagons = hex_scores.sum()
hex_scores = hex_scores * modifiedAreaFromScen / modifiedAreaFromHexagons
# return hex_scores
# Now adding nb of exploitants concerned by the transition
nbexpl_byscen = scen.mul(patches['id_expl']).apply(lambda s: patches_inter_hexagons[['ID_PARCEL','FID']].merge(s.pipe(lambda x:x[x>0]).rename('count',copy=False).reset_index(), on='ID_PARCEL').groupby('FID')['count'].nunique(),axis=1).fillna(0).astype(int)
hex_nbexpl = nbexpl_byscen.mul(scores_by_scen,axis=0).sum(axis=0)/(nbexpl_byscen/nbexpl_byscen).mul(scores_by_scen, axis=0).sum()
result = pd.concat([hex_scores.to_frame(), hex_nbexpl.to_frame(key+'_nb_expl')], axis=1, sort=False)
result.fillna(0, inplace=True)
# Converting from m² to ha
return result / 10000
def buildHexagons_mp(combi, cult, cultToDecrease, onlyFirstDecile=False):
global patches
global scenarios
global scores
global patches_inter_hexagons
return buildHexagons(combi, cult, cultToDecrease, patches, scenarios, scores, patches_inter_hexagons, onlyFirstDecile)
def post_process(outputdir, nbprocess, nblevels=2):
global patches
global scenarios
global scores
global patches_inter_hexagons
patches = patutils.load_pat_patches('../output/PAT_patches/PAT_patches.shp')
scenarios = pd.read_hdf(sorted(glob.glob(outputdir + '/cult_iter_*.h5'))[-1])
scores = pd.read_csv(sorted(glob.glob(outputdir + '/scores_iter_*.csv'))[-1], index_col=0)
hexagons = gpd.GeoDataFrame.from_file('../output/hexagons/hexagons1000PAT.shp')
# Caching the intersection between hexagons and municipalities
hexagons.set_index('FID', inplace=True)
patches_inter_hexagons = gpd.GeoDataFrame.from_file('../output/hexagons/patches_inter_hexagons.shp')
# For cultural types to decrease, we will count the patches that have leave the considered cultural type
targetPAT = pd.read_csv('../resources/targetPAT.csv',delimiter=';',index_col=0).rename(index=patutils.code_cultgeopat)
targetRatio = (targetPAT['2050']-targetPAT['2016'])/targetPAT['2016']
cultToDecrease = targetRatio[targetRatio<0].index
combis = list(itertools.product(range(nblevels),repeat=len(patutils.indicators)))
if nbprocess==0:
for cult in tqdm(patutils.reverse_code_cultgeopat):
for combi in combis:
result = buildHexagons(combi, cult, cultToDecrease, patches, scenarios, scores, patches_inter_hexagons, onlyFirstDecile=True)
for cult in tqdm(patutils.reverse_code_cultgeopat):
buildHexForCult = partial(buildHexagons_mp, cult=cult, cultToDecrease=cultToDecrease, onlyFirstDecile=True)
with mp.Pool(nbprocess) as p:
# result =, combis)
result = list(tqdm(p.imap(buildHexForCult, combis), total=len(combis)))
gdf = gpd.GeoDataFrame(pd.merge(pd.concat(hexscores, axis=1), hexagons, how='left', on='FID'))
# normalization, skipping the columns 'FID' (index 0) and 'geometry' (index -1)
for c in gdf.columns[1:-1]:
if int(c[0]) in cultToDecrease:
gdf[c] = - gdf[c]
gdf = joinInitialHexData(gdf)
# FIXME geometry seems broken after joining
gdf.reset_index(inplace=True) # for having FID in geojson output
# TODO reproject from 2154 to 3857: crs parameter of fonction to_file (from fiona)
gdf.to_file(outputdir+'/hexscores.geojson', driver='GeoJSON') #, crs={'init': 'epsg:3857', 'no_defs': True})
def joinInitialHexData(hexScores):
hexInit = gpd.GeoDataFrame.from_file('../output/hexagons/hexagons1000PAT.shp')
hexInit.set_index('FID', inplace=True)
return hexScores.join(hexInit)
if __name__ == '__main__':
mp.set_start_method('fork', force=True)
import sys
if len(sys.argv) < 2:
print('Post-Process optimization results and generate hexagons map')
print('Arguments required: [output directory] [optional: nb process]')
if len(sys.argv) == 2:
post_process(sys.argv[1], nbprocess=0)
post_process(sys.argv[1], nbprocess=int(sys.argv[2]))
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