I have a shapefile with hundreds of caribou movement corridors. using this data, I am trying to create a map displaying the frequency of use to determine important movement and migration pathways for different herds.

The original shapefile has several years of data of individual collared caribou movements derived using Brownian Bridge Movement Modelling. From this file, I have successfully created my desired outcome by: 1. Using the union tool 2. Multipart to singlepart 3. Spatial join to itself 4. Change shapefile symbology. Example of intended final product The problem I am now encountering is that I have extended my study dates and included several years into a single model. When I try to follow the above steps with the large dataset, I a thrown an error code 'insufficient memory' at the spatial join step (3.). Clearly, my computer does not have enough ram to run this calculation, so I have been trying to divide my shapefile into 4 equal parts, using a fishnet grid, which I plan to merge back together after I run the spatial join on each parcel. I have successfully created the grid and again using the Union tool, created 4 separate shapefiles of data. Now, when I try to 'split' the caribou polygons from the grid attribute, the operation fails. All polygons with fishnet grid prior to splitting of grids I understand that ArcGIS pro and perhaps a patch for arcmap could help my computer access more of its ram to run such operation, however, I do not currently have access to these, so am looking for some kind of alternative method of achieving my desired results.

I will try to add some images shortly.

  • What happens when the split operation fails? Do you get an error message? – csk Apr 5 '19 at 17:52
  • Thanks for responding. No there is no error message, it just times out after about an hour of processing, and in the Results Status it simply says 'failed'. I should also add, that I was unable to use the Union tool in arcamp to divide the polygons into equal parts. Instead I used QGIS for this particular action. – user62765 Apr 5 '19 at 18:01
  • Have a look at this question, for a possible solution in QGIS that doesn't require any geoprocessing at all. You just use blending modes to create a heatmap. – csk Apr 5 '19 at 18:34
  • Try running this operation using ArcGIS Pro, which utilizes 64 bit architecture. Please let me know if you can successfully run it on that platform. – Aaron Apr 5 '19 at 18:34
  • 1
    You can install [64 bit geoprocessing] ( desktop.arcgis.com/en/arcmap/10.3/analyze/executing-tools/…) for arcmap but when you do you run the tools in the background. – Hornbydd Apr 5 '19 at 19:05

I think the simplest way to solve that is divide polygons into non-overlapping sets, converting individual sets into rasters and use cell statistics on these rasters.

So, run polygon neighbours tool on your corridors:

arcpy.PolygonNeighbors_analysis(in_features="SUBCATCHMENTS", out_table="C:/SCRATCH/SCRATCH.gdb/NEIGHBOURS", in_fields="", area_overlap="AREA_OVERLAP", both_sides="NO_BOTH_SIDES")

making sure that "Include area overlap" is checked and run script below from mxd. Script creates new field "PART_NO" and populates it. First parameter of the script is your polygons layer (shapefile!), second parameter is polygon neighbours output.

creates non-adjacent groups of polygons
import arcpy
import networkx as nx

infc = arcpy.GetParameterAsText(0)
table = arcpy.GetParameterAsText(1)
for f,t in fromto:
d = nx.coloring.greedy_color(G, strategy=nx.coloring.strategy_largest_first)
arcpy.AddField_management(infc, "PART_NO","Short")
with arcpy.da.UpdateCursor(infc,("FID","PART_NO")) as cursor:
    for row in cursor:

Output shows that original layer can be presented as 5 for non-adjusting sets. There are multiple overlaps in your data, so expect slightly more than 5 groups to deal with:

enter image description here

import os
import sys
import arcpy

arcpy.env.overwriteOutput = True
# import pythonaddins

#######User Selection 1
poly_lyr = arcpy.GetParameterAsText(0)
# poly_lyr = r"C:\Users\EB\Desktop\polygons\Export_OutputPro.shp"

#######User Selection 2
# num_out_polys = 10
num_out_polys = int(arcpy.GetParameterAsText(1))
#map units (eg meters) and the difference in area between the largest and the smallest polygons
#0.005 - 0.02%; 0.01 - 0.03%; 0.05 - 0.1%; 0.1 - 0.3%;
step_value = 1

#######User Selection 3
# orientation = 'NS' #'WE' / 'NS'
orientation = arcpy.GetParameterAsText(2)

if orientation != 'NS' or orientation != 'WE':
    orientation == 'NS'

# number of splits
splits = [round(float(100)/float(num_out_polys), 2)] * num_out_polys

#spatial reference of the output fc will be of the polygon layer
sr = arcpy.SpatialReference(arcpy.Describe(poly_lyr).spatialReference.factoryCode)

#source polygon fields
fields = [f.name for f in arcpy.ListFields(poly_lyr) if not f.required]

if int(arcpy.GetCount_management(poly_lyr).getOutput(0)) != 1:
    arcpy.AddMessage('Need to have exactly one feature selected', 'Error')

#get polygon geometry and extent property
with arcpy.da.SearchCursor(poly_lyr, fields + ["SHAPE@"]) as cur:
    for row in cur:
        attributes = list(row[:-1])
        polygon = row[-1]
        extent = polygon.extent

#orient lines either North-South (up-down) or West-East (left to right)
if orientation == 'NS':
    x_max = extent.XMax + step_value
    x_min = extent.XMin + step_value
    y_max = extent.YMax
    y_min = extent.YMin

if orientation == 'WE':
    x_max = extent.XMax
    x_min = extent.XMin
    y_max = extent.YMax - step_value
    y_min = extent.YMin

cut_poly = polygon
# name of output shapefile
outputshape_name = arcpy.GetParameterAsText(3)
#output feature class create/clean
# arcpy.env.scratchGDB
mem_path = os.path.join(arcpy.Describe(poly_lyr).path, str(outputshape_name)+".shp")
if arcpy.Exists(mem_path):
mem = arcpy.CopyFeatures_management(poly_lyr, mem_path)

lines = []

with arcpy.da.InsertCursor(mem, fields + ["SHAPE@"]) as icur:
    for i in splits[:-1]: #need to get all but the last item
        tolerance = 0
        while tolerance < i:
            pnt_arr = arcpy.Array()
            if orientation == 'NS':
                #construct North-South oriented line
                pnt_arr.add(arcpy.Point(x_min, y_max))
                pnt_arr.add(arcpy.Point(x_min, y_min))

            if orientation == 'WE':
                #construct West-East oriented line
                pnt_arr.add(arcpy.Point(x_min, y_max))
                pnt_arr.add(arcpy.Point(x_max, y_max))

            line = arcpy.Polyline(pnt_arr, sr)

            #cut polygon and get split-parts
            cut_list = cut_poly.cut(line)
            if orientation == 'NS':
                tolerance = 100 * cut_list[1].area / polygon.area
                x_min += step_value

            if orientation == 'WE':
                tolerance = 100 * cut_list[0].area / polygon.area
                y_max -= step_value

        # part 0 is on the right side and part 1 is on the left side of the cut
        if orientation == 'NS':
            cut_poly = cut_list[0]
            icur.insertRow(attributes + [cut_list[1]])

        if orientation == 'WE':
            cut_poly = cut_list[1]
            icur.insertRow(attributes + [cut_list[0]])

    #insert last cut remainder
    if orientation == 'NS':
        icur.insertRow(attributes + [cut_list[0]])

    if orientation == 'WE':
        icur.insertRow(attributes + [cut_list[1]])

#for illustration purposes only
arcpy.CopyFeatures_management(lines, 'in_memory/lines')

#evaluation of the areas error
done_polys = [f[0] for f in arcpy.da.SearchCursor(mem_path, 'SHAPE@AREA')]

#the % of the smallest and the largest areas
# arcpy.AddMessage("{} Precision error".format(round(100 - 100 * (min(done_polys) / max(done_polys)), 2)))
arcpy.AddWarning("OutPut Path is :"+mem_path)

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.