Commit 65006363 authored by Dumoulin Nicolas's avatar Dumoulin Nicolas
Browse files

Hexagons used in post-processing

parent 60ed93dc
......@@ -131,6 +131,35 @@ def build_PAT_municipalities(pat_municipalities_dir, bdv_filename, bdv_fixes_fil
municipalities = municipalities.merge(admin_express_epci, how='inner', on='CODE_EPCI')
return municipalities
def buildHexagons(patches, muns, hexagonsFilename):
muns.crs={'init':'epsg:2154'}
hexagons = gpd.GeoDataFrame.from_file(hexagonsFilename)
hexagonsPAT=gpd.sjoin(hexagons, muns[['geometry']], op='intersects')
# drop duplicates on index
hexagonsPAT.reset_index(inplace=True)
hexagonsPAT.drop_duplicates(subset='index', inplace=True)
# keep only geometry
hexagonsPAT.drop(['index','left','bottom','right','top','index_right'], axis=1, inplace=True)
hexagonsPAT.reset_index(drop=True, inplace=True)
# Get indexes as a columns for keeping indexes after 'overlay' call
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.to_file('output/hexagons/patches_inter_hexagons.shp')
# computes area of intersection
patches_inter_hexagons['area'] = patches_inter_hexagons.area
# 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']
# 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()
hexagonsPAT.to_file('output/hexagons/hexagons1000PAT.shp')
return hexagonsPAT
if __name__ == '__main__':
import yaml
try:
......@@ -159,3 +188,4 @@ if __name__ == '__main__':
patches.to_file('output/PAT_patches', encoding='utf-8')
gpd.GeoDataFrame(municipalities).to_file('output/municipalities', encoding='utf-8')
buildHexagons(patches, municipalities, resources['hexagons_filename'])
......@@ -70,3 +70,8 @@ MAJIC_parcelles:
file: private/MAJIC_parcelles/MAJIC_parcelles.shp
description: MAJIC landed property data, used for extracting "revenu cadastral"
website: https://www.demarches-simplifiees.fr/dossiers
hexagons:
variable: hexagons_filename
file: hexagons/grid3400.shp
description: Hexagons layer created in QGIS with MMQGIS plugin, with Y spacing = 3400, giving hexagons of almost 1kha
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