topic Re: faster script for making polygon from points in Python Questions

faster script for making polygon from points

LindseyWood — Tue, 12 Mar 2013 16:29:00 GMT

I am just curious if any one has any ideas on how this may perform faster as some of these tables can get quite large and at the rate I am going it is not looking good: about 700 recs/hr

1) I have an input .csv file that has four xy or lat lon points per record (Make a copy into gdb for joining to shapefile after conversion)
2) Open and plot 4 points into a feature class
3) Create Min bounding geometry for shapefile
4) append to shapefile so all record geometries in the end are stored into one shapefile.
5) join attributes from table in gdb back to shapefile by oids.

Any input appreciated... originally this was an "ave" script and it is much faster so I am going to take a look at that next for ideas.
Thanks!

import sys, os, csv, datetime, time, string, traceback, arcpy from os.path import splitdrive scriptPath = sys.path[0] (drive,tail) = splitdrive(scriptPath)  # arcpy environments and liscences from arcpy import env from arcpy.sa import * arcpy.CheckOutExtension("spatial") arcpy.CheckOutExtension("management") arcpy.env.overwriteOutput = True   # Set the XY Domain to  #   xmin of -180 #   ymin of -90 #   xmax of 180 #   ymax of 90 arcpy.env.XYDomain ="-180 -90 180 90" arcpy.SpatialReference= spRef = "D:/ArcGIS/Desktop10.0/Coordinate Systems/" +\                                     "Geographic Coordinate Systems/World/WGS 1984.prj"   if __name__ == '__main__':          ##Steps to create polygon shapefile from 4 or 2 points hm....     ## 1) get list of 4 points and store in an array by x,y     print drive, tail     arcpy.env.workspace = scriptPath          outFolder = scriptPath     gdb = os.path.join(scriptPath,"temp.gdb")     pointFC = "sample1Test.shp"     csvFile = os.path.join(scriptPath,'test.csv')          #arcpy.CreateFeatureclass_management(outFolder, pointFC, "POINT", "","","","")     #Input table to database     arcpy.TableToTable_conversion(csvFile, gdb, "vegtable")       fieldNames = ['"UL_LON"','"UL_LAT"','"LL_LON"','"LL_LAT"','"UR_LON"','"UR_LAT"','"LR_LON"','"LR_LAT"','"PHOTO_ID"']          print fieldNames[0]     print csvFile     inFile = open(csvFile, "r")     #inFile = open("D:/CreateShape/temp.gdb/vegtable", "r")     headerLine = inFile.readline()     valueList = headerLine.strip().split(",")     print valueList     latIndex = valueList.index(fieldNames[0])     lonIndex = valueList.index(fieldNames[1])     llLatIndex = valueList.index(fieldNames[2])     llLonIndex = valueList.index(fieldNames[3])     urLatIndex = valueList.index(fieldNames[4])     urLonIndex = valueList.index(fieldNames[5])     lrLatIndex = valueList.index(fieldNames[6])     lrLonIndex = valueList.index(fieldNames[7])     photoId = fieldNames[8]     feats = arcpy.ListFeatureClasses()     try:         for feat in feats:             print feat             arcpy.Delete_management(feat)         arcpy.env.workspace = gdb         for feat in feats:             print feat             arcpy.Delete_management(feat)     except:         print "Couldnt be deleted"               #Read each line in the csv file     starttime = time.time()     localtime = time.asctime( time.localtime(time.time()) )     i = 0     #cursor = arcpy.InsertCursor(pointFC)     for line in inFile.readlines():         #print line         fields = line.split(",")         #print fields[0]                  points = "points%s.shp" %(i)         print points         arcpy.CreateFeatureclass_management(outFolder, points, "POINT", "","","","")         #arcpy.MakeFeatureLayer_management()         cursor = arcpy.InsertCursor(points)                  field = line.split(",")         ## 2) create points and append or merge points together in one feature class         latValue = field[latIndex]         lonValue = field[lonIndex]         #print "first point: (%s,%s)" %(latValue,lonValue)         point = arcpy.CreateObject("Point")         point.X = latValue         point.Y = lonValue         feature = cursor.newRow()         feature.shape = point         cursor.insertRow(feature)          lat2Val=field[llLatIndex]         lon2Val = field[llLonIndex]         #print "Second point: (%s)" %lat2Val,lon2Val         point = arcpy.CreateObject("Point")         point.X = lat2Val         point.Y = lon2Val         feature = cursor.newRow()         feature.shape = point         cursor.insertRow(feature)                  lat3Val=field[urLatIndex]         lon3Val = field[urLonIndex]         #print "Third point: (%s)" %lat3Val,lon3Val         point = arcpy.CreateObject("Point")         point.X = lat3Val         point.Y = lon3Val         feature = cursor.newRow()         feature.shape = point         cursor.insertRow(feature)          lat4Val=field[lrLatIndex]         lon4Val = field[lrLonIndex]         #print "Fourth point: (%s)" %lat4Val,lon4Val         point = arcpy.CreateObject("Point")         point.X = lat4Val         point.Y = lon4Val         feature = cursor.newRow()         feature.shape = point         cursor.insertRow(feature)         shpGdb = "%s/temp.gdb/box%s" %(scriptPath, i)         allMerged = "%s/temp.gdb/merged2" %(scriptPath)                  ## 3) perfom minimum bounding box tool by rectangle :)         arcpy.MinimumBoundingGeometry_management(points,shpGdb,"RECTANGLE_BY_AREA", "ALL")         #arcpy.AddField_management(shpGdb, str(photoId), "TEXT", "", "", 25)         #expression = ("\"%s\"") %(field[0])         #arcpy.CalculateField_management(shpGdb, str(photoId), expression, "PYTHON")                  if i == 0:                          arcpy.CopyFeatures_management(shpGdb, allMerged)                                  if i >0:             try:                 arcpy.Append_management(shpGdb,allMerged)                 print "appended?"             except:                 print arcpy.GetMessages()         arcpy.Delete_management(shpGdb)         arcpy.Delete_management(points)                                           i+=1      #arcpy.AddJoin_management(allMerged,"OBJECTID","vegetable", "OBJECTID")     arcpy.JoinField_management(allMerged, "OBJECTID", "vegtable", "OBJECTID")        inFile.close()     endtime = time.time()     totaltime = endtime-starttime     print "\nScript took " + str(totaltime/60) + " minutes to run"

Re: faster script for making polygon from points

MathewCoyle — Tue, 12 Mar 2013 16:39:05 GMT

If you are just creating a polygon from your coordinates out of your CSV there is no need to create point objects to then create a polygon from, you can add them to an array and input them all at once skipping the need for the minimum bounding geometry tool completely. Unless there is some other reason you are going this route.

http://resources.arcgis.com/en/help/main/10.1/index.html#//018z00000061000000

Re: faster script for making polygon from points

LindseyWood — Tue, 12 Mar 2013 17:25:10 GMT

Thanks much this is perfect! I just had not ever had to do this before so I plotted the points first.

Re: faster script for making polygon from points

LindseyWood — Sat, 11 Dec 2021 19:49:10 GMT

The new script is very fast indeed 🙂
However it makes bow ties instead of rectangles... I assume its because of the order that the lat and longs are coming in but I am not sure how to sort it to make it draw a rectangle? I have a pic that the purple is from my first script and green is the bowtie.
Thanks again!
[ATTACH=CONFIG]22555[/ATTACH]

    coordList = []
    point = arcpy.Point()
    array = arcpy.Array()
    featureList = []
        

    #Read each line in the csv file
    starttime = time.time()
    localtime = time.asctime( time.localtime(time.time()) )
    ##i = 0
    #cursor = arcpy.InsertCursor(pointFC)
    for line in inFile.readlines():
        #print line
        field = line.split(",")
        #print fields[0]
        points = [[field[latIndex],field[lonIndex]],[field[llLatIndex],field[llLonIndex]],[field[urLatIndex],field[urLonIndex]],[field[lrLatIndex],field[lrLonIndex]]] 
        coordList.append (points)
        inFile.close()
        
    

    for feature in coordList:
        # For each coordinate pair, set the x,y properties and add to the 
        #  Array object.
        print feature
        #
        for coordPair in feature:
            print coordPair[0],coordPair[1]
            point.X = coordPair[0]
            point.Y = coordPair[1]
            array.add(point)

        # Add the first point of the array in to close off the polygon
        #
        array.add(array.getObject(0))

        # Create a Polygon object based on the array of points
        #
        polygon = arcpy.Polygon(array)

        # Clear the array for future use
        #
        array.removeAll()

        # Append to the list of Polygon objects
        #
        featureList.append(polygon)

    # Create a copy of the Polygon objects, by using featureList as input to 
    #  the CopyFeatures tool.
    #
    arcpy.CopyFeatures_management(featureList, "polygons3.shp")
    
    endtime = time.time()
    totaltime = endtime-starttime
    print "\nScript took " + str(totaltime/60) + " minutes to run"

Re: faster script for making polygon from points

MathewCoyle — Tue, 12 Mar 2013 19:34:47 GMT

Yes it's your order, switch the last two

[field[urLatIndex],field[urLonIndex]],[field[lrLatIndex],field[lrLonIndex]]

[field[lrLatIndex],field[lrLonIndex], [field[urLatIndex],field[urLonIndex]]]

Re: faster script for making polygon from points

LindseyWood — Tue, 12 Mar 2013 19:45:58 GMT

😮 That fixed it! Thanks so much again I can go home happy now.