J'ai une couche d'entités polygonale que je voudrais numéroter automatiquement mais en spirale (En escargot). De facon à avoir une continuité par voisinage. je suis un débutant en python, merci de m'aider. Merci

J'ai une couche d'entités polygonale que je voudrais numéroter automatiquement mais en spirale (En escargot). De facon à avoir une continuité par voisinage. je suis un débutant en python, merci de m'aider. Merci

See Dan Patterson's comment in this thread: numbering polygon in clockwise. Although the message doesn't have a code example, it gives a suggestion on how this problem might be approached.

- 2 people found this helpful
`def angle_2pnts(p0, p1):`

"""Two point angle

: angle = atan2(vector2.y, vector2.x) - atan2(vector1.y, vector1.x);

: Accepted answer from the poly_angles link

"""

ba = p1 - p0

ang_ab = np.arctan2(*ba[::-1])

return np.rad2deg(ang_ab % (2 * np.pi))

p1 = np.array([1, 1])

p0 = np.array([0, 0])

angle_2pnts(p0, p1)

45.0Sadly for some, I use numpy to perform all my calculations, but this may help

p0... consider that your center point

p1.... pN are all your other points.

If you had a list of points, I could modify my example or you could translate this into straight python using elementary geometry

In short, the angles are relative to the X-axis ranging from 180 to -180... sort in descending order to give you a clockwise direction, sort in ascending order to get counterclockwise.

If you have data in a particular form, I can provide a field calculator expression, and I have a toolbox somewhere that does radial, lexicographic and other sorting types if needed (I have probably done a blog post at some time as well)

Ok... I found some field calculations.

If you have a featureclass/shapefile, this field calculator will do the calculations for you.

You must have projected coordinates!

If you copy that file into a text editor and save it with as .... azimuth_to.cal .... you can load it up into the field calculator. It is the same principle, but it does the azimuth to... but you can then sort it using

**Sort**.... since you are only sorting on a regular field.Normally the Sort tool (in Map or PRO) is only available at the advanced license if you want to sort on the Shape field or to use multiple fields.

`# -*- coding: UTF-8 -*-`

"""-----------------------------------------

Input shape field: returns angle between 0 and <360 based upon the first and last point

azimuth_to(!Shape!,from_x, from_y, from_north=True)

ie azimuth_to(!Shape!, 300050, 5000050, True)

"""

import math

def azimuth_to(shape, from_x, from_y, from_north):

x = shape.centroid.X

y = shape.centroid.Y

radian = math.atan2((y - from_y), (x - from_x))

angle = math.degrees(radian)

if from_north:

angle = (450 - angle) % 360

return angle

__esri_field_calculator_splitter__

azimuth_to(!Shape!, 300050, 5000050, True)NOW see the last line... the 300050 and 5000050 are the centre X and Y values... you will need to use the python parser in the field calculator and determine and enter your central points. If you want them from North, then entre True, otherwise you get the angles from the X-axis.

If set to False, then you have to sort in descending order.

if set to True, then you will be sorting in ascending order.

Hope this is clear as mud... draw a picture if it helps.

Hi Ned Charles ,

If you want to number them as a "spiral", I assume that distance from the first (center) also influences.the resulting numbering as well as the start angle:

Can you provide a screenshot or sample of your data to see what it looks like?

- 20/5000image in attachment
The parcels are in red and the black line is the neighborhood numbering direction.

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

Dear Dan,

Based on our last discussions(Spiral matrix !) , I see that you understand the problem very well. Your beginnings of solution are very good but I would like to be able to apply them to a layers of polygons. Thank you for explaining to me how to number this layer of polygons.Thank you so much, we are close to the solution

- 1 person found this helpful
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.

- 2 people found this helpful
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.

- 1 person found this helpful
Looks good! If distance and direction to the parcel centroids are needed, I can provide those pieces to supplement this.

The minimum spiral tightness to ensure complete selection of all the parcels could be an iterative process, but I suspect 'corner-cases' can be found for any implementation.

- 2 people found this helpful
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()- 2 people found this helpful
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

but I couldn't find it. Would have made things a lot easier for me.*'contains the endpoint of'*

Dear

**Xander Bakker****The "rectangular spiral" seems to me the most suitable.**Since I am a beginner, I can not run the script by copying it to the Python window on ArcGIS. I would like, if possible, to have a procedure to be able to execute it correctly. For this, I send you a small group of parcels (168) in a geodatabase. Please help me understand the parts of the script that I have to modify later for other lots of parcels.

Many thanks.PS: Some print screen could help too.

- 2 people found this helpful
Ned a good example of what is going on conceptually is given here https://stackoverflow.com/questions/36834505/creating-a-spiral-array-in-python

You can think of your problem in terms of blocks.

Pictorally it makes more sense of the blocks have an odd number of rows and columns.

The only difference in what is shown here is that the

**dx**(change in x) can be different in the**dy**(change in y) perstepSo here is a sample of a study area with each 'cell' numbered from 1 to 49

`# ---- 7 rows by 7 cols, number from 1 to 49`

a = np.arange(1, 7*7+1).reshape(7, 7)

array([[ 1, 2, 3, 4, 5, 6, 7],

[ 8, 9, 10, 11, 12, 13, 14],

[15, 16, 17, 18, 19, 20, 21],

[22, 23, 24, 25, 26, 27, 28],

[29, 30, 31, 32, 33, 34, 35],

[36, 37, 38, 39, 40, 41, 42],

[43, 44, 45, 46, 47, 48, 49]])Conceptually, if you wanted to move in the above array in a clockwise order starting at the top left, your steps would be in the form of

`spiral_cw(a)`

array([ 1, 2, 3, 4, 5, 6, 7, 14, 21, 28, 35, 42, 49, 48, 47, 46, 45, 44,

43, 36, 29, 22, 15, 8, 9, 10, 11, 12, 13, 20, 27, 34, 41, 40, 39, 38,

37, 30, 23, 16, 17, 18, 19, 26, 33, 32, 31, 24, 25])where you begin at position 1, move to 7, drop down to 14 then 21, 28, 35, 42, 49 then left to 43, then up through 36, 29, 22, 15, 8 etc etc

This information (see the code implementation) enables you to construct a clockwise spiral.

And you want to translate it into this form of a clockwise set of block

`print(to_spiral(a)) # using the code from the reference`

[[43 44 45 46 47 48 49]

[42 21 22 23 24 25 26]

[41 20 7 8 9 10 27]

[40 19 6 1 2 11 28]

[39 18 5 4 3 12 29]

[38 17 16 15 14 13 30]

[37 36 35 34 33 32 31]]So you can either build your study area divisions by knowing the coordinates of the top left 'cell' and proceeding in a logical order changing your position at the end of the line...

OR ... by beginning in the middle and spiraling outward in a clockwise direction (

So how the block coordinates are build can be done many ways. All you need is a start coordinated and a dx, dy to make your step moves.

- 2 people found this helpful
So I just downloaded your data and had a go at it and this is the result:

To do that I changed the following settings in the first script:

`# settings`

fc_parcels = r'C:\GeoNet\Spiral\Python.gdb\Data\PARCELS'

fc = r'C:\GeoNet\Spiral\Python.gdb\rectangle_spiral_v01'

sr = arcpy.Describe(fc_parcels).spatialReference

spacing = 10.0- change the path to the input fc_parcels and define the output path to the fc with the "rectangle spirals"
- I defined a small spacing, but you could use a larger value. Just make sure that every parcels is crossed.(use a select by location to validate this)

After this there are two steps that I did not include in the code.

First intersect the parcels with the rectangle spiral

The I did a multi part to single part to make sure that every part is a single feature:

For the second script I defined the following settings:

`fc = r'C:\GeoNet\Spiral\Python.gdb\intersect_spiral_parcels_mp2sp_v01'`

fld_id = 'FID_Parcels'

fc_spiral = r'C:\GeoNet\Spiral\Python.gdb\rectangle_spiral_v01'

fc_parcels = r'C:\GeoNet\Spiral\Python.gdb\Data\PARCELS'

fld_out = 'Result'- change the paths to the input single part feature class, point to the correct spiral featureclass (the one before the intersect), point to the parcels featureclass and define the name for the output field.
- Make sure that Parcels featureclass is not loaded (locked) by ArcMap

See the result attached below.

Thank you very much dear Xander.

Life is worth living, it is because there are people like you!

Many thanks!

Another question: Is there a way to create a spiral from an extent of a polygon?

Thank you for informing me if possible.- 2 people found this helpful
The answer to your question would be yes, since you can simply read out the extent of a polygon instead of the extend of the featureclass. However, after looking at the image I get the idea that you are looking for something else ... which could require rotation of each "side" and could involve a parallelogram rather than a rectangle:

This will be more complex to create...

- 2 people found this helpful
Just tried something using an internal multiring buffer so irregular shapes can be used in combination of the start angle and taking into account in which ring the maximum area of the parcel falls:

However this still creates some odd result (though explainable due to the angle). See labels 5 to 10:

Might be interesting to know what you exactly trying to achieve and why this is necessary.

Translation: I have a layer of polygonal entities that I would like to automatically number but spiral. In order to have continuity by neighborhood. I am a beginner in python, thank you for helping me. Thank you

Question: Do you want the spiral in a clockwise or counterclockwise direction? ( Voulez-vous la spirale dans le sens des aiguilles d'une montre ou dans le sens contraire?)

Tagging Python