nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

points to polygons

SiranErysian

I???m trying to delineate polygons along a river channel delineated by points representing different n values. I was wondering if anyone knew of a tool that would allow me to create polygons based on the same input value instead of what I am doing now which is connecting the dots. I tried Thiessen polygons but it gave me shapes that went way outside the boundaries of the points. Anyone know of such a function? See attachment for screen capture.

points2polygons.docx

Find more posts tagged with

Accepted answers

All comments

T__WayneWhitley

Right, so you must not be using ArcInfo (the Advanced license) - there are certain tools you must have the advanced license to run, Feature To Polygon is one of them.

EDIT: If you cannot borrow an ArcInfo license or you do not have time to study/modify the Python code to work around the licensing disadvantage, I suggest you download/install a free utility such as ET GeoWizards at ian-ko.com. The Build Polygon Wizard is completely free and works well.

SiranErysian

The tool ran fine I think because it created the routes, outer poly, channel and break polys. However, attached is a list of errors that came up after the tool completed. I will still have to manually add the n-values field to those polygons to see where they change. This information is in the point files. Aside from that, I think this should work for my purposes. I'm not sure if correcting the errors will give us those values I still need but I can enter them in a new field in the polygon feature classes.

linear_ref tool errors.docx

T__WayneWhitley

Okay, here's the deal -- I haven't heard back from you and I cannot use your incomplete attachment. But what I'll do is send you back a more automated 'package'....unzip the file attachment provided (making sure all files go to the same root location you choose; it is 'relative pathed' so should still find file refs) and run the script tool (open the toolbox and double-click the script tool 'linref' inside...

It is now geared to run off of your shapefile (the one you sent, also included on the provided attachment in case yours has changed).
When the tool opens, insert the shapefile - you'll see it alongside as in the illustration provided (also attached).

Provided you have ArcGIS Desktop 10 license (Advanced), then it should run in about 20 sec (as tested on my machine) to overwrite all you see contained in the file gdb. I left the test run feature classes in the gdb...but when you run it, it should overwrite everything. No 'addmessages' statements were added to the script, so it's very plain.

Good luck - oh, one more thing, for anyone interested, to 'adapt' this for a script tool, I did nothing more than change a few lines and create a shapefile input parameter in creating a toolbox script tool within Catalog...below are the critical line changes at the beginning of the script:

import arcpy, sys, os
scriptPath = sys.argv(0)
origData = arcpy.GetParameterAsText(0)
arcpy.env.workspace = scriptPath.replace(scriptPath.split(os.sep)[-1], 'stage.gdb')

That's all! Enjoy. Make sure the errors get posted next time, if there are any.

-Wayne

packagedToScriptTool.jpg

packagedToScriptTool.zip

T__WayneWhitley

Hello Siran...think your attachment may not have uploaded properly. I'll take a look probably tomorrow and it's likely a pathname or data type is set wrong. I don't know how well you understand the code, but certainly having output to examine will be essential for you. Maybe it'll help to make a zipped folder containing only relative pathnames so it will run without changing a thing.

Anyway, do you have any colleagues looking at this so far to give you a sense of confidence in the final output?....always good to do a reality check! Hope all is well.

SiranErysian

I tried running the script and got the following error messages attached. The polygons were not created. I would like help in tracing the errors to correct them and run the script again, if possible. What could be wrong according to the red errors? I think I Uploaded the files associated with this project to check if field names, etc are correct.

test_polys.py.zip

T__WayneWhitley

...sorry, Siran, I've been a little busy with other problems. I took a brief look at your attachments - so far so good except it doesn't look like you were able to achieve a finished channel polygon result? If you run the script as I have last posted with a similar line file as you have posted earlier, it should run to completion (you have ArcInfo, or the advanced ArcGIS Desktop, correct?). So, I have 2 questions/comments pertaining to the scripting:

- Did you have any error messages in the script run or do you have any specific related questions?
- You should attach your original script file -- do not copy to a Word document as this serves no purpose (and can create indention errors, failure in Python - so leave as is in the py file and attach that).

And I strongly suggest this (and you may have already followed up accordingly) -- check back with your associates familiar with the model output/techniques to be certain that your solution fits the problem. I am only familiar with the scripting techniques...not whether they are appropriate for your scenario. Consider that a final disclaimer. Just want you to know I enjoyed the challenge you posed, thank you, and if you are at a point you are satisfied with my 'contribution', please mark this post 'answered'.

Thanks again - if you have additional questions I can answer (scripting related), please be specific and I'll get back to you. Otherwise, looks good so far, so keep trying.

-Wayne

SiranErysian

I modified the attached code for my file names and have also attached the shapefile format of the feature classes I used. I would like you to check this code to see if it makes sense. I think I will only need you to look at the portion below the dashed line ----------------
as I already have the point events and route events as feature classes.
Let me know what you think--if I am on the right track.

Code-points to polygon-rev.docx

test.zip

SiranErysian

Thank you. I"m going to give this a try and will let you know the outcome.

T__WayneWhitley

I suppose I should add something here in an attempt to summarize better, at least what I've introduced, namely:
- the linear referencing to gain access to point measures on lines
- ordering point objects to create polygons
- a crude technique to 'break' the polygons at changing N values from the orig line file (I'm sure this technique will need work.)

A rudimentary script file has been included as an attachment - I would urge you to use this script at your own discretion. I'm not sure if the output will work in satisfactorily applying a solution for your original problem - however, this will at least introduce you to some geoprocessing techniques.

Good luck! I'll try to paste the code in a window below as well as include as attachment.
Note that if you test the code (which is in very simple format, no error trapping, etc.), the code works with your provided line file (so you'll need to change the shapefile pathname depending on where this source is and you'll need to create a file gdb workspace to which to write outputs (and change that pathname accordingly as well). There's nothing really special about this code - just a 'brute force' method going step-by-step as we've pretty much discussed in earlier posts.

origData = r'C:\Documents and Settings\whitley-wayne\Desktop\pointsToPolySerysian2\orig\Reach1B.shp'
import arcpy
arcpy.env.workspace = r'C:\Documents and Settings\whitley-wayne\Desktop\stage.gdb'
arcpy.env.outputCoordinateSystem = origData
arcpy.env.overwriteOutput = True

# Copy Features
# http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//001700000035000000
arcpy.CopyFeatures_management(origData, 'Reach1B_lines')

# Create Routes
# http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//003m00000005000000.htm
#CreateRoutes_lr (in_line_features, route_id_field, out_feature_class, measure_source, {from_measure_field}, {to_measure_field}, {coordinate_priority}, {measure_factor}, {measure_offset}, {ignore_gaps}, {build_index})
arcpy.CreateRoutes_lr('Reach1B_lines', 'ProfileM', 'Reach1B_routes', 'TWO_FIELDS', 'OID_', 'Shape_Leng')

fields = ['LeftBank_ft', 'RightBank_ft']
 
# Add Field
# http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//001700000047000000
# AddField_management (in_table, field_name, field_type, {field_precision}, {field_scale}, {field_length}, {field_alias}, {field_is_nullable}, {field_is_required}, {field_domain})

for field in fields:
     arcpy.AddField_management('Reach1B_lines', field, 'DOUBLE')

# Calculate Field
# http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//00170000004m000000
# CalculateField_management (in_table, field, expression, {expression_type}, {code_block})

for field in fields:
     exp = '!' + field.split('_')[0] + '!*!Shape_Leng!'
     arcpy.CalculateField_management('Reach1B_lines', field, exp, 'PYTHON')

# Copy Rows
# http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//0017000000n4000000
# CopyRows_management (in_rows, out_table, {config_keyword})

arcpy.CopyRows_management('Reach1B_lines', 'in_table')

# Make Route Event Layer
# http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Make_Route_Event_Layer/003m00000008000000/
# MakeRouteEventLayer_lr (in_routes, route_id_field, in_table, in_event_properties, out_layer, {offset_field}, {add_error_field}, {add_angle_field}, {angle_type}, {complement_angle}, {offset_direction}, {point_event_type})

in_event_properties = 'ProfileM POINT LeftBank_ft'
arcpy.MakeRouteEventLayer_lr('Reach1B_routes', 'ProfileM', 'in_table', in_event_properties, 'LeftBankPts')
arcpy.CopyFeatures_management('LeftBankPts', 'LeftBankPoints')

in_event_properties = 'ProfileM POINT RightBank_ft'
arcpy.MakeRouteEventLayer_lr('Reach1B_routes', 'ProfileM', 'in_table', in_event_properties, 'RightBankPts')
arcpy.CopyFeatures_management('RightBankPts', 'RightBankPoints')


# ...the channel poly building component:

# Array
# http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Array/000v0000005r000000/
polyArray = arcpy.Array()

inputPts = 'LeftBankPoints'

# SearchCursor (dataset, {where_clause}, {spatial_reference}, {fields}, {sort_fields})
# http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//000v00000039000000
rows = arcpy.SearchCursor(inputPts, '', '', '', 'ProfileM A')

for row in rows:
     pnt = row.Shape.getPart()
     polyArray.add(pnt)

inputPts = 'RightBankPoints'

rows = arcpy.SearchCursor(inputPts, '', '', '', 'ProfileM D')

for row in rows:
     pnt = row.Shape.getPart()
     polyArray.add(pnt)

# Polygon
# http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//000v000000n1000000

polygon = arcpy.Polygon(polyArray)
arcpy.CopyFeatures_management(polygon, 'channel')




# ...the first/last point extraction from lines, outer poly building component:

polyArray.removeAll()
inputLns = 'Reach1B_lines'

# Polyline
# http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#/Polyline/000v000000n2000000/

rows = arcpy.SearchCursor(inputLns, '', '', '', 'ProfileM A')

for row in rows:
     pnt = row.Shape.firstPoint
     polyArray.add(pnt)

rows = arcpy.SearchCursor(inputLns, '', '', '', 'ProfileM D')

for row in rows:
     pnt = row.Shape.lastPoint
     polyArray.add(pnt)

polygon = arcpy.Polygon(polyArray)
arcpy.CopyFeatures_management(polygon, 'outer_poly')


# ...the check for change in n-vals:

rows = arcpy.SearchCursor(inputLns, '', '', '', 'ProfileM A')
row = rows.next()

Nleft = row.Leftbank_N
Nchan = row.channel_N2
Nrght = row.rightbank_
row = rows.next()

# list to break polys
breakList = []

while row:
     if row.Leftbank_N != Nleft or row.channel_N2 != Nchan or row.rightbank_ != Nrght:
          breakList.append(row.ProfileM)
          Nleft = row.Leftbank_N
          Nchan = row.channel_N2
          Nrght = row.rightbank_
     row = rows.next()

del row, rows

# form feature layer where_clause
where_clause = '('

for each in breakList:
     where_clause = where_clause + str(each) + ', '

# AddFieldDelimiters
# http://help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//000v0000004n000000
# AddFieldDelimiters (datasource, field)

where_clause = arcpy.AddFieldDelimiters(inputLns, 'ProfileM') + ' IN ' + where_clause[0:-2] + ')'

# MakeFeatureLayer_management (in_features, out_layer, {where_clause}, {workspace}, {field_info})
arcpy.MakeFeatureLayer_management(inputLns, 'changeN', where_clause)
arcpy.CopyFeatures_management('changeN','breakLines')

# FeatureToPolygon_management (in_features, out_feature_class, {cluster_tolerance}, {attributes}, {label_features})
inFeatures = ['breakLines','channel','outer_poly']
arcpy.FeatureToPolygon_management(inFeatures,'break_polys')

To answer your question: yes, the endpoints of the lines should be where the outer boundary of the polygons would be delineated.

I now have n-values for each of the points based on the join I did with the values in the excel table and the same values for the cross sections. I selected all the points from the Reach_banks file that had the same n-values as the corresponding cross section. That gave me a field to calculate the values in that point file.

Yes, the join will be the input for the script for the entire river ( I have other reaches but we can practice on this one) The change in the N-value is where the break in the polygons should occur.

I have attached the cross sections and the banks point file. At this point I would need a script that uses the values in the cross section file to create polygons where the same n-values occur and breaks where they change for the left, channel and right. Will you be able to use the line file or the point file?

I will be out of the office until the 19th of November. We can touch base then

Reach1_Banks.zip

Reach1B.zip

T__WayneWhitley

Siran,
I'm still following this interesting thread.... I do not have any further product to post now but have a question about the outermost polygon boundaries (of the 2 outer polys) - is this delineated by the endpoints of the lines in your line shapefile or are there more values in the table I am not aware of? I should have the time tomorrow to post a complete script file solution. I'm updating you now with what I understand so it gets incorporated properly into the script.

If you've had time to read the last post, I suppose you now know that the channel points for the 'inner' polygon can be derived from your line shapefile via the linear referencing technique I outlined. It looks like in your attached data, you preprocessed a join from Excel to your line shapefile, ProfileM = Station, is that correct? ...and if so, is that join going to be the script input to use for the entire river as in the subset you provided earlier? If so, then I'm correct to assume the change in the sum of the 3 N values is where you want to 'break' the channel polygon? And do you also want to break the outer polygons at the same location, so that the final result is a break that runs the complete cross section line?

Wayne

T__WayneWhitley

...more time to look again at this - confirmed to my satisfaction that this part of the problem is really a 'linear referencing' exercise, see the attached (I only processed the right bank numbers). The red circles represent the points you sent in the shapefile and the points displayed within the circles represent the result of 'Make Route Event Layer'. Looks like I screwed up a bit on portions - this is simply where the measurement is from the 'wrong end', so that'll need a little tweaking. Once you get the points you can feed them into the script already provided to create polys from points.

So you'll probably need a brief description how I made the route event --

- 1st use 'Create Routes' on your line fc, Reach1B, and I just used your FID as a route identifier since each line is essentially a route.
(I don't think FID will be allowed as input to the tool, so create a 'dummy' numeric field, say FID1, and use the Field Calculator FID1 = FID)

- Use 'Make Route Event Layer' to create your points...you'll be using a table output from Reach1B since it contains the necessary RightBank and LeftBank measures. Also, remember we established that these measure values are actually proportion of the length of the respective lines, so you can establish the number of feet by multiplying the measure value by your val in Shape_Leng field (add a double field, use the Field Calculator).

EDIT - You have to use the 2 pt method in the 'Create Routes' process to get the correct result as shown in the attached linearRef_corrected.pdf.... See the attached jpg on how I filled out the dialog of 'Create Routes'. Basically you need a start measure of zero (0) --- I simply used a field with 0 already in the table, but if you need to create a 'dummy' numeric field, that'll work.

I needed to sharpen my skills on Linear Referencing anyway, so hope this helps you as well. You may need to 'Copy Features' of your route event output to make the input points for the script...try it out and see.

linearRef_corrected.pdf

linearRef.pdf

linRefCreateRtes_2point.jpg

T__WayneWhitley

Just thought I'd add this - I haven't had time to work on this further but I wanted to pass along that since RightBank and LeftBank are apparently 'measures', this has to do with 'linear referencing'...there are tools to generate points at those 'measure' locations, just I haven't used them in a while. I can also bypass that step using Python, but if you are familiar with the linear referencing tools, then you can also generate the points for the channel (I'm fairly certain that's how the point shapefile you have was generated) and use the script already made to build the polygons from the points.

I'll probably work more on this if I get time tomorrow, but keep me posted about this interesting problem.

SiranErysian

Thanks for touching base. This sounds great! Yes the N values pertain to flow. I would have never thought to look at the shape length. That would be great if the points could be computed with those values. Get back to me once you have tested your hunch.

T__WayneWhitley

OK, I've had a chance to look at this again... unfortunately it's toward the end of the day and I need to test a hunch since I'm not accustomed to this type of data. I'm just touching base with you before proceeding further. I don't think we need to compute anything , at least not for polygon geometry, with the fields pertaining to N values (is this not pertaining to flow?) -- I'm talking about these fields, according to what you attached:

channel_N2
leftbandN
rightbank_

However, I may have discovered a pattern with these fields (I have to test to see if my hunch is true though.):

LeftBank
RightBank
Shape_Leng

Looking at a single 'strand' cross-section line, say FID 171, with the values LeftBank 0.3254, RightBank 0.41533, Shape_Leng 1641.142808, it seems you can locate the points along the line with:

0.3254(1641.142808) = 534.028 (approx.)
0.41533(1641.142808) = 681.616 (approx.)

Then 681.616 - 534.028 = 147.588

...and if I carefully measure along that line FID 171 between the points in the shapefile you provided, using ArcMap's measuring tool, I get 147.58792.

So I'd say we have a winner algorithm for computing the points - I'll test further by 'chopping' your lines with the points for the segment lengths, make the computation as in the above, and make sure the results closely match the GIS computed length geom segments...pretty sure this isn't coincidental.

From there, I believe there are python methods for getting the 3 polygon sets you ask. Well, it's looking up, what do you think?

SiranErysian

In response to your message: Paragraph 1: I don't have point geometry for all the channel. It was derived from a model only for that section for some reason. So I have to use the values in the cross section lines. Forget about the points. Paragraph 3: The lines represent the cross section of not just the channel but the left and right bank values for flow.(N) The polygons will represent the values from line to line connecting the same values as polygons--for the channel, left bank and right bank. So there will be 3 sets of polygons. You are correct that the points represented the channel but I don't have that data for the rest of the river. Paragraph 4: How do we derive the inner channel or any of the polygons if we don't have xy locations? In the previous correspondence, I was assuming that we could assign coordinated based on some value in the cross section line file? If it is not possible to form polygons for all three areas, left, right and channel, maybe we can just do the channel.

T__WayneWhitley

What I'm saying is, if you have it, you can directly use point geometry (as I last demonstrated) or you can use XY coordinate values if you have that in a table, Excel, txt, or whatever. Or if you can derive those points somehow from the lines - I didn't see how, maybe I'm just being dense, but you know your data better than I.

Only real difference in using coord vals is you'd need to fetch the coord pairs from the table to assign to point objects before loading them in the array (the last sample simply fetches the point objects themselves from the point shapefile).

I'm slightly confused about 1 thing: what exactly you want polygons to represent.
The channel cross-sections bank-to-bank is what the lines represent, correct? But the points shown before were, as I understood it visually, where the 'normal' or main channel is? So you do not want a polygon bank-to-bank but the actual channel? Correct me if I'm wrong in these assumptions, and tell me how these 'channel points' are determined and did I somehow miss that from the table?

At this point (if the above assumptions are correct), I'm thinking if the 'inner' channel points are related to the measure of the line cross-sections, then the point locations you want to form the channel polygon out of can be derived. Just curious, why do you only have points for part of the channel?

Anyway, I can look at this again probably tomorrow and figure it out -- fill me in when you can and we'll solve this quicker.

SiranErysian

Thanks for your last message and the code. I hope I haven't been sending you around in circles. Here is where I am: I have values in Excel that correspond with the left, channel and right banks of the river. I joined this table to the line file based on a common field. (I am not going to use the point file since the values exist only for a short section of the river and I need the whole). So are you saying that the code will assign coordinates based on those values to a line file with no need for xy values? See attachment of values along the cross sections. Top label is left, middle is channel and bottom is right bank. As I mentioned earlier, where the values change is where the break point would be. Can I substitute the shp file in the code with another and can I substitute lines for points?
Is there anything else on the code that would have to be revised?

cross section N values.docx

T__WayneWhitley

IDLE is simply the interactive editor/debugger app I use to test commands as I am typing each one, sort of like the command window you can now open within ArcMap to execute Python. Basically, if you take out the command prompts " >>> " and a few extraneous print commands, etc., get the indentation right, then you have a working script.

Not sure I yet understand what form all your data is in...are you saying you have XY coordinate values in Excel? Is there a reason why you have data in multiple forms, or is this because you just haven't assimilated it yet, or is all part of your data conversion process?

Yes, if you have coord vals in Excel, a search cursor can access that directly and assign each coord pair to a point object (in the code provided, the search cursor accessed the point shapefile and captured each point object directly - no need for the XY values).

Note: Not sure what you mean by creating polys "where values change" but if you happen to have ID values in Excel as in the input file used in the sample of the Writing Geometries help file, you can use this as your means to loop on a specific set or range of values, breaking the polygons where you wish.
Also, about the spatial join, I was using that because I thought you were missing a reliable common field on which to make an attribute join. I didn't need this to create the points although it appears this extra data may be needed to define where to 'break'.

SiranErysian

Yes I get the code. Will need to study it further. Not sure what IDLE is. I would like to use an excel spreadsheet I have with all the values I need: left bank, channel, right bank, as you mentioned above. Will that work? Need to create polygons where values change at each of these locations. I only had points for a short section. None were created for the other sections. I am assuming that no spatial join is needed if I have all the data in an attribute table? Just join the excel table to the line file? Thinking out loud...will get back to you once I study the code a bit more and see if I can make it work on some of the data.

T__WayneWhitley

I meant to get back to you with this. Basically, it's a demo of IDLE commands that can be refined as you wish into a finished script. All I did was this:

- Took a subset of your point shapefile as defined in a query on the search cursor.
- Loaded the point geometry objects into a python list.
- According to the order you wanted based on your input, used the 'device' already described above for 'odd' and 'even' sides.
- Resorted the point objects to form the perimeter of a polygon.
- Created a new row via insertcursor in an existing polygon shapefile into which the point objects are to define new poly geom.
- Loaded the result into an object array (thought I could load the py object list directly into the insert row, but that didn't work).
- Loaded the object array into the Shape field for the new row. (See the pdf view of the result attached.)

Here's the rough code (again, this was done by trial and error within IDLE), see below. If you want to create 1 huge poly out of the point shapefile, you can - simply remove the query. Otherwise, if you want to create individual polygons out of progressive sets of points from the input shapefile, you can set up a looping cursor either based on a common ID code in the point table or redefining your query. Hope that helps.

>>> import arcpy
>>> arcpy.env.workspace = r'C:\Documents and Settings\whitley-wayne\Desktop\pointsToPolySerysian\Reach1_Banks'
>>> polyArray = arcpy.Array()
>>> where_clause = '"FID" > 1131 AND "FID" < 1148'
>>> inputPts = "Reach1_Banks.shp"

>>> pointsLeft = []
>>> pointsRight = []
>>> points = []

>>> rows = arcpy.SearchCursor(inputPts, where_clause)
>>> for row in rows:
 pnt = row.Shape.getPart()
 points.append(pnt)

>>> for i in range(0,len(points)-1,2):
 pointsLeft.append(points)
 pointsRight.append(points[i+1])

>>> len(pointsRight)
8

>>> pointsRight2 = []

>>> for i in range(len(pointsRight)-1,-1,-1):
 pointsRight2.append(pointsRight)

>>> len(pointsRight2)
8

>>> for eachPoint in pointsRight2:
 pointsLeft.append(eachPoint)
 
>>> len(pointsLeft)
16

>>> cur = arcpy.InsertCursor("testPolys.shp")
>>> feat = cur.newRow()

>>> for eachPoint in pointsLeft:
 polyArray.add(eachPoint)

 
>>> feat.Shape = polyArray
>>> cur.insertRow(feat)
>>> del cur
>>>

demoSpatial_.pdf

T__WayneWhitley

Okay, I've had a chance to get back to this - I haven't done much, just taken a look and attached a couple of files. Looks more like a 2-part problem, that maybe you'd like to 'see' data in a single table? Otherwise the problem statement is as you stated from the beginning to create a polygon based on ordering of FIDs, odd on one side, even on the other. Is there any change to that?

Pertaining to the data table, it seems to me you could do a spatial join, to 'combine' or join the line attribute data to the points. Make sure this data is accurate, see the attached. I used the one-to-many option of Spatial Join making the points the target so that you'd be able to 'look' at the data in a single table. Still it is redundant though - you could condense the table further to see a single record for each line (reversing the join, more or less). Anyway, I don't see how this really affects your orig problem creating polygons - have I missed something?

If all is as I understand it, I'll run the solution and attach it tomorrow.

SpatialJoinResults.pdf

spatialJoin.zip

T__WayneWhitley

OK, I looked briefly at your attachments and will look further into it when I have more time...probably later today.

I would suggest 1st that if you know how to better combine your data, I would write everything to a temp workspace - that way, your orig data is preserved and you hopefully will have intermediate output that's simply easier to work with (or at least appears less complicated).

Sounds like you may benefit from either joining, sorting, or otherwise housing the data in a single table? That way, when you go building the geom, you then don't have to go to multiple lookup tables, etc., for attributes if it isn't necessary. I haven't looked into it yet at that detail, but sometimes it's more expedient to work with what scripts you've already experimented with and go ahead if at all possible to 'clean up' or adapt the data accordingly... I'll get back to this shortly. Meanwhile, if you have a breakthrough, certainly post it!

SiranErysian

I have attached 2 files, one showing points Reach1B_Banks) but only for part of the stretch of river I am trying to calculate, and the other a line file (Reach1B). Values were not included from the spreadsheet but are included in the cross section line file. I hope the script can handle values in this file. What I need to do is combine the N_values for the left bank, channel and right bank which are in the attribute table as Left_Bank_N, Channel_N and Right_bank_N. These values will determine the changes in flow. The screen shot shows how I connected the dots but Reach 1B point does not have values so I want to use the red lines which have values.

Does this make sense? Will it work?

I will get back to you either today or Tuesday witha sample data set.

T__WayneWhitley

Certainly. Can you attach a sample of your data? I ask because at this point it would probably be easiest to write the code specifically adapted to your input, whether that be a point feature class with these n-values in the attribute table or other format, such as a text input file containing these values in addition to XY point coords. It would be easier to explain it with comments in the code which I can post here. In short, I can give you the below explanation (along with that above for sorting the input).

If you already have a point feature class, you can use a sorted search cursor to access each table row and transfer point object geometry directly to an object array. Or, if you have a text file (or just a table, no geometry yet), you can read each line in the file (or row in the table), create point geometry from the read XY coordinates, and load the resultant points into an object array.

So how are the sets defined during iteration?...i.e., when do you stop loading the object array to load the point object set into a feature polygon geometry? By detecting change in the aforementioned n-values read with each row (or, if txt file input, with each line) - when the n-value changes, that's the end of the set. The array load is then interrupted to load the array into the polygon feature's geometry which in turn is inserted into a destination polygon feature class (your output feature table).

The iteration continues on in similar fashion with each successive set... there won't be any n-value change detected on the last set, so the final set will of course load the array but will not be loaded into polygon feature geometry within the loop - this final polygon must be completed outside the loop.

SiranErysian

Thanks. I have been off on another project but will get back to this in the next day or so.

You mentioned closing the polygons: I'd like to do this where the n-values change (manning's coefficient value) which is in the attribute table as n-value. Would this be possible to insert somewhere in the code or I can just split the polygon created at those locations? See attachment. Upper number in label is the FID and lower number is the n-value.

N-value and FID.docx

T__WayneWhitley

Well, then, once you've digested that, getting your point objects into the list (I guess you'll be doing this with a searchcursor), the rest isn't too difficult. This may not be the way you go exactly, but could be an example, if not too premature to load you up with more to read---

# these are point objects (no longer text)
points = [point1, point2, point3, point4, point5, point6]

array = arcpy.Array()

for point in points:
     array.add(point)
     
# close the poly
array.add(points[0])

# load the array into geometry object
polygon = arcpy.Polygon(array)


# Here, you'll probably want an insertcursor to loop and load geometry for multiple features...
# But for demo on a single feature, you can use copy features:
arcpy.CopyFeatures_management(polygon, r'C:\JustTesting.shp')

SiranErysian

Wow, Thanks. I'll have to read this over several times to absorb it. Not a real expert on Python but its practical stuff and makes sense when I read it. Will get back to you.

T__WayneWhitley

So you really need to connect in an order as in, for example:
['point1', 'point3', 'point5', 'point6', 'point4', 'point2']

Basically, your solution lies in that with arcpy geom methods you can load points into an array to in turn load into a line or polygon geometry object. The 'puzzle', if you will, is what sequence you feed the points in... it's late, so I don't have the finished code sample at the moment, but consider this (and make sure it is based on accurate assumptions):

You said you had points in a sequence, odd on one side, even on the other. I've seen your picture - just not sure yet how the 'polygon' closes, but this is not as critical as it is to at least get an algorithm going for 'connecting the dots'...

So what you need is a device for loading odd points, then even points...
-- (This may not be optimal code, but it'll work.)---

# For example say instead of the string list shown, the following is a python list of point objects:
>>> points = ['point1', 'point2', 'point3', 'point4', 'point5', 'point6']

# Define 2 new lists to load for left and right sides:
>>> pointsLeft = []
>>> pointsRight = []

Now you need a looping mechanism to load the empty lists:
>>> for i in range(0,len(points)-1,2):
 pointsLeft.append(points)
 pointsRight.append(points[i+1])

# Now the lists are 'segregated' odd and even:
>>> print pointsLeft
['point1', 'point3', 'point5']

>>> print pointsRight
['point2', 'point4', 'point6']

# Reverse the order, say on the right, to set up 'connecting' the 2 lines:
>>> pointsRight.sort(reverse=True)

>>> print pointsRight
['point6', 'point4', 'point2']

# Now the list can be 'mashed' together:
>>> for eachPoint in pointsRight:
 pointsLeft.append(eachPoint)

# Still called 'pointsLeft', but it's really a sequence around both sides... 
>>> print pointsLeft
['point1', 'point3', 'point5', 'point6', 'point4', 'point2']
>>>

Join the Conversation

Signed in members can post, follow updates, and more. New here? Register a free account.

Quick Links

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.