Suite

Script Python autonome pour diviser une polyligne avec une couche de points


J'ai une couche de polylignes et une couche de points qui se trouvent au-dessus des polylignes. Je n'ai pas la licence avancée ArcGis, je ne peux donc pas utiliser l'outil Split Line at Point. Quelqu'un peut-il m'aider avec un script Python autonome qui divisera ma polyligne par ma couche de points ? (Au fait, je sais qu'il existe un outil comme celui-ci dans ET Geowizards, mais je n'ai pas la licence pour cela, donc je ne peux pas faire de cet outil un script autonome).


Voici un script qui pourrait le faire pour vous. La classe native arcpy.Geometry a une méthode appelée "cut" qui coupera n'importe quelle entité à l'aide d'une autre polyligne. Malheureusement, puisque vous utilisez des "points", nous devons créer de "fausses" lignes à partir de ces points. J'ai essentiellement fait une polyligne à gratter avec les points[(Point.X+10, Point.Y+10), (Point.X-10, Point.Y-10)]- par exemple. une ligne diagonale pour chaque point pour effectuer chaque coupe. Lorsque vous avez plus d'un point coupant une ligne, vous devez recouper les lignes de manière récursive de plus en plus petites jusqu'à ce qu'elles ne puissent plus être coupées (sans produire de géométries de longueur nulle).

Quelques mises en garde que je dois mentionner :

  • Le script suppose qu'il n'y a pas de géométries en plusieurs parties dans les deux classes d'entités. j'ai toujours utiliséPoint.premierPointpour chaque point, de sorte que tous les multi-points supplémentaires par la suite seront omis. Je n'ai pas testé l'effet de la coupe d'une ligne en plusieurs parties, mais je suppose que les descendants resteront également en plusieurs parties.
  • Lorsque vous avez un segment parfaitement diagonal d'une polyligne, il ne sera pas coupé car notre "fausse ligne de coupe" est également diagonale. Vous devrez peut-être les diviser manuellement, s'ils existent.
  • La sortie n'aura pas les attributs des polylignes d'origine, mais une jointure spatiale peut s'en occuper.
  • Enfin, bien sûr, assurez-vous que les deux classes d'entités sont dans le même système de coordonnées.

Il suffit de faire l'outil squelette dans la boite à outils qui accepte 3 paramètres :

  1. Classe d'entités polyligne (Direction = Entrée, Type = Classe d'entités, Filtre = Ligne)
  2. Classe d'objets point (Direction = Entrée, Type = Classe d'objets, Filtre = Point)
  3. Classe d'entités en sortie (Direction = Sortie, Type = Classe d'entités)

Voici quelques exemples:

Et le code :

__author__ = "John K. Tran" __contact__ = "[email protected]" __version__ = "4.0" __created__ = "7/16/15" __credits__ = "http://gis.stackexchange.com/questions/154708/ standalone-python-script-to-split-a-polyline-with-a-point-layer" """ 'Cuts', 'splits' or 'dices' une classe d'entités polyligne à l'aide d'une classe d'entités ponctuelles. Le schéma de sortie ( c'est-à-dire la liste des champs) sera vide, mais vous pouvez utiliser une jointure spatiale pour repeupler les attributs dans la classe d'entités en sortie. """ import arcpy import os import sys arcpy.env.overwriteOutput = True arcpy.SetProgressor("default", "Firing up script… ") linefc = arcpy.GetParameterAsText(0) # Votre classe d'entités polyligne pointfc = arcpy.GetParameterAsText(1) # Votre classe d'entités ponctuelles output = arcpy.GetParameterAsText(2) # Nouvelle classe d'entités en sortie. spatialref1 = arcpy.Describe(pointfc).spatialReference spatialref2 = arcpy.Describe(linefc).spatialReference assert spatialref1.name == spatialref2.name, "Assurez-vous que les deux classes d'entités ont le même système de coordonnées projetées" # Si la licence Advanced est disponible, utilisez simplement les outils de géotraitement normaux. if arcpy.ProductInfo() dans [u'ArcInfo', u'ArcServer'] : arcpy.SetProgressorLabel("Splitting lines at points") arcpy.SplitLineAtPoint_management(linefc, pointfc, output, 1.0) arcpy.SetProgressorLabel("Suppression des tranches en double ") outshapefieldname = arcpy.Describe(output).shapeFieldName arcpy.DeleteIdentical_management(output, [outshapefieldname]) sys.exit(0) # Sinon, continuez avec le script. # Configurez quelques variables préliminaires pour commencer. arcpy.SetProgressorLabel("Regrouper les géométries") points = [row[0] pour la ligne dans arcpy.da.SearchCursor(pointfc, "[email protected]")] lines = [row[0] pour la ligne dans arcpy.da.SearchCursor(linefc , "[email protected]")] cutlines = list() iterations = 0 # Définir une fonction pour la coupe nous permettra de la réutiliser pour les coupes suivantes. def CutLines() : arcpy.SetProgressorLabel("Running cut : Iteration {0}".format(str(iterations))) for line in lines[:] : iscut = "Invalid" if line.length > 0.0 : # Assurez-vous ce n'est pas une géométrie vide. iscut = "Not Cut" pour le point dans les points : if line.distanceTo(point) < 1.0 : # Même les points "coïncidants" peuvent apparaître comme spatialement non coïncidants dans leurs valeurs XY à virgule flottante, nous définissons donc une tolérance. snappoint = line.snapToLine(point).firstPoint # Pour assurer une coïncidence, accrochez le point à la ligne avant de continuer. sinon (snappoint.equals(line.lastPoint) et snappoint.equals(line.firstPoint)): # Assurez-vous que le point n'est pas sur une extrémité de ligne, sinon la coupe produira une géométrie vide. cutline1, cutline2 = line.cut(arcpy.Polyline(arcpy.Array([arcpy.Point(snappoint.X+10.0, snappoint.Y+10.0), arcpy.Point(snappoint.X-10.0, snappoint.Y-10.0) ]), spatialref1)) # Couper la ligne. if cutline1.length > 0.0 et cutline2.length > 0.0 : # Assurez-vous que les deux descendants ont une géométrie non nulle. lines.insert(0, cutline1) # Renvoyer les lignes de coupe dans la liste des "lignes" à recouper. lines.insert(0, cutline2) # La boucle coupée ne sortira que lorsque toutes les lignes ne pourront pas être recoupées de plus en plus petites (sans produire des géométries de longueur nulle). iscut = "Cut" if iscut == "Not Cut": cutlines.insert(0, line) lines.remove(line) # Effectue la boucle Cut pseudo-récurive jusqu'à ce que la polyligne FC ne puisse plus être coupée. while lines : CutLines() itérations += 1 if iterations > 500 : # Fail-safe pour arrêter une boucle infinie si quelque chose ne va pas (ou si plus de 500 points coupent une seule ligne). break # Crée la classe d'entités en sortie. arcpy.SetProgressorLabel("Création d'une classe d'entités en sortie") arcpy.CreateFeatureclass_management(os.path.dirname(output), os.path.basename(output), "POLYLINE", None, None, None, spatialref2) # Insérez chaque ligne de coupe dans la nouvelle classe d'entités. arcpy.SetProgressorLabel("Insertion de lignes de coupe") avec arcpy.da.InsertCursor(output, "[email protected]") comme insertcursor : pour la ligne de coupe dans les lignes de coupe : insertcursor.insertRow((cutline,)) # Supprimez les doublons en comparant leurs objets géométriques. arcpy.SetProgressorLabel("Suppression des doublons") crows = [ligne pour ligne dans arcpy.da.SearchCursor(output, ["[email protected]", "[email protected]"])] deleteOIDs = set() # Conserver les OID pour suppression, effectué dans le curseur suivant. pour drow in crows : # Comparez chaque géométrie à toutes les autres géométries dans le FC de sortie. dOID = drow[0] dgeom = drow[1] si dOID pas dans deleteOIDs : pour srow dans corbeaux : sOID = srow[0] sgeom = srow[1] si sOID != dOID et sOID pas dans deleteOIDs : # Empêcher les comparaisons redondantes , ne comparez que les nouvelles combinaisons. dextent = dgeom.extent sextent = sgeom.extent if dextent.XMin == sextent.XMin et dextent.XMax == sextent.XMax et dextent.YMin == sextent.YMin et dextent.YMax == sextent.YMax : # Voir si leurs étendues correspondent en premier. if dgeom.equals(sgeom): # Ensuite, voyez s'ils sont égaux. deleteOIDs.add(sOID) # Si c'est le cas, ajoutez l'OID à la liste des OID à supprimer. avec arcpy.da.UpdateCursor(output, "[email protected]") comme deletecursor : # Supprimez les enregistrements ici en utilisant votre liste de "DeleteOIDs". pour delrow dans deletecursor : delOID = delrow[0] si delOID dans deleteOIDs : deletecursor.deleteRow() # Nettoyer les variables restantes. namespace = dir() pour var dans ["linecursor", "pointcursor", "insertcursor", "deletecursor"] : si var dans l'espace de noms : exec "del " + var arcpy.ResetProgressor()

Si vous obtenez des erreurs, faites le moi savoir.


Voir la vidéo: Exemple de script Python sous QGis calculant la position de Chacun des Pixels dun Raster 33 (Octobre 2021).