Cut polygons with grid

[UlrichStroetz]

I have a several polygons with certain statistical values. I want to cut the polygon with a fishnet (1000m*1000m). Eventually I want that the tiles of the grid has the percentual values of the polygon and I want to convert that to a raster.

Any ideas how to do this? Picture is attached.

[ATTACH=CONFIG]18319[/ATTACH]

[T__WayneWhitley · Accepted solution]

Sounds like a job for 'Make Feature Layer' and 'Union' ?  ...then one of the 'To Raster' conversion tools?  (or can you use a conversion tool directly?)

The idea with preprocessing with 'Make Feature Layer' is to apply a query if you need it and a percentage ratio based on the original geometry area before feeding into 'Union' with the grid.  This 'ratio' is applied in the 'Field Info' section of the Make Feature Layer tool - you specify the field which you want to 'Use Ratio Policy' - simple example is use a double field originally calculated to 1.0 (or 100.0) for per cent) but you probably want to use this function on other values?  Feeding the MFL output into 'Union' (it will be in the 'drop-down' on inputs), the output of 'Union' the policy will be applied (any portions cut off within grid bounds, say half in one grid half in the adjacent grid, the values with ratio policy applied should be also split in half).  Just an example - I'm not convinced this is appropriate for what you want though.

[MikeAlonzo]

Hello,

This thread also mostly addresses a question I had about the split policy with MakeFeatureLayer. I know that it is the right tool for me (as I distribute population data across intersected polygons) but I'm having trouble with the python syntax and the help item just refers the user to: " See the usages for more information." (and I don't understand where this is telling me to go.)

So, could I get someone to copy in the syntax that you used to specify the "field_info" field and then specify that you wanted ratio splitting? Also, if I'm missing some "usages" section of the help text, please let me know.  Thanks much!

Mike


........
Mike Alonzo
UCSB Geography

[T__WayneWhitley]

This should help you out - Make Feature Layer allows access to a fieldInfo object which in turn allows making use of the ratio policy:

�?�A split policy can be set by using the Field Info control Use Ratio Policy option.

Make Feature Layer (Data Management)
Desktop » Geoprocessing » Tool reference » Data Management toolbox
http://resources.arcgis.com/en/help/main/10.2/index.html#//00170000006p000000


There's a nice sample script from the webhelp page for FieldInfo (bottom of page) that reads or 'gets' the 4 fieldInfo properties - the input field name, the new (output) field name, the split rule (ratio policy), and the visibility of fields.

FieldInfo (arcpy)
Desktop » Geoprocessing » ArcPy » ArcPy classes
http://resources.arcgis.com/en/help/main/10.2/index.html#/FieldInfo/018z0000004v000000/

"RATIO �??The attributes of resulting features are a ratio of the original feature's value. The ratio is based on the division of the original geometry. If the geometry is divided equally, each new feature's attribute gets one-half of the value of the original object's attribute."


So what you could do after creating your fieldInfo obj is to find the field you mean to 'set' the ratio policy for, with:

findFieldByName (field_name)  #Finds the field index by field name

...then use the following setSplitRule method with the proper index to set the rule (which in your case would be to enable it, the text 'RATIO')
setSplitRule (index, rule)


That's probably all the info you need, but if you cannot follow that lead I can post more specific code...

Hope that helps,
Wayne


EDIT:
This is untested code below, but should do the trick - once the 'layer' object ("temp_layer") from MFL is 'stamped' with the modified fieldInfo control, you should be able to use it in another overlay process to employ the RATIO policy.

import arcpy

feature_class = "c:/Data/wells.shp"
layer = "temp_layer"
arcpy.MakeFeatureLayer_management(feature_class, layer)

# Create a describe object
desc = arcpy.Describe(layer)

# Create a fieldinfo object
field_info = desc.fieldInfo

# Get index of desired field for which to set ratio policy
# Let's say your desired field is called 'area'
fldIndex = field_info.findFieldByName("area")

# Set the policy
field_info.setSplitRule(fldIndex, "RATIO")

[MikeAlonzo]

Thanks, Wayne. This all worked out well based on your instructions.

Mike

[MikeAlonzo]

Hi Wayne --

Sorry to keep this going but I've had trouble making this work in a stand-alone setting. I used the code above, which seemed to work fine. However, when I proceeded the intersect the ratio split policy wasn't honored. I'm wondering if I'm missing a step to connect the policy more explicitly back to the layer that is input to the intersect.  Code below.

Thanks for any thoughts,
Mike

inFeaturesAlb = "blocks_01_alb.shp"

layer = "temp_layer"
arcpy.MakeFeatureLayer_management(inFeaturesAlb, layer)

# Create a describe object
desc = arcpy.Describe(layer)

# Create a fieldinfo object
field_info = desc.fieldInfo

# Get index of desired field for which to set ratio policy
# Let's say your desired field is called 'area'
fldIndex = field_info.findFieldByName("POP10")

# Set the policy
field_info.setSplitRule(fldIndex, "RATIO")  

######################################################################    
#Do the intersection
#####################################################################     
intOut = "State_01_int_3.shp"    
inFeatures_int = [USA_layer, layer]

arcpy.Intersect_analysis(inFeatures_int, intOut)

[T__WayneWhitley]

Ratio policy should be honored by Intersect; however, I see you're feeding multiple input layers (as a list) and I don't see any ratio policy being defined for one of them in your code.  Without testing, my guess is that as it's being fed from the list, either the policy you intended to set isn't being honored because it's missing for the other layer or it is simply lost in translation.

If you cannot get it to work feeding in the layers in a single Intersect processing, to at least see if you can get this to work, why not feed in separate executions of Intersect for each of your inputs?  That'll be fine particularly if you're only interested in the ratio on one of the layers.  ...but I'd be interested to know if it is as simple as making certain ratios are set for both layers.  Also, this is probably a given and you've already got it covered, but this policy only applies to numerical fields.

I'll test this later if you're still having problems.

Hope you get it running,
Wayne

Also, see this blog below, pretty interesting and relevant...search the page for Ratio Policy and you'll see there's a topic where the policy is applied in a model - you could also very quickly set up a model for testing purposes, then proceed to committing to script.

More adventures in overlay: splitting polygons with cartographic spaghetti
by Dale Honeycutt on November 9, 2012
[Apportioning the NTrees attribute]
http://blogs.esri.com/esri/arcgis/2012/11/09/splitting_polygons/


EDIT:
One more possible issue about your method - I noticed you appear to be trying to 'apportion' what looks to be a population field -- you have to be very careful about this because population is related to density and apportioning by area will incorrectly skew a density figure.

[RichardFairhurst]

Hi Wayne -- 

Sorry to keep this going but I've had trouble making this work in a stand-alone setting. I used the code above, which seemed to work fine. However, when I proceeded the intersect the ratio split policy wasn't honored. I'm wondering if I'm missing a step to connect the policy more explicitly back to the layer that is input to the intersect. Code below. 

Thanks for any thoughts, 
Mike 

inFeaturesAlb = "blocks_01_alb.shp"

layer = "temp_layer"
arcpy.MakeFeatureLayer_management(inFeaturesAlb, layer)

# Create a describe object
desc = arcpy.Describe(layer)

# Create a fieldinfo object
field_info = desc.fieldInfo

# Get index of desired field for which to set ratio policy
# Let's say your desired field is called 'area'
fldIndex = field_info.findFieldByName("POP10")

# Set the policy
field_info.setSplitRule(fldIndex, "RATIO")  

######################################################################    
#Do the intersection
#####################################################################     
intOut = "State_01_int_3.shp"    
inFeatures_int = [USA_layer, layer]

arcpy.Intersect_analysis(inFeatures_int, intOut)

The help has a warning that "Geoprocessing tools do not honor geodatabase feature class or table field split policies." The help for tools like Intersect and Union also say "However, if the input is a layer or layers created by the Make Feature Layer tool and a field's Use Ratio Policy is checked, then a ratio of the input attribute value is calculated for the output attribute value." (Frankly, whoever wrote the help on this topic could have done a much better job, since the help layout never clearly links the warning and solution together).

I would set the split policy rule with the MakeFeatureLayer_Managment tool directly. I believe the help means to indicate that geoprocessing requires the use of this form of the split policy for other tools to honor it. If you build the Make Feature Layer Management tool in the ModleBuilder environment and check the split policy flag for the field(s) you want and export the model to a python script you will get the syntax you need. Here is an example of what ModelBuilder outputs:

arcpy.MakeFeatureLayer_management(inFeaturesAlb, layer, "", "", "OBJECTID OBJECTID VISIBLE NONE;APN APN VISIBLE NONE;POP10 POP10 VISIBLE RATIO")

The field info portion of the tool is a single string that contains a semi-colon separated list of the layer fields, where each field is defined with a space separated list of four items:

The Original Feature Class/Table field name
The Layer display field name
The Visibility flag (either VISIBLE or NONE)
The Split Ratio Policy flag (either RATIO or NONE)

Additionally, I don't think shapefiles support a split ratio policy at all. Only geodatabase feature classes and tables support a split policy. Therefore your script will never work as long as your input is a shapefile. Your script might work as you have written it if you changed the feature class input from a shapefile to a geodatabase feature class.

Finally, only numeric fields honor the split policy and only Double and Float fields will honor the policy with an acceptable degree of accuracy.

[MikeAlonzo]

Thanks for the additional insight. The testing with model builder was a good call. The code generated from that (that worked fine) below.

It is still somewhat unclear to me why my original code doesn't honor the split though. I don't believe the shapefile is the issue since I am creating a feature layer using MakeFeatureLayer_management. And, when I check the the field_info properties using getsplitrule after creation I see, for my numeric field in question, that "RATIO" is listed.

I do appreciate the explicit connection between the created layer and the split policy that is evidenced in the model builder code version though.

As to Wayne's point about apportioning population: I'm with you about the population density issues. However in this case this is census block data and the spatial resolution is fine enough where we're not worrying about minor spurious re-apportionment. It is actually useful to have sliver polygons that are output from the intersection have populations of 0 or 1 so we can catch the error easier.

import arcpy

# Local variables:
blocks_01_alb_shp = "...\\blocks_01_alb.shp"
USA_int_131111 = "USA_int_131111"
State_01_int_7_shp = "...\\State_01_int_7.shp"
blocks_01_alb_Layer1 = "blocks_01_alb_Layer1"

# Process: Make Feature Layer
arcpy.MakeFeatureLayer_management(blocks_01_alb_shp, blocks_01_alb_Layer1, "", "", "FID FID VISIBLE NONE;Shape Shape VISIBLE NONE;STATEFP10 STATEFP10 VISIBLE NONE;POP10 POP10 VISIBLE RATIO")

# Process: Intersect
arcpy.Intersect_analysis("blocks_01_alb_Layer1 #;USA_int_131111 #", State_01_int_7_shp, "ALL", "", "INPUT")

[T__WayneWhitley]

Ah, I think now I see the problem and I have just 2 points to make:
- When the describe and field info objects are created, it is apparently separate from the layer it was created from.
- There needs to be a better distinction between 'split' policy and 'ratio' policy in the webhelp - the ratio policy is what we are manipulating here; the split policy is associated for example with an attribute domain that governs how values are populated while editing, etc...so split policies set up in that way are not supported by geoprocessing tools (although I suppose you could mimic with the ratio policy).

Back to your script, the trick I believe is to 'stamp' your layer only after getting and modifying the field info object - below is your code I modified to work (hopefully):

inFeaturesAlb = "blocks_01_alb.shp"

# Create a describe object
desc = arcpy.Describe(inFeaturesAlb)

# Create a fieldinfo object
field_info = desc.fieldInfo

# Get the field index and set ratio policy
field_info.setSplitRule(field_info.findFieldByName("POP10"), "RATIO")  

layer = "temp_layer"
arcpy.MakeFeatureLayer_management(inFeaturesAlb, layer, '', '', field_info)

######################################################################    
#Do the intersection
#####################################################################     
intOut = "State_01_int_3.shp"    
inFeatures_int = [USA_layer, layer]

arcpy.Intersect_analysis(inFeatures_int, intOut)

Hope that clears up this issue a bit.


Wayne