0

I have an Esri polygon feature class in a grid that covers a geographic region, and one field, which we'll call 'Label' has an arbitrary label that consists of two digits for north/south, and two digits for east/west location. So when completely labeled a grid may look like this:

1010 1011 1012 1013
0910 0911 0912 0913
0810 0811 0811 0813

As you move east/west, the second pair of numbers increases/decreases respectively. As you move north/south, the first pair of numbers increases/decreases respectively.

The grids I'm dealing with may only have a handful populated, i.e.:

0000 1011 0000 1013
0000 0000 0000 0000
0000 0000 0811 0000 (with the zeroes just being null)

Or it may look like this: Grid Example (However, that it is never perfectly square like this.)

I've attempted to iterate through each polygon with arcpy, but that road I went down was too memory-intensive.

Can I save the 'Labels' field out into a Pandas or geopandas dataframe and manipulate it, and import back into the same field, keeping the geographic references/locations the same? I'm at a complete loss as to how to reference north/south/east/west (or a very literal up/down/left/right) sense of direction within a dataframe. How do I attack this problem

1
  • Perhaps show screenshot.
    – FelixIP
    Commented Jun 4, 2018 at 19:40

3 Answers 3

3

You are asking two question (your question is still to broad):

  1. Populate missing fields/labels for a grid of polygons with python?
  2. Can I save the 'Labels' field out into a Pandas or geopandas dataframe and manipulate it, and import back into the same field, keeping the geographic references/locations the same?

This is an answer to 2: You can use pd.DataFrame.from_records and the da.SearchCursor to create a pandas dataframe with centroid as one column and labels as the other. Then do what you want i pandas and use a Dictionary of centroids as keys and new labels as values and da.UpdateCursor to update the feature class:

import arcpy
import pandas as pd

fc = r'C:\database.gdb\feature_class'
labelfield = 'OldLabel'
field_to_calculate = 'NewLabel'

df = pd.DataFrame.from_records(data=(i for i in arcpy.da.SearchCursor(fc, ['SHAPE@XY',labelfield])),
                                     columns=['XY',labelfield])
#do something in pandas with the labels

d = df.set_index('XY')[labelfield].to_dict()

with arcpy.da.UpdateCursor(fc,['SHAPE@XY',field_to_calculate]) as cursor:
    for row in cursor:
        row[1] = d[row[0]]
        cursor.updateRow(row)

You might need to round the coordinates, not sure if pandas will keep the decimals exactly the same.

1

I don't know a thing about Pandas, but I'd approach labeling these grids from a different direction. The direction I'd go is iterating through polygons with labels, selecting of neighbor polygons, and performing calculations based on extent coordinate distances.

I have a similar grid data set handy that I was able to test on:

enter image description here

GRIDNO is a text field that will be populated with new values.

Determining neighbor numbers will involve checking distances between each populated grid's extent and the neighbor. For example, if a neighboring polygon has a similar Y-min for its extent but a substantial X-min then the polygon is the eastern neighbor. Carry similar logic over to all sides, and then repeat the process on newly populated polygons until no new polygons get populated.

Grids representing real-world features tend to not perfectly align with X- and Y- grids or be perfectly uniform. Because of this it's useful to have a minimum dimension of a polygon. If a polygon's neighbor's extent has shifted 50 units north it is probably an eastern or western neighbor, as opposed to 3000 units, which probably indicates a northern or southern neighbor. My polygon's units are feet and are all over 2000 feet along wide and over 3000 feet long. enter image description here

enter image description here

So in my process a neighboring polygon has to be a similar distance to the east or west to be a neighbor along the X-axis, and a similar distance away to the north or south to be a neighbor along the Y-axis.

My grid also doesn't have perfectly neighboring polygons, so I include a maximum distance apart of 50 feet to be counted as a neighbor.

The first iteration of the process would produce this: enter image description here

Then repeat, repeat, repeat until all polygons are updated.

#input feature class
inFc = r"C:\Path\to\your.gdb\featureclass"
#grid field
fld = "GRIDNO"

#max distance between neighbors as string (for selection)
dist = "50 FEET"

#minimum distance between neighbors for X and Y (in your feature class's unit of measuring length)
minXDist = 1700
minYDist = 2500

#-----

print "importing"
import arcpy
import os

#sql to create feature layer of populated grids
source = os.path.basename (inFc)
delimFld = arcpy.AddFieldDelimiters (source, fld)

#this process will repeat until a certain criteria is met
more = True
matchedOids = set ()

#continue pattern while variable 'more' set to True
while more:
    #set 'more' to False - will be reset to True if new polygons found
    more = False
    #layer of populated polygons
    print "creating first layer"
    sql = "{} IS NOT NULL".format (delimFld)
    result = arcpy.MakeFeatureLayer_management (inFc, "PopLyr", sql)

    #layer of unpopulated polygons
    print "creating second layer"
    sql = "{0} IS NULL OR {0} = '0000'".format (delimFld)
    result = arcpy.MakeFeatureLayer_management (inFc, "UnpopLyr", sql)

    print "iterating populated polygons"
    print arcpy.GetCount_management ("PopLyr").getOutput (0), "polygons to check"
    with arcpy.da.SearchCursor ("PopLyr", ["SHAPE@", fld, "OID@"]) as curs:
        for geom, nums, oid in curs:

            #skip polygons that have already been checked
            if oid in matchedOids: continue
            #add polygons to set of checked polygons
            matchedOids.add (oid)

            #get current polygon number
            yNum = int (nums[:2]) #y Num
            xNum = int (nums[2:]) #x Num

            #get extent object
            ext = geom.extent
            #get x min
            xMin = ext.XMin
            yMin = ext.YMin

            #select bordering polygons sharing line segment
            arcpy.SelectLayerByLocation_management ("UnpopLyr", "", geom, dist)
            #no selection? continue
            if not arcpy.Describe ("UnpopLyr").FIDSet:
                continue

            #found more polygons - set 'more' to True
            more = True

            #update cursor - iterate selected neighbors
            with arcpy.da.UpdateCursor ("UnpopLyr", ["SHAPE@", fld, "OID@"]) as uCurs:
                for uGeom, uNums, uOid in uCurs:
                    ext = uGeom.extent
                    uXMin = ext.XMin
                    uYMin = ext.YMin

                    #get difference between x mins and y mins
                    xDiff = xMin - uXMin
                    yDiff = yMin - uYMin

                    #skip cater-corner plots
                    if abs (xDiff) > minXDist and abs (yDiff) > minYDist:
                        continue

                    #logic time

                    #check polygon is west of polygon
                    if xDiff > 0 and abs (xDiff) > minXDist:
                        uXNum = xNum - 1
                        uYNum = yNum

                    #check polygon is east of polygon
                    elif xDiff < 0 and abs (xDiff) > minXDist:
                        uXNum = xNum + 1
                        uYNum = yNum

                    #check polygon is south of polygon
                    elif yDiff > 0 and abs (yDiff) > minYDist:
                        uXNum = xNum
                        uYNum = yNum - 1

                    #check polygon is north of polygon
                    elif yDiff < 0 and abs (yDiff) > minYDist:
                        uXNum = xNum
                        uYNum = yNum + 1
                    outYNum = str (uYNum).zfill (2)
                    outXNum = str (uXNum).zfill (2)
                    outVal = outYNum + outXNum
                    row = (uGeom, outVal, uOid)
                    uCurs.updateRow (row)

    #delete layers
    print "deleting layers"
    arcpy.Delete_management ("PopLyr")
    arcpy.Delete_management ("UnpopLyr")

print "done"

Result:

enter image description here

1

Here is a start of answer to question 1:

1) get the XY coordinates of your centroids (see Identify X Y coordinates of polygons and add them to the attribute table in QGIS? )

2) subtract the coordinates of the origin of your grid

3) compute the grid indices in the field calculator. The X and Y_spacings are the width and the heigth of your grid cells.

In arcpy

'{0:02d}{1:02d}'.format(int(!relative_Y_field_name!/Y_spacing), relative_X_field_name!/X_Spacing)

In QGIS (both rcpy and QGIS tags are present, I suggest you to split your question), it would be

 lpad(  to_string( to_int(relative_Y_field_name/Y_spacing)), 2,'0') +  lpad(  to_string( to_int(relative_X_field_name/X_spacing)), 2,'0') 

Warning : for security, do it in a new field or after selecting the 0000 in your feature class

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