Elevation data added for patches (with BDAlti and module rasterio)

fix in dataframe merge for preserving CRS

PATBuilder.py

@@ -5,6 +5,7 @@ import numpy as np
 import pandas as pd
 import geopandas as gpd
 import os, glob
+import rasterio
 
 def load_MAJIC(majic_dir):
     patches = gpd.GeoDataFrame()
@@ -12,8 +13,10 @@ 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,
-    pat_municipalities_dir, bdv_filename, bdv_fixes_filename,
+    bdalti_dir, bdalti_tiles, pat_municipalities_dir, bdv_filename, bdv_fixes_filename,
     adminexpress_com_filename, adminexpress_epci_filename):
     '''
     This method build the initial patches set by consolidating the data sources
@@ -29,7 +32,7 @@ def build_initial_PAT(rpg_parcelles_filename, rpg_ilots_filename, pat_cultural_c
     '''
     municipalities = build_PAT_municipalities(pat_municipalities_dir, bdv_filename, bdv_fixes_filename,
         adminexpress_com_filename, adminexpress_epci_filename)
-    municipalities = gpd.GeoDataFrame(municipalities[['INSEE_COM','Nom','CODE_EPCI','geometry','POPULATION','BV2012','LIBBV2012']])
+    municipalities = municipalities[['INSEE_COM','Nom','CODE_EPCI','NOM_EPCI', 'geometry','POPULATION','BV2012','LIBBV2012']]
     rpg = gpd.GeoDataFrame.from_file(rpg_parcelles_filename, encoding='utf-8')[['ID_PARCEL','CODE_CULTU','geometry']]
     # we will use the centroid for doing the sjoin with municipalities
     rpg['patches'] = rpg.geometry
@@ -78,8 +81,19 @@ def build_initial_PAT(rpg_parcelles_filename, rpg_ilots_filename, pat_cultural_c
         'Classe GéoPAT':'cultgeopat',
         'BV2012':'Bdv'
         })
-    # delete unconsidered rows
-    # patches_expl = patches_expl[patches_expl['cultgeopat']!='Non Considérée']
+    # Adding elevation from BD Alti
+    rasters = [] # reading tiles
+    for tile in bdalti_tiles:
+        dataset = rasterio.open(bdalti_dir+tile+'.asc')
+        rasters.append([dataset, dataset.read(1)])
+    def getElevation(point): # Function for getting the elevation from a point
+        # Looping on every loaded rasters
+        for raster in rasters:
+            coords = raster[0].index(point.x, point.y)
+            # If coords are in bounding box, returning the elevation
+            if coords[0] >= 0 and coords[1] >= 0 and coords[0] < 1000 and coords[1] < 1000:
+                return raster[1][coords]
+    patches_expl['elevation'] = patches_expl.geometry.centroid.apply(getElevation)
     # TODO join with MAJIC
     patches_expl['SURF_PARC'] = patches_expl.area
     return [patches_expl, municipalities]
@@ -114,7 +128,7 @@ def build_PAT_municipalities(pat_municipalities_dir, bdv_filename, bdv_fixes_fil
     del municipalities['Code Insee']
     # load EPCI
     admin_express_epci = gpd.GeoDataFrame.from_file(adminexpress_epci_filename)[['CODE_EPCI','NOM_EPCI']]
-    municipalities = pd.merge(admin_express_epci,municipalities,how='inner',left_on='CODE_EPCI',right_on='CODE_EPCI')
+    municipalities = municipalities.merge(admin_express_epci, how='inner', on='CODE_EPCI')
     return municipalities
 
 if __name__ == '__main__':
@@ -124,13 +138,15 @@ if __name__ == '__main__':
     except ImportError:
         from yaml import Loader
     resources={}
-    for k,v in yaml.load(open('resources/INDEX.yml','r'), Loader=Loader).items():
+    config = yaml.load(open('resources/INDEX.yml','r'), Loader=Loader)
+    for k,v in config.items():
         resources[v['variable']] = 'resources/'+v['file'] if 'file' in v else k
     [patches, municipalities] = build_initial_PAT(
         resources['rpg_parcelles_filename'],
         resources['rpg_ilots_filename'],
         resources['pat_cultural_classes_filename'],
         resources['majic_filename'],
+        resources['bdalti_dir'], config['BDALTI']['tiles'],
         resources['pat_municipalities_dir'],
         resources['bdv_filename'], resources['bdv_fixes'], resources['adminexpress_com_filename'],
         resources['adminexpress_epci_filename'])

pip-requires.txt

@@ -5,3 +5,4 @@ pandas
 tables #PyTables for hdf5 store
 tqdm
 overrides
+rasterio

resources/INDEX.yml

 # Datasources in this directory
+BDALTI:
+  description: BD Alti 75m (format ASC) MNT elevation
+  file: dl/BDALTIV2_2-0_75M_ASC_LAMB93-IGN69_FRANCE_2018-01-15/BDALTIV2/1_DONNEES_LIVRAISON_2018-01-00245/BDALTIV2_MNT_75M_ASC_LAMB93_IGN69_FRANCE/
+  variable: bdalti_dir
+  tiles:
+    - BDALTIV2_75M_FXX_0750_6600_MNT_LAMB93_IGN69
+    - BDALTIV2_75M_FXX_0675_6525_MNT_LAMB93_IGN69
+    - BDALTIV2_75M_FXX_0750_6525_MNT_LAMB93_IGN69
+    - BDALTIV2_75M_FXX_0675_6600_MNT_LAMB93_IGN69
+    - BDALTIV2_75M_FXX_0675_6450_MNT_LAMB93_IGN69
+    - BDALTIV2_75M_FXX_0750_6450_MNT_LAMB93_IGN69
+  webpage: http://professionnels.ign.fr/bdalti
+  url: https://wxs-telechargement.ign.fr/jvam1hsjm11u8voorw81v2xb/telechargement/prepackage/BDALTI-75M_PACK_FXX_2018-01-24$BDALTIV2_2-0_75M_ASC_LAMB93-IGN69_FRANCE_2018-01-15/file/BDALTIV2_2-0_75M_ASC_LAMB93-IGN69_FRANCE_2018-01-15.7z
 
 BV2012_au_01-01-2018_V2:
   file: dl/BV2012_au_01-01-2018_V2/BV2012 au 01-01-2018_V2.xls