Recherche d'un processus d'historisation de base de données d'occupation des sols appliqué au référentiel géographique forestier de l'IGN

VI.A - « identite_1.py »

# User : Romain Louvet, stagiaire

# 28/08/2013

# IGN, equipe produit forêt environnement, projet OCS & mise a jour du RGFor

# DESCRIPTION :

# test d'algorithme d'appariement attribuant automatiquement une identite pre-existante a un

nouvel objet.

# Les criteres d'appariement sont :

# - CRITERE 1 : le successeur possede la meme nomenclature que l'anterieur

# - CRITERE 2 : le successeur est compris dans l'anterieur, le successeur recupere l'identite de

l'anterieur

# - CRITERE 3 : le successeur contient l'anterieur et l'aire de l'anterieur est > ou = a 40% de l'aire du

successeur

# - CRITERE 4 : le successeur intersecte l'anterieur et l'aire du successeur est > ou = a 40% de l'aire du

anterieur

# Si CRITERE 1 ET (CRITERE 2 ou 3 ou 4) sont vrais, l'entite successeur = l'entite anterieur

# Sinon, le successeur est une nouvelle entite.

# IMPORTANT :

# - ouvrir session de mise a jour avant de charger

# - charger dans fenetre python d'ArcGis sous arcmap (mode immediat d'utilisation de python avec

ArcGis)

# - charger indentite_1, activer mode edition, charger intendite_2

# - si erreur avec 'locks' au moment de sauvegarder, fermer arcmap, redemarrer, effacer feature

class cree, relancer OU

# ouvrir dossier de la gdb rechercher "lock" et supprimer lock sur "APRES..." et "AVANT..." (solution

non teste)

import arcpy

from datetime import datetime import sys

### I - VARIABLES a definir

fc = " C:/Users/Romain/Desktop/IGN/test_histoRGFOR_OCS/test_rgfor65.gdb/test65/RGFOR65_test"

# adresse de la table d'actualites

date = """ date '2010-08-01' """# date de la PVA utilisee pour la mise a jour

date_fmt = '%Y-%m-%d'

gdb = " C:/Users/Romain/Desktop/IGN/test_histoRGFOR_OCS/test_rgfor65.gdb/test65/" # adresse

d'enregistrement

AVANT = "AVANT2010" # nom de la table des objets valide jusqu'a la date de PVA

APRES = "APRES2010" # nom de la table des objets valide depuis la date de PVA

### II - creation couche AVANT - APRES date de la PVA a partir d'une selection

arcpy.MakeFeatureLayer_management (fc, "actulyr")

x

arcpy.SelectLayerByAttribute_management ("actulyr", "NEW_SELECTION", """ "DATE_V_TO" = """+date)

arcpy.CopyFeatures_management ("actulyr",gdb+AVANT)

arcpy.SelectLayerByAttribute_management ("actulyr", "NEW_SELECTION", """ "DATE_V_FROM" =

"""+date)

arcpy.CopyFeatures_management ("actulyr",gdb+APRES)

identite_2.py

identite_3.py

evenements_test.py

xi

VI.B - « identite_2.py » import arcpy

from datetime import datetime import sys

### I - VARIABLES a definir

fc = " C:/Users/Romain/Desktop/IGN/test_histoRGFOR_OCS/test_rgfor65.gdb/test65/RGFOR65_test"

# adresse de la table d'actualites

date = '2010-08-01'# date de la PVA utilisee pour la mise a jour

date_fmt = '%Y-%m-%d'

gdb = " C:/Users/Romain/Desktop/IGN/test_histoRGFOR_OCS/test_rgfor65.gdb/test65/" # adresse

d'enregistrement

AVANT = "AVANT2010" # nom de la table des objets valide jusqu'a la date de PVA

APRES = "APRES2010" # nom de la table des objets valide depuis la date de PVA

# liste des codes d'entites

code_ent = open(" C:/Users/Romain/Desktop/IGN/test_histoRGFOR_OCS/code_RGFOR065.txt","r") code_ent = code_ent.read()

code_ent = code_ent.split("\n")

### III - selection par localisation dans AVANT depuis APRES

# creation d'une liste contenant les identifiants absolus des lignes listeOID = []

# creation d'un curseur de recherche pour parcourir la table "APRES" lignes = arcpy.SearchCursor(APRES)

# boucle recherche ligne par ligne de la table "APRES" et ajout des valeurs a la liste "listeOID" for ligne in lignes:

# recupere la valeur du champ OBJECTID, champ d'identifiant absolu

value = ligne.getValue("OBJECTID")

# ajoute a la liste "listeOID" listeOID.append(value)

# fait de meme avec "CLEABS" value = ligne.getValue("CLEABS")

# mesure "listeOID" pour determiner le nombre d'iteration lgr = len(listeOID)

# trouve la derniere valeur "ENTITE_ID"

# servira a appeler nouveau code pour nouvelle entite

arcpy.SelectLayerByAttribute_management ("actulyr", "NEW_SELECTION", """ "ENTITE_ID" = (SELECT

MAX ("ENTITE_ID") FROM RGFOR65_test) """)

resultat = arcpy.SearchCursor("actulyr")

for cmax in resultat:

cle = cmax.getValue("OBJECTID")

xii

# boucle

for i in range(lgr):

# selection par attribut dans la table APRES de la ligne d'index "i" dans "listeOID"

OID = str(listeOID[i])

select = """ "OBJECTID" = """+OID

arcpy.SelectLayerByAttribute_management (APRES, "NEW_SELECTION", select)

# creation d'un curseur de recherche dans la table "APRES" lignes = arcpy.SearchCursor(APRES)

# condition si l'objet est une nouvelle entite new = 0

# boucle de recherche ligne par lignee de la table "APRES" des valeurs "NOMEN", "SHAPE_Area" et

"DATE_V_CREA"

for ligne in lignes:

tfvAP = ligne.getValue("NOMEN")

aireAP = ligne.getValue("SHAPE_Area")

# CRITERES 3 et 4 : calcul du % de l'aire dans "APRES" aireAP100a = aireAP*40/100

aireAP100b = aireAP/1.4

# CRITERE 2 : selection par emplacement dans la couche "AVANT" des polygones # qui contiennent polygones dans la couche "APRES"

arcpy.SelectLayerByLocation_management (AVANT, "CONTAINS", APRES)

# creation d'un curseur pour la couche "AVANT" lignes2 = arcpy.SearchCursor(AVANT)

# condition si selection ne fonctionne pas entIDAV = 0

# pour toutes les lignes de la selection for ligne2 in lignes2:

# recupere valeur du champ "NOMEN" et "DATE_CREA"

tfvAV = ligne2.getValue("NOMEN")

# CRITERE 1 : si valeur "NOMEN" est identique a celle de l'objet dans "APRES"

if tfvAP == tfvAV:

# recupere valeur du champ "ENTITE_ID", "PARTIE" et "DATE_V_CREA" de "AVANT" et

l'attribue a celui de "APRES"

entIDAV = ligne2.getValue("ENTITE_ID")

entIDAV = str(entIDAV)

entiteID = '"'+""+entIDAV+""+'"'

partieAV = ligne2.getValue("PARTIE") partieAV = str(partieAV)

partie = '"'+""+partieAV+""+'"'

crea = ligne2.getValue("DATE_V_CREA")

xiii

crea = '"'+""+str(crea)+""+'"'

arcpy.CalculateField_management (APRES,"ENTITE_IDtest", entiteID, "PYTHON") arcpy.CalculateField_management (APRES,"PARTIE_test", partie, "PYTHON") arcpy.CalculateField_management (APRES,"DATE_V_CREA_test", crea, "PYTHON")

# l'objet n'est pas une nouvelle entite new = 1

# si CONTAINS ne donne rien : if entIDAV == 0:

# CRITERE 3 : selection par emplacement dans la couche "AVANT" des polygones # qui sont contenus par des polygones dans la couche "APRES" arcpy.SelectLayerByLocation_management (AVANT, "WITHIN", APRES)

# creation d'un curseur pour la couche "AVANT" lignes2 = arcpy.SearchCursor(AVANT)

# condition si selection ne fonctionne pas entIDAV = 0

# listes pour la selection liAIREAV = [] liENTITEAV = [] liPARTIE = [] liCREA = []

# pour toutes les lignes dans la selection for ligne2 in lignes2:

# recupere valeur du champ "NOMEN" tfvAV = ligne2.getValue("NOMEN")

# CRITERE 1 : si valeur "NOMEN" de "AVANT" est identique a celle de l'objet

# dans "APRES" ajoute valeur de son champ "SHAPE_Area" et "ENTITE_ID" a une liste

if tfvAP == tfvAV:

aireAV = ligne2.getValue("SHAPE_Area")

entIDAV = ligne2.getValue("ENTITE_ID")

partieAV = ligne2.getValue("PARTIE")

crea = ligne2.getValue("DATE_V_CREA")

liAIREAV.append(aireAV) liENTITEAV.append(entIDAV) liPARTIE.append(partieAV) liCREA.append(crea)

# trouve l'objet le plus grand et son "ENTITE_ID" try:

airemax = 0

airemax = max(liAIREAV)

xiv

ind = liAIREAV.index(airemax)

except (SyntaxError,ValueError,IndexError): pass

# si aire "AVANT" remplie condition par rapport a aire "APRES" if airemax > aireAP100a:

# recupere valeur du champ "ENTITE_ID", "PARTIE" et "DATE_V_CREA" de "AVANT" et

l'attribue a celui de "APRES"

entIDAV = liENTITEAV[ind]

partieAV = liPARTIE[ind]

crea = liCREA[ind]

entIDAV = str(entIDAV)

entiteID = '"'+""+entIDAV+""+'"' partieAV = str(partieAV)

partie = '"'+""+partieAV+""+'"' crea = '"'+""+str(crea)+""+'"'

arcpy.CalculateField_management (APRES,"ENTITE_IDtest", entiteID, "PYTHON") arcpy.CalculateField_management (APRES,"PARTIE_test", partie, "PYTHON") arcpy.CalculateField_management (APRES,"DATE_V_CREA_test", crea, "PYTHON")

# l'objet n'est pas une nouvelle entite new = 1

# si WITHIN ne donne rien non plus if entIDAV == 0:

# CRITERE 4 : selection par emplacement dans la couche "AVANT" des polygones # qui possedent une intersection avec les polygones dans la couche "APRES" arcpy.SelectLayerByLocation_management (AVANT, "INTERSECT", APRES)

# creation d'un curseur pour la couche "AVANT" lignes2 = arcpy.SearchCursor(AVANT)

# listes pour la selection liAIREAV = [] liENTITEAV = [] liPARTIE = [] liCREA = []

for ligne2 in lignes2:

# recupere valeur du champ "NOMEN" tfvAV = ligne2.getValue("NOMEN")

# CRITERE 1 : si valeur "TFV" de "AVANT" est identique a celle de l'objet

# dans "APRES" ajoute valeur de son champ "SHAPE_Area" et "ENTITE_ID" a une liste

if tfvAP == tfvAV:

aireAV = ligne2.getValue("SHAPE_Area")

entIDAV = ligne2.getValue("ENTITE_ID")

xv

partieAV = ligne2.getValue("PARTIE") crea = ligne2.getValue("DATE_V_CREA")

liAIREAV.append(aireAV) liENTITEAV.append(entIDAV) liPARTIE.append(partieAV) liCREA.append(crea)

# trouve l'objet le plus grand et son "ENTITE_ID"

try:

airemax = 0

airemax = max(liAIREAV)

ind = liAIREAV.index(airemax)

except (SyntaxError,ValueError,IndexError):

pass

# si aire "AVANT" remplie condition par rapport a aire "APRES" if airemax >= aireAP100b:

# recupere valeur du champ "ENTITE_ID", "PARTIE" et "DATE_V_CREA" de "AVANT" et

l'attribue a celui de "APRES"

entIDAV = liENTITEAV[ind]

partieAV = liPARTIE[ind]

crea = liCREA[ind]

entIDAV = str(entIDAV)

entiteID = '"'+""+entIDAV+""+'"' partieAV = str(partieAV)

partie = '"'+""+partieAV+""+'"' crea = '"'+""+str(crea)+""+'"'

arcpy.CalculateField_management (APRES,"ENTITE_IDtest", entiteID, "PYTHON") arcpy.CalculateField_management (APRES,"PARTIE_test", partie, "PYTHON") arcpy.CalculateField_management (APRES,"DATE_V_CREA_test", crea, "PYTHON")

# l'objet n'est pas une nouvelle entite new = 1

# si l'objet est une nouvelle entite if new == 0:

# cherche un nouveau code ENTITE_ID dans la liste entiteID = '"'+""+str(code_ent[cle])+""+'"'

# attribue premier element de la liste de code pour PARTIE partie = '"'+""+"aaa"+""+'"'

# attribue nouvelle date de creation crea = datetime.strptime(date,date_fmt) crea = '"'+""+str(crea)+""+'"'

xvi

arcpy.CalculateField_management (APRES,"ENTITE_IDtest", entiteID, "PYTHON") arcpy.CalculateField_management (APRES,"PARTIE_test", partie, "PYTHON") arcpy.CalculateField_management (APRES,"DATE_V_CREA_test", crea, "PYTHON")

# prochain code cle = cle + 1

xvii

VI.C - « identite_3.py » import arcpy from datetime import datetime import sys

### I - VARIABLES a definir

fc = " C:/Users/Romain/Desktop/IGN/test_histoRGFOR_OCS/test_rgfor65.gdb/test65/RGFOR65_test"

# adresse de la table d'actualites

date = '2010-08-01'# date de la PVA utilisee pour la mise a jour

date_fmt = '%Y-%m-%d'

gdb = " C:/Users/Romain/Desktop/IGN/test_histoRGFOR_OCS/test_rgfor65.gdb/test65/" # adresse

d'enregistrement

AVANT = "AVANT2010" # nom de la table des objets valide jusqu'a la date de PVA

APRES = "APRES2010" # nom de la table des objets valide depuis la date de PVA

### IV- Partie

### Si plusieurs objets ont le meme code ENTITE_ID, ils doivent avoir differents codes PARTIE

### Variables

# liste des codes PARTIE

code_partie =

open(" C:/Users/Romain/Desktop/IGN/test_histoRGFOR_OCS/code_PARTIE065.txt","r") code_partie = code_partie.read()

code_partie = code_partie.split("\n")

# Variables pour IVA

liOID = [] # liste OBJECTID liENT = [] # liste ENTITE

liPAR = [] # liste PARTIE depart

liENPAR = [] # liste id d'entite complet

liARE = [] # liste SHAPE_area

# Variable pour IVB

liTTL = [] # liste totale index liPAR2 = [] # liste PARTIE arrivee

### IVA - Listes des valeurs des champs

# creation d'un curseur de recherche pour parcourir la table "APRES" lignes = arcpy.SearchCursor(APRES)

# boucle recherche ligne par ligne de la table "APRES"

for ligne in lignes:

# recupere la valeur du champ "OBJECTID", "ENTITE_IDtest", "PARTIE_test" et "SHAPE_Area"

vOID = ligne.getValue("OBJECTID")

vENT = ligne.getValue("ENTITE_IDtest")

vPAR = ligne.getValue("PARTIE_test")

vARE = ligne.getValue("SHAPE_Area")

# ajoute aux listes

xviii

liOID.append(vOID) liENT.append(vENT) liPAR.append(vPAR) liENPAR.append(vENT+vPAR) liARE.append(vARE)

liPAR2 = liPAR

lgr = len(liOID) # longueur de la table

### IVB - Boucle de detection de nouvelle partie

for i in range(lgr):

# creation et remise a zero des variables calculees par la boucle liIND = [] # liste index des lignes avec meme ENTITE_ID + PARTIE liARE2 = [] # valeur aire de ces lignes

liPAR1 = [] # liste PARTIE en cours

# si cette ligne n'a pas deja ete traitee

# essaye de trouver i dans liste totale index try:

test = liTTL.index(i)

# si ne fonctionne pas, lance la recherche except(ValueError):

# cherche si code ENTITE_ID plus PARTIE # existe plusieurs fois

enpar = liENPAR[i]

# pour chacun des codes ENTITE_ID + PARTIE for indx, ID in enumerate(liENPAR):

# si le code est egal au code recherche if ID == enpar:

# son index est ajoute a la liste liIND.append(indx)

# si cette liste est > 1, il existe plusieurs codes identique lgi = len(liIND)

if lgi > 1:

# avec la liste d'index je trouve la valeur max du code PARTIE et de l'airepour cette ENTITE_ID for j in liIND:

liARE2.append(liARE[j]) liPAR1.append(liPAR[j])

mxaire = max(liARE2) mxliPAR1 = max(liPAR1)

# trouve valeur PARTIE la plus grande dans liste num = code_partie.index(mxliPAR1)

xix

# pour toutes les lignes, sauf celle avec la plus grande aire

# le champ PARTIE possede un nouveau code

for k in liIND:

aire = liARE[k]

if aire < mxaire:

# donne index du code suivant dans la liste de code PARTIE

num = num+1

# remplace ancien par nouveau code liPAR2[k]= code_partie[num]

# ajoute les lignes deja faites liTTL.extend(liIND)

### IVC - boucle calcul du champ PARTIE

for i in range(lgr):

arr = liPAR2[i]

# selection de la ligne a modifier

OID = str(liOID[i])

select = """ "OBJECTID" = """+OID

arcpy.SelectLayerByAttribute_management (APRES, "NEW_SELECTION", select)

lignes = arcpy.SearchCursor(APRES)

# modification du champ "PARTIE" avec le nouveau code

arr = '"'+""+str(arr)+""+'"'

arcpy.CalculateField_management (APRES,"PARTIE_test", arr, "PYTHON")

xx