3

I want to create polygons based on two separate text files:

  • the first one (as parents) lists the Polygon ID (ID_el) and three to four Points ID (IDpt1, IDpt2, IDpt3..)
  • the second contains the Points ID and their xyz coordinates (RW, HW,HH).

I want for each ID in the first file to set an Array (3 to 4 rows) of the coordinates in the second file

enter image description here

I found different scripts to do some of the parts of the process, but I really need help to figure out how to get the Points ID and coordinates for each polygon into an Array, keeping also the attributes of the parent element (ElTyp, ID_el).

Here is the bit of code I used to extract these 2 text files:

with open(arcpy.GetParameterAsText (0),'r') as f1:
    with open(r'Nodes.txt', 'a') as f2:
        f2.write("ID RW HW HH" + "\n")
        for line in f1:
            if line.startswith("ND"):
                f2.write(line.replace("ND", "").lstrip())
with open (arcpy.GetParameterAsText (0),'r') as f1:
    with open (r'Polygonen.txt', 'a') as f3:
        f3.write("ELtyp ID_el IDpt1 IDpt2 IDpt3 IDpt4 Mat" + "\n")
        for line in f1:
            coord = " ".join(line.split() [0:-1])
            mat = " ".join(line.split() [-1:])
            if line.startswith ("E3T"):
                f3.write ((coord.lstrip())+ "  " + mat+ "\n")
            elif line.startswith ("E4Q"):
                f3.write(line.lstrip())
  • Is the format of both text files shown in the image is standard? Also, from my understanding of your question that you want to create separate polygons like first line 26748, 26286, 26747 refers to coordinates from second text file and want to keep ELtype and ID_el in respective polygon's attributes. Please clear both points. – Surya Mar 9 '15 at 10:50
  • Yes, the Format is standard: both files are outputs of previous steps in my script, the data (in a 2dm file) were extracted by starts of lines (startswith). I do want to create separate polygons for each line (ID_el), keeping the parents attributes (ElTyp, ID_el,IDpt1, IDpt2, Idpt3...) – PMessager Mar 9 '15 at 10:58
  • What would be the spatial reference? – Surya Mar 9 '15 at 11:12
  • Usually the 2dm data is not projected as such, but here they are in "DHDN_3_Degree_Gauss_Zone_4 - WKID: 31468. " – PMessager Mar 9 '15 at 11:31
  • It is doable. Give me some time I will upload a working script in evening (IST) when I reach home. – Surya Mar 9 '15 at 11:34
4

A similar approach to what has Erica published, but with more details:

  1. Export text files to a file geodatabase table (useful because you will get 0 for the vertex ID for those polygons which have only 3 points - thus no need to handle this later on);
  2. Convert those tables into dictionaries;

polygons {41880: (26287, 26286, 26748, 26747), 41879: (26748, 26286, 26747, 0)}

vertices {26747: (19, 20, 21), 26748: (12, 14, 15), 26286: (13, 16, 17), 26287: (22, 24, 25)}

  1. Building an array of polygons (each has points with XYZ values);
  2. Pre-create the polygon feature class with the Z coordinates stored (important); doable with the GP tool and adding the PolyID field.
  3. When loading polygons, important to specify the has_z attribute, otherwise it won't store Z-coordinates. You can see the Z-value for every vertex of the polygon when in the Editing session and having the Edit Sketch Properties window open (while editing the polygon with vertices shown);
  4. As the last step, use the Join Field GP tool to transfer all other attribute fields from the Polygons file GDB table (based on the PolyID). I didn't want to take them into the dictionary just to keep things more clear.

Useful Esri Help links: Polygon class; Writing geometries; da.Insert cursor

The ready-to-use code:

import arcpy
import os
folder = "C:\GIS\Temp"
os.chdir(folder)
polygonsFile = "polygons.txt"
verticesFile = "vertices.txt"

outfileGDB = r"C:\GIS\Temp\test.gdb"

if not arcpy.Exists(r"C:\GIS\Temp\test.gdb\polygons"):
    arcpy.TableToTable_conversion(in_rows=polygonsFile,
                                  out_path="C:\GIS\Temp\test.gdb",
                                  out_name="polygons",
                                  where_clause="")

    arcpy.TableToTable_conversion(in_rows=verticesFile,
                                  out_path="C:\GIS\Temp\test.gdb",
                                  out_name="vertices",
                                  where_clause="")

arcpy.env.workspace = outfileGDB
with arcpy.da.SearchCursor("polygons","*") as poly_cur:
    print "polygons"
    poly_dict = {x[1]: x[2:] for x in poly_cur}
    print poly_dict

with arcpy.da.SearchCursor("vertices","*") as vertex_cur:
    print "vertices"
    vertex_dict = {x[1]: x[2:] for x in vertex_cur}
    print vertex_dict

polyArray = {}
for polykey in poly_dict:
    polynodes = poly_dict[polykey]
    coordsList = []
    for polynode in polynodes:
        coords = (v for k,v in vertex_dict.iteritems() if k == polynode)
        for coord in coords:
            coordsList.append(coord)
    polyArray[polykey] = coordsList

print polyArray
#>>>polyArray
#>>>{41880: [(22, 24, 25), (13, 16, 17),(12, 14, 15), (19, 20, 21)],
#>>>41879: [(12, 14, 15), (13, 16, 17), (19, 20, 21)]}

features = []
fc = r"C:\GIS\Temp\test.gdb\PolygonFC"

for key,values in polyArray.iteritems():
    features.append((key,arcpy.Polygon(arcpy.Array([arcpy.Point(*value) for value in values]),
                                       arcpy.SpatialReference(4326),True))) #True - for has_z parameter

with arcpy.da.InsertCursor(fc,["PolyID","SHAPE@"]) as cur:
    for feature in features:
        cur.insertRow(feature)
  • Alex, after few minor adjustements, everything worked fine. I have not understood all of it, I need to do some more digging. Would you know how to get the WKID from a SR given thourgh GetParameter as Text? – PMessager Mar 10 '15 at 12:44
  • @PMessager, you could make comments about what is wrong/need to be adjusted, so I could edit my answer for other users. To get WKID, please look at this answer: gis.stackexchange.com/a/40959/14435. Just use the loadFromString method of the SpatialReference class. Then get the 'factoryCode` property of the class instance. – Alex Tereshenkov Mar 10 '15 at 14:48
  • After creating the dbf from the text, the number of fields increases so I adjusted the Indexes in the poly_curs part into "poly_dict = {x[2]: x[3:] for x in poly_cur}", change "PolyID" to "ID_el" and add a JoinField to update the attributes. Nothing much really. – PMessager Mar 11 '15 at 9:26
  • OK, this is great. – Alex Tereshenkov Mar 11 '15 at 9:41
2

This is an outline/pseudocode of one possible method...

  1. Create dictionaries for the point ID and coordinates. Populate as you loop through the table.

    dictRW = {}
    dictHW = {}
    dictHH = {}
    
    for row in arcpy.da.SearchCursor(tblCoordinates):
        ID = row[0]
        dictRW[ID] = row[1]
        dictHW[ID] = row[2]
        dictHH[ID] = row[3]
    
  2. Create dictionaries using the polygon ID and its associated point IDs. Populate as you loop through the table.

    listPolygons = []
    ELType = {}
    pt1 = {}
    pt2 = {}
    pt3 = {}
    pt4 = {}
    
    for row in arcpy.da.SearchCursor(tblPolygon):
        polyID = row[1]
        ELType[polyID] = row[0]
        pt1[polyID] = row[2]
        pt2[polyID] = row[3]
        pt3[polyID] = row[4]
        pt4[polyID] = row[5]
        listPolygons.append(polyID)
    
  3. Work through the list of polygons and points to insert the data into a feature class. (The following assumes the feature class exists and has the fields that you want. If you want to create it in-script, refer to arcpy.CreateFeatureclass_management)

    for polyID in listPolygons:
        featureData = []
    
        featureData.append(polyID) # first attribute is the polygon ID
    
        ptID1 = pt1[polyID]
        ptID2 = pt2[polyID]
        ptID3 = pt3[polyID]
        ptID4 = pt4[polyID]
    
        if not ptID4 is None: # this is an attempt to deal with 3-point polygons
            array = arcpy.Array([arcpy.Point(dictRW[ptID1], dictHW[ptID1], dictHH[ptID1]),
                                 arcpy.Point(dictRW[ptID2], dictHW[ptID2], dictHH[ptID2]),
                                 arcpy.Point(dictRW[ptID3], dictHW[ptID3], dictHH[ptID3]),
                                 arcpy.Point(dictRW[ptID4], dictHW[ptID4], dictHH[ptID4])])
        else:
            array = arcpy.Array([arcpy.Point(dictRW[ptID1], dictHW[ptID1], dictHH[ptID1]),
                                 arcpy.Point(dictRW[ptID2], dictHW[ptID2], dictHH[ptID2]),
                                 arcpy.Point(dictRW[ptID3], dictHW[ptID3], dictHH[ptID3])])
    
        polygon = arcpy.Polygon(array)
        featureData.append(polygon)
    
        with arcpy.da.InsertCursor(outFC, ["PolyID", "SHAPE@"]) as cursor:
            cursor.insertRow(featureData)
    

Note that the InsertCursor fields should be adapted if you want to preserve/insert additional attributes (e.g., ELType). Just make sure you also include a featureData.append() statement in order -- that is building the list of data which the insertRow command uses.

  • 1
    funny that after looking at your code I've found the PolyID, this is exactly how I am used to name this kind of field. :) – Alex Tereshenkov Mar 9 '15 at 13:19
  • Well, i got it to run after small changes (ID_el instead of PolyID, for eaxample). Many thx, Erica, now I also got more understanding on dictionnaries and array. – PMessager Mar 10 '15 at 12:47
2

As my understanding of your problem, I wrote a script below. I assumed that you have X,Y and Z values for every point.

Script:

# Example of PrentID text file
#ELtyp ID_el IDpt1 IDpt2 IDpt3
#E3T 41879 26748 26286 26747
#E3T 41880 26287 26286 26748
#E3T 41881 26749 26287 26748
#E3T 41882 26749 26750 26287

#-----------------------------------------------------------
# Example of PointsID text file
#ID RW HW HH
#26748 4.40384415e+006 5.79104121e+006 7.85830000e+001
#26749 4.40386009e+006 5.79106456e+006 7.83740000e+001
#26750 4.40386214e+006 5.79108649e+006 7.80960000e+001
#26286 4.40381493e+006 5.79104349e+006 7.85940000e+001
#26287 4.40383655e+006 5.79107131e+006 7.84350000e+001
#26747 4.40383190e+006 5.79101331e+006 7.87290000e+001

#-----------------------------------------------------------
import arcpy

arcpy.env.overwriteOutput = True
def main():
    # Polygon ID text file path
    parent = r"D:\Python\Scratch\ParentID.txt"
    # Points ID text file path
    point = r"D:\Python\Scratch\Points.txt"
    # spatial reference wkid
    wkid = 31468
    # output path (could be gdb or folder which creates
    # feature class or shapefile respectively)
    output = r"D:\Python\ScratchDatabase\Geodatabase.gdb\polygon_123"


    points = {}
    fields = []
    # create point id and coordinates dictionary
    with open(point, "rt") as point_file:
        for line in point_file:
            if line.upper().startswith('ID'.upper()):
                continue
            points.update({line.split()[0]: [float(pt) for pt in line.split()[1:]]})
    sr = arcpy.SpatialReference(wkid)
    ftr = arcpy.CreateFeatureclass_management('in_memory', 'Test', "POLYGON", spatial_reference=sr)
    flds = [fld.name for fld in arcpy.ListFields(ftr)]
    with open(parent, 'rt') as parent_file:
        for line in parent_file:
            if line.upper().startswith("ELtyp".upper()):
                for field in line.split():
                    if field not in flds:
                        # Add Fields for attributes
                        arcpy.AddField_management(ftr, field, "TEXT", field_length=10, field_alias=field)
                    fields.append(field)
                fields.append("SHAPE@")
            else:
                # insert values
                with arcpy.da.InsertCursor(ftr, fields) as cursor:
                    row = [val for val in line.split()]
                    pts = [arcpy.Point(points.get(val)[0],
                                       points.get(val)[1],
                                       points.get(val)[2]) for val in line.split()[2:]]
                    row.append(arcpy.Polygon(arcpy.Array(pts), sr, True))
                    cursor.insertRow(row)
    arcpy.CopyFeatures_management(ftr, output)
    print "Done"


if __name__ == '__main__':
    main()

May be this will help.

  • @ Surya, many thanks for your script!. I ran it and got a Runtime error on line 41 "TypeError: 'NoneType' object has no attribute 'getitem'". The polygon file is created with all required fields but without geometrie. Thanks for giving some great lines here. I have learn a lot. – PMessager Mar 10 '15 at 12:37
  • @PMessager: Thanks for feedback. That may be due to the dictionary which I made to store X,Y,Z values. I checked the script before posting that worked for me. Glad to know that I helped a little bit. :) – Surya Mar 10 '15 at 13:10
  • One last thing, I believe your work is related to TIN model. How did you generated those two text files? – Surya Mar 10 '15 at 13:23
  • I have updated the bit of code used to extract the files. It's originally a 2dm file (part of a SMS hydraulic model representing a mesh) – PMessager Mar 11 '15 at 9:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.