J'ai une couche d'entités polygonale que je voudrais numeroter automatiquement mais en spiral (En escargot). De facon à avoir une continuité par voisinage. je suis un debutant en python, merci de m'aider. Merci
I see, you want to flip between max Y and max X and min Y and min X pairs. A Peano sort won't cut it. I suspect you will have to generate the spiral from the centre with a repeating angle and an increasing distance increment. The rate at which the distance from the centre changes (ie what function, if any, sin, cos etc) will determine the shape of the spiral.
I am sure it has a name, but the market for spiral generation outside of mazes and shellfish is somewhat limited 
Although a reversed squished Archimede's spiral looks promising https://en.wikipedia.org/wiki/Spiral
rigid solutions exist already in rosetta code https://rosettacode.org/wiki/Spiral_matrix#Python the trick is to oblate the spiral and smooth the steps
Interesting what you can get if you don't create a "rectangular spiral" as in the image the OP provided, but just follow the line of a somewhat oval spiral (since the input area is a rectangle and not a square):
However, this is not what the OP wants...
I do think that the spacing used between the lines of the spiral will influence the resulting order. This might be something to look at.
Just to add some explanation of what I did. I created two scripts, one to create the spiral and the other one to order the polygons.
The first script to create the spiral looks like this:
def main():
# formulas:
# x = phi cos phi
# y = phi sin phi
import arcpy
from math import pi, sin, cos
arcpy.env.overwriteOutput = True
# settings
fc_parcels = r'C:\GeoNet\Spiral\Parcels.shp'
fc = r'C:\GeoNet\Spiral\data.gdb\oval_spiral_v01'
sr = arcpy.Describe(fc_parcels).spatialReference
scale_extent = 1.25 # size it to 125%
pnts_cnt = 1000
pnts_div = 40.0
# determine output extent
ext = arcpy.Describe(fc_parcels).extent
ext = ScaleExtent(ext, scale_extent)
# run 1, create the normal spiral
lst_x, lst_y, lst_xy = [], [],[]
lst = range(pnts_cnt)
for i in reversed(lst):
phi = float(i) / pnts_div * pi
x = phi * cos(phi)
y = phi * sin(phi)
lst_x.append(x)
lst_y.append(y)
lst_xy.append((x, y))
# determine scaling of spiral
x_factor, y_factor, x_offset, y_offset = DefineScaling(lst_x, lst_y, ext)
# run 2, translate the spiral
pnts = []
for xy in lst_xy:
x = xy[0] * x_factor + x_offset
y = xy[1] * y_factor + y_offset
pnt = arcpy.Point(x, y)
pnts.append(pnt)
polyline = arcpy.Polyline(arcpy.Array(pnts), sr)
arcpy.CopyFeatures_management([polyline], fc)
def ScaleExtent(ext, scale_extent):
new_width = ext.width * scale_extent
new_height = ext.height * scale_extent
x_cc = (ext.XMin + ext.XMax) / 2.0
y_cc = (ext.YMin + ext.YMax) / 2.0
return arcpy.Extent(x_cc - new_width / 2.0, y_cc - new_height / 2.0,
x_cc + new_width / 2.0, y_cc + new_height / 2.0)
def DefineScaling(lst_x, lst_y, ext_out):
xmin, xmax = min(lst_x), max(lst_x)
ymin, ymax = min(lst_y), max(lst_y)
width_in = xmax - xmin
height_in = ymax - ymin
x_factor = ext_out.width / width_in
y_factor = ext_out.height / height_in
x_offset = ext_out.XMin - xmin * x_factor
y_offset = ext_out.YMin - ymin * y_factor
return x_factor, y_factor, x_offset, y_offset
if __name__ == '__main__':
main()
And it creates the spiral line.
Just to validate if it crosses all the polygons I manually did a select by location (but you can easily script this too).In case it doesn't select all the parcels, one should play with the values for pnts_cnt and pnts_div which define the number of spins.
In case it does select all the second step is to define the order of the parcels. This is done by first performing an intersect between the spiral line and the parcels. This too I did manually, but you can script this also. The resulting line featureclass will have all the intersecting segments. Next you will have to do a Multipart to Singlepart to ceate a feature for each part.
The order is then defined by the second script that validate the center of the segment on its position on the spiral. The lower this position (0-100%) the lower the assign value for the order.
def main():
import arcpy
fc = r'C:\GeoNet\Spiral\data.gdb\spiral_parcel_intersect_mp2sp_v01'
fld_id = 'FID_Parcels'
fc_spiral = r'C:\GeoNet\Spiral\data.gdb\oval_spiral_v01'
fc_parcels = r'C:\GeoNet\Spiral\Parcels.shp'
fld_out = 'Result'
# get data into dcts
dct_fc ={r[0]:[r[1], r[2]] for r in arcpy.da.SearchCursor(fc, ('OID@', 'SHAPE@', fld_id))}
spiral = arcpy.da.SearchCursor(fc_spiral, ('SHAPE@')).next()[0]
# create the dictionary with the position
dct_parcels = {}
for oid, lst in dct_fc.items():
polyline = lst[0]
fid = lst[1]
pntg = polyline.positionAlongLine(0.5, True)
pos = spiral.measureOnLine(pntg, True)
if fid in dct_parcels:
# there can be multiple intersects
min_pos = dct_parcels[fid]
if pos < min_pos:
dct_parcels[fid] = pos
else:
dct_parcels[fid] = pos
# assign the order based on position
dct_res = {}
cnt = 0
for fid, pos in sorted(dct_parcels.items(), key=lambda x: x[1]):
cnt += 1
dct_res[fid] = cnt
# add output field to the parcels
if len(arcpy.ListFields(fc_parcels, fld_out)) == 0:
arcpy.AddField_management(fc_parcels, fld_out, "SHORT")
# update the parcels
flds = ('OID@', fld_out)
with arcpy.da.UpdateCursor(fc_parcels, flds) as curs:
for row in curs:
fid = row[0]
if fid in dct_res:
row[1] = dct_res[fid]
curs.updateRow(row)
if __name__ == '__main__':
main()The result will look like this:
Since your spiral is more rectangular, the code would have to generate the "spiral" in a different way.
This could be a result using the "rectangular spiral":
Results depend highly on the spacing size...
def main():
import arcpy
arcpy.env.overwriteOutput = True
# settings
fc_parcels = r'C:\GeoNet\Spiral\Parcels.shp'
fc = r'C:\GeoNet\Spiral\data.gdb\rectangle_spiral_v10'
sr = arcpy.Describe(fc_parcels).spatialReference
spacing = 250.0
# determine output extent
ext = arcpy.Describe(fc_parcels).extent
x = ext.XMin
y = ext.YMax - spacing
pnt = arcpy.Point(x, y)
# loop
finished = False
pnts = [pnt]
xmin = x
xmax = ext.XMax
ymax = y
ymin = ext.YMin
done = []
while finished == False:
xmax -= spacing
while x < xmax:
xy = "{}#{}".format(x, y)
if not xy in done:
done.append(xy)
pnts.append(arcpy.Point(x, y))
x += spacing
ymin += spacing
while y > ymin:
xy = "{}#{}".format(x, y)
if not xy in done:
done.append(xy)
pnts.append(arcpy.Point(x, y))
y -= spacing
xmin += spacing
while x > xmin:
xy = "{}#{}".format(x, y)
if not xy in done:
done.append(xy)
pnts.append(arcpy.Point(x, y))
x -= spacing
ymax -= spacing
while y < ymax:
xy = "{}#{}".format(x, y)
if not xy in done:
done.append(xy)
pnts.append(arcpy.Point(x, y))
y += spacing
if (xmax - xmin) < spacing and (ymax - ymin) < spacing:
finished = True
# create the rectangle spiral
polyline = arcpy.Polyline(arcpy.Array(pnts), sr)
arcpy.CopyFeatures_management([polyline], fc)
if __name__ == '__main__':
main()I was playing with that idea in a slightly different fashion.
If you make a polyline of the extent, then select by location those parcels that touch it, you can order the outer 'spiral' beginning at the top left and proceeding clockwise until you reach the parcel just below your start parcel . You remove those parcels, form a new extent from the remaining parcels and proceed again. It was getting too slow, so I quit.
I think your ideas are good and I would keep the multipart intersection of the parcels with the spiral and the data can be carried over.
Which lines up well with radiating lines.
I could have sworn that there was a select by location option that accepted 'contains the endpoint of' but I couldn't find it. Would have made things a lot easier for me.
