zonal_stats.org 16.1 KB

Construction d’un data frame avec les stats zonales

Init

Petit bout de python vide pour initialiser la session. Ca permet d’avoir la version de python utilisée

import geopandas as gp
import pandas as pd

Simplification du fichier TERLAB

Eliminer les géométries qui ne contiennent pas d’information de rendement.

Combien de polygones en tout?

source /work/OT/theia/oso/jordi/OTB/otb_superbuild/otb_superbuild-otb_develop-Release-install/config_otb_tf_py35.sh
ogrinfo -sql "select COUNT(*) from \"SURFACES-2017-PARCELLES-GRAPHIQUES-CONSTATEES_031_20180210\"" /work/OT/theia/oso/shapes/TERLAB/SURFACES-2017-PARCELLES-GRAPHIQUES-CONSTATEES_031_20180210.shp 

Combien contiennent du rendement?

source /work/OT/theia/oso/jordi/OTB/otb_superbuild/otb_superbuild-otb_develop-Release-install/config_otb_tf_py35.sh
ogrinfo -sql "select COUNT(*) from \"SURFACES-2017-PARCELLES-GRAPHIQUES-CONSTATEES_031_20180210\" where RENDNORME > 0" /work/OT/theia/oso/shapes/TERLAB/SURFACES-2017-PARCELLES-GRAPHIQUES-CONSTATEES_031_20180210.shp 

Eliminer ceux qui n’en contiennent pas et créer un nouveau fichier

source /work/OT/theia/oso/jordi/OTB/otb_superbuild/otb_superbuild-otb_develop-Release-install/config_otb_tf_py35.sh
ogr2ogr -sql "select * from \"SURFACES-2017-PARCELLES-GRAPHIQUES-CONSTATEES_031_20180210\" where RENDNORME > 0" -f SQLite /work/OT/theia/oso/shapes/TERLAB/terlab31_rend.shp /work/OT/theia/oso/shapes/TERLAB/SURFACES-2017-PARCELLES-GRAPHIQUES-CONSTATEES_031_20180210.shp 
source /work/OT/theia/oso/jordi/OTB/otb_superbuild/otb_superbuild-otb_develop-Release-install/config_otb_tf_py35.sh
ogrinfo -sql "select COUNT(*) from \"surfaces_2017_parcelles_graphiques_constatees_031_20180210\" where RENDNORME > 0" /work/OT/theia/oso/shapes/TERLAB/terlab31_rend.shp 

Il faut faire attention au fait que les fichiers shx et prj n’existent pas pour le fichier qui ne contient que les polys de renement. Peut-être qu’il faut donner la projection pour le fichier généré par ogr2ogr?

Génération des fichiers de stats zonales

On pourrait commencer par ajouter un identifiant unique à chaque polygone du fichier TERLAB de façon à pouvoir faire le lien entre les différents fichiers de stats locales (reflectances, DEM, masques) générés à partir du même fichier TERLAB.

/work/OT/theia/oso/jordi/src/maassp/jobs/addid.pbs

[2019-02-05 Tue] le job add id tourne Attention : j’ai dû modifier le script dans vectortools. Il vaudra mieux le copier ici.

Le code est ici /work/OT/theia/oso/jordi/src/maassp/jobs/zonal_stats.pbs

A voir si cette partie sera nécessaire une fois que les polygones sans rendement sont éliminés.

Comme l’image pour les stats zonales contient 720 bandes, le fichier vecteur contient > 720 * 4 attributs, ce qui dépasse le nb de colonnes dans un fichier sqlite. Le shapefile n’a pas ce problème, mais est limité à 4GB, ce qui est trop petit.

J’ai ensuite séparé la série en 2 demi-séries pour pouvoir générer 2 shapefiles, mais on a ensuite le pb de faire le lien entre les polygones des 2 shapefiles générés, parce qu’on perd l’identifiant.

Import des fichiers dans des DF

Le fichier TERLAB qui contient le dept 31 est:

terlabdir = '/work/OT/theia/oso/shapes/TERLAB/'
terlab31 = terlabdir+'terlab31.shp'
df_terlab = gp.read_file(terlab31)
print(df_terlab.head())

Les fichiers shape avec les stats sont

zs_dir='/work/OT/theia/oso/jordi/maassp_data'
zs_rf1=zs_dir+'/zonal_stats_refl_1.shp'
zs_rf2=zs_dir+'/zonal_stats_refl_2.shp'
zs_dem=zs_dir+'/zonal_stats_dem.shp'
zs_mask=zs_dir+'/zonal_stats_mask.shp'

Lire les fichiers dans un DF pandas

df_rf1 = gp.read_file(zs_rf1)
df_rf2 = gp.read_file(zs_rf2)
df_dem = gp.read_file(zs_dem).to_crs
df_mask = gp.read_file(zs_mask)
print(len(df_rf1))
print(len(df_rf2))
print(len(df_dem))
print(len(df_mask))

Il paraît que le fichier de DEM a perdu quelques polygones?

Vérification des projections :

print(df_dem.crs)
print(df_mask.crs)
print(df_rf1.crs)
print(df_rf1.crs)

Reprojection du DEM

df_dem = df_dem.to_crs(df_mask.crs)
print(df_dem.crs)

Fusion dans un seul DF

for c in df_terlab.columns:
    print(c)

for c in df_mask.columns:
    print(c)

print(df_dem.geometry.head())
print(df_mask.geometry.head())
print(df_rf1.geometry.head())
print(df_rf2.geometry.head())

On fait un spatial join.

df_dem_mask = gp.sjoin(df_mask, df_mask, op='intersects')
print(df_dem_mask.head())
print(len(df_dem_mask))

On a beaucoup trop d’éléments dans le nouveau DF, ce qui veut dire qu’il y a trop d’intersections. Si on le fait entre les 2 fichiers de réflectances, qui sont issus de la même donnée :

df_dem_mask = gp.sjoin(df_rf2, df_rf1, op='intersects')
print(df_dem_mask.head())
print(len(df_dem_mask))

Essayons de merger sur plusieurs colonnes

df_tmp = df_mask
df_dem_mask = df_tmp.merge(df_rf1, on=['CODE_CULTU', 'SURF_ADM'])
print(df_dem_mask.head())
print(len(df_dem_mask))

Templates


import matplotlib.pyplot as plt
# plt.figure(figsize=(10,5))
# df31.hist(column='RENDNORME')
# plt.savefig(matplot_lib_filename)
matplot_lib_filename