Select to view content in your preferred language

points to polygons

6925
32
10-12-2012 12:29 PM
SiranErysian
Deactivated User
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.
0 Kudos
32 Replies
T__WayneWhitley
Honored Contributor
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.
0 Kudos
T__WayneWhitley
Honored Contributor
...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.
0 Kudos
T__WayneWhitley
Honored Contributor
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
0 Kudos
SiranErysian
Deactivated User
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
0 Kudos
T__WayneWhitley
Honored Contributor
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')
0 Kudos
SiranErysian
Deactivated User
Thank you. I"m going to give this a try and will let you know the outcome.
0 Kudos
SiranErysian
Deactivated User
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.
0 Kudos
T__WayneWhitley
Honored Contributor
...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
0 Kudos
SiranErysian
Deactivated User
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.
0 Kudos
T__WayneWhitley
Honored Contributor
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.
0 Kudos