Commit 576b4682 authored by Dumoulin Nicolas's avatar Dumoulin Nicolas
Browse files

Surfaces and additionnal data added for hexagons

parent 65006363
......@@ -6,6 +6,7 @@ import pandas as pd
import geopandas as gpd
import os, glob
import rasterio
import json
def load_MAJIC(majic_dir):
patches = gpd.GeoDataFrame()
......@@ -13,8 +14,6 @@ def load_MAJIC(majic_dir):
patches.append(gpd.GeoDataFrame.from_file(f))
return patches
def build_initial_PAT(rpg_parcelles_filename, rpg_ilots_filename, pat_cultural_classes_filename, majic_dir,
bdalti_dir, bdalti_tiles, pat_municipalities_dir, bdv_filename, bdv_fixes_filename,
adminexpress_com_filename, adminexpress_epci_filename):
......@@ -145,19 +144,33 @@ def buildHexagons(patches, muns, hexagonsFilename):
hexagonsPAT.rename_axis('FID', inplace=True)
hexagonsPAT.reset_index(inplace=True)
# Intersection between hexagons and patches
patches_inter_hexagons = gpd.overlay(patches.reset_index()[['ID_PARCEL','geometry','id_expl','elevation']], hexagonsPAT, how="intersection")
patches_inter_hexagons = gpd.overlay(patches[['ID_PARCEL','geometry','id_expl','elevation', 'cultgeopat','INSEE_COM']], hexagonsPAT.reset_index(), how="intersection")
patches_inter_hexagons.to_file('output/hexagons/patches_inter_hexagons.shp')
# computes area of intersection
patches_inter_hexagons['area'] = patches_inter_hexagons.area
patches_inter_hexagons['area_ha'] = patches_inter_hexagons.area / 10000
# Adding elevation mean
patches_inter_hexagons['tmp'] = patches_inter_hexagons['elevation']*patches_inter_hexagons['area']
grouped = patches_inter_hexagons.groupby('FID')[['tmp','area']].sum()
hexagonsPAT['elevation_mean'] = grouped['tmp']/grouped['area']
patches_inter_hexagons['tmp'] = patches_inter_hexagons['elevation']*patches_inter_hexagons['area_ha']
grouped = patches_inter_hexagons.groupby('FID')[['tmp','area_ha']].sum()
hexagonsPAT['elevation_mean'] = grouped['tmp']/grouped['area_ha']
# remove hexagons without intersection with patches
hexagonsPAT = hexagonsPAT[~hexagonsPAT['elevation_mean'].isnull()]
# Adding number of exploitants
hexagonsPAT['nb_expl'] = patches_inter_hexagons.groupby('FID')['id_expl'].nunique()
# Adding initial cultural surfaces
surfaces = patches_inter_hexagons.groupby('FID').apply(lambda x:x.groupby('cultgeopat')['area_ha'].sum())
surfaces = surfaces.reset_index().pivot(index='FID', columns='cultgeopat', values='area_ha')
surfaces.columns = ['init_surf_cult' + str(col) for col in surfaces.columns]
surfaces.fillna(0, inplace=True)
hexagonsPAT = hexagonsPAT.join(surfaces)
hexagonsPAT.to_file('output/hexagons/hexagons1000PAT.shp')
# Extracting list of municipalities included in each hexagon and the cumulated area of patches
areasByHexByMunDataFrame = patches_inter_hexagons.groupby('FID').apply(lambda h: h.groupby('INSEE_COM').apply(lambda p:p.area.sum() / 10000)).to_frame('area_ha')
# Converting as a dict
areasDict = areasByHexByMunDataFrame.groupby(level=0).apply(lambda dff:dff.xs(dff.name)['area_ha'].to_dict()).to_dict()
areasDict['__README']='This dict gives area of patches (in ha) for each municipality in each hexagon. The first key is the ID of the hexagon.'
# Writing in file
with open("output/hexagons/hexagons_muns.json", "w") as output:
json.dump(areasDict, output)
return hexagonsPAT
if __name__ == '__main__':
......
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