3

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.

Improve this question
8
  • What happens when the split operation fails? Do you get an error message? Commented Apr 5, 2019 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. Commented Apr 5, 2019 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. Commented Apr 5, 2019 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. Commented Apr 5, 2019 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. Commented Apr 5, 2019 at 19:05

2 Answers 2

Reset to default
1

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) fromto=arcpy.da.TableToNumPyArray(table,("src_FID","nbr_FID")) G=nx.Graph() for f,t in fromto: G.add_edge(f,t,weight=1) 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: row[1]=d[row[0]]+1 cursor.updateRow(row) arcpy.AddMessage("Done")

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

Improve this answer
-1 
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') sys.exit(0) #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): arcpy.Delete_management(mem_path) mem = arcpy.CopyFeatures_management(poly_lyr, mem_path) arcpy.DeleteFeatures_management(mem) 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) lines.append(line) #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)
Improve this answer

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.