0

I'm trying to run a loop based on a SearchCursor that selects layers from an input polygon and polyline. I want to select a feature from the polygon to use as a mask for running Euclidean Allocation on selected lines within the polygon feature. The goal is to have multiple areas of Euclidean Allocation run on a subset of polylines, each masked by a corresponding polygon.

Here's my original code. Note that this function is part of a much larger, otherwise functional script.

The input parameters for the three functions below are inline = a polyline feature class, inpoly = a polygon feature class, and valfield = a single field which corresponds to both inline and inpoly: i.e., each value of valfield = x in inline corresponds to lines that fall within inpoly features of the same value. All of these parameters are derived from in-memory feature classes fed in by earlier functions in the script--they're not located directly on disk.

def eucal_by_poly(inline, inpoly, valfield): polylist = [] with arcpy.da.SearchCursor(inpoly, ("OID@", "SHAPE@", valfield)) as cursor: vallist = sorted({row[2] for row in cursor}) for v in vallist: lyr = "lyr" selexp = "{0} = {1}".format(arcpy.AddFieldDelimiters(lyr, valfield), v) mask_lyr = arcpy.MakeFeatureLayer_management(inpoly, lyr, selexp) arcpy.env.mask = mask_lyr line_lyr = arcpy.MakeFeatureLayer_management(inline, lyr, selexp) eucalsub = arcpy.sa.EucAllocation(line_lyr, "", "", "", "LAT_ID") eucalsubpoly = arcpy.RasterToPolygon_conversion(eucalsub, "eucalsubpoly", "SIMPLIFY", "Value") polylist.append(eucalsubpoly) polyeucal = arcpy.Merge_management(polylist, "polyeucal") del cursor return polyeucal

When I run this code I get this error:

Traceback (most recent call last): File "G:\tgf_python\tools\py_scripts_working\transect_working_rasversion.py", line 406, in <module> polyeucal = eucal_by_poly(cliptrans, eucal_poly, IDfield) File "G:\tgf_python\tools\py_scripts_working\transect_working_rasversion.py", line 397, in eucal_by_poly arcpy.env.mask = mask_lyr File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\geoprocessing\_base.py", line 529, in set_ self[env] = val File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\geoprocessing\_base.py", line 581, in __setitem__ ret_ = setattr(self._gp, item, value) RuntimeError: Object: Error in accessing environment <mask>

I've tried several variations, including using arcpy.env.mask = "SHAPE@", to no avail. I also tried saving the poly layer to a real file on disk, both as a feature class and a rasterized version. In all cases I get the same "Error in accessing environment" message when trying to apply the mask.

Does anyone have a solution? I'm using ArcGIS 10.1 with an Advanced License and Python 2.7 on Windows 7. (Any other suggestions on improving code are welcome.)

  • Please extract your function out into a code snippet that works up to where you are stuck. There is an edit button beneath your question that you can use to improve it in this way. – PolyGeo Jan 30 '16 at 23:06
  • Sorry--I'm not sure how to make it more clear. The specific line where the code fails is: 'arcpy.env.mask = mask_lyr' Don't you need to see the entire function to understand the context? Please clarify how I can improve the question--I would really like some solutions. – qui_bono Jan 31 '16 at 23:59
  • You are feeding three parameter values inline, inpoly and valfield into your function so just provide hardwired values for them and a description of the data that they represent so that we can use that to try and test. Our volunteers are usually happy to try and help debug code snippets but your code and its context is your responsibility. – PolyGeo Feb 1 '16 at 0:07
  • I get that it's my responsibility--just wasn't sure what was needed. I've updated above with further explanation. Let me know if more is needed. – qui_bono Feb 1 '16 at 2:08
  • I haven't gotten any responses, so I'm wondering if I'm still posting my code incorrectly. Any advice? – qui_bono Feb 10 '16 at 2:09
0

My previous answer was rushed, so I’d like to try again, now that I’ve had some time to reconsider.

If each polygon intersects only one line, the simplest strategy is an Intersect:

def linesToPoly(inLines, inPoly, outName):
        result = arcpy.Intersect_analysis([inLines, inPoly], outName)
        return result

If multiple lines intersect each polygon, then an interpolator-type tool is needed to create value-areas around the lines. The major issue is that none of the standard interpolators (IDW, Spline, EBK, etc.) take lines as inputs, only points. I’m guessing that’s why Euclidian Allocation was used, because it takes lines as inputs and performs an interpolation-like operation. After playing with some mock data, I also discovered that Euclidian Allocation behaves badly (for reasons unknown to me) with line-file inputs covering larger areas. I assume this is also why masks were used to make the area smaller, instead of running the Euclidean Allocation tool all in one go.

While the original approach might work, but it’s likely pretty slow, because it runs Euclidian Allocation again and again for every polygon-group. There might also be trouble if any polygons overlap. Merge doesn’t have conflict resolution if overlapping areas have different values, as calculated by the interpolator. Finally, running an interpolator inside of a mask rather than globally means that it lacks awareness of nearby neighbors, which has notable effects at the borders. Some of these effects might be desired, but I would take a different approach to solve this problem.

If lines are the major stumbling block, then my first step would be to make points out of lines. There are multiple ways to accomplish this, but I used Densify + Create Points from Vertices. Other methods are discussed in following link:

Creating point file along line (trail) at regular distances (0.20 miles) using ArcGIS Desktop?

With points, any interpolator can be used, not just Euclidian Allocation. I also found that while Euclidian Allocation performs poorly with large area line files, it works amazingly well with points. So, run any interpolator on the point file, convert its output to raster, and intersect with the original polygons:

def EAToPoly(inLines, inPoly, valfield):
        arcpy.Densify_edit(inLines, "DISTANCE", 100)
        arcpy.FeatureVerticesToPoints_management(inLines, "NewPoints", "ALL")
        arcpy.sa.EucAllocation("NewPoints", "", "", "", valfield, "", "EAOutRaster")
        arcpy.RasterToPolygon_conversion("EAOutRaster", "EAPoly", "SIMPLIFY", valfield)
        result = arcpy.Intersect_analysis(["EAPoly", inPoly], "EAPoly_intersect")
        return result

Note: no loops, no masks, no make feature layers, no merge, no cursors, no set comprehensions, and no conflicts with overlapping areas, all the while incorporating awareness of neighbors. And, you can substitute ANY interpolator method in the place of Euclidian Allocation.

Here’s a graphic of my output with Euclidian Allocation used as the interpolator. Pay attention especially to the handling of overlapping polygons, and the influence of nearby areas. enter image description here

  • Thanks for all the input. My company ended up implementing an FME-based solution for this problem but these responses, especially Priscilla's, were great and helped me modify existing code for other uses. – qui_bono May 16 '16 at 15:14
  • I'm so glad it helped. I remember spending a good chunk of time on this 2nd answer, after submitting such a bad 1st answer. :) – Priscilla May 19 '16 at 20:33
1

This might not be your entire problem, but MakeFeatureLayer doesn't return the layer name - you've already provided it. I don't know that it has a defined return value, so you're trying to set the mask environment variable to an unknown (from our perspective) value. Instead, you want to set it to the layer name, so change that line to arcpy.env.mask = lyr. That might not be the only issue if it's telling you it can't access the mask environment variable, but start there.

-1

This problem would be clearer if we could get a sense of your spatial problem, not just code. What are you trying to get arcgis to do? With only lines and polygons, I have the sense that you're probably using the wrong tools for the job.

Euclidean allocation usually takes rasters, but why do you even need rasters here with lines and polygons? With lines and polygons, you should probably be using clip, spatial select, or intersect, and not bother with masks. The intersect tool with a join might eliminate the need for this entire function, but I can't say because you haven't articulated your spatial problem.

Merge tool is incorrectly applied. Merge is used to combine two datasets of the same type. You're combining a python list with a... I'm not sure what because "polyeucal" doesn't even exist before you call it with Merge. Even if "polyeucal" was a real feature class, you still can't use Merge to combine it with a python list. That's more of a field calculator operation.

I think maybe you should take a step back and get a better understanding of some fundamental arcgis tools to solve this problem. Here are some helpful links:

Euclidean Allocation: http://pro.arcgis.com/en/pro-app/tool-reference/spatial-analyst/euclidean-allocation.htm

Clip: http://pro.arcgis.com/en/pro-app/tool-reference/analysis/clip.htm

Spatial Select: http://desktop.arcgis.com/en/arcmap/10.3/tools/data-management-toolbox/select-layer-by-location.htm

Intersect: https://pro.arcgis.com/en/pro-app/tool-reference/analysis/intersect.htm

Merge: http://pro.arcgis.com/en/pro-app/tool-reference/data-management/merge.htm

Calculate Field: https://pro.arcgis.com/en/pro-app/tool-reference/data-management/calculate-field.htm

  • 2
    arcpy.Merge_management can accept a list, so long as the list contains feature classes of the same type (which it looks like it does in this case). – nicksan Feb 11 '16 at 0:40
  • @nicksan You just blew my mind a bit. I didn't realize that lists could contain feature classes. I originally interpreted the list as holding the values from the polygons, not the polygons themselves. I see this problem differently now. – Priscilla Feb 11 '16 at 4:00
  • just to be a little clearer than I was before - it should be a list of paths to feature classes. I think in this case, it's still not happening properly because RasterToFeatureClass_conversion doesn't output the name as its return value either, but in theory it works (and I've done it). I have stored feature data directly in lists too, after reading them directly from the shape field, but I'm not certain that merge would work on them (I suspect it would not). – nicksan Feb 11 '16 at 20:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.