Is it possible to divide a raster into equal parts based off the sum of its cells?
Hi Brandon Jumonville ,
This is what I came up with...
For starters, having the REST url did not give me access to the actual data, so you might have to rethink part of what I am going to share.
So I captures a screenshot, georeferenced with using the base map at Baton Rouge 
Then I have to do some post processing and resampled values ranging from 1 to 10 corresponding (not the range of values self) to your original screenshots.
What I did was converting the raster to points, so each point will have the "time" value:
Then I picked 4 locations to start the zones:
The next step was to use Generate Near Table—Help | ArcGIS Desktop to generate a table that has the distance between each zone point and all the raster point (time values).
I used this in the following script as input:
def main():
import arcpy
# points (converted from raster, each cell is a point)
fc_pnt = r'C:\GeoNet\ZonasTimeCalls\data.gdb\points'
fld_oid = 'OID@'
fld_value = "grid_code"
fld_zone = "Zone"
dct_pnts = {r[0]: r[1] for r in arcpy.da.SearchCursor(fc_pnt, (fld_oid, fld_value))}
max_val = max(dct_pnts.values())
print max_val
zone_score = sum(dct_pnts.values()) / 4 + max_val
print zone_score
# near table (created with generate near table, using 4 points and the raster points)
tbl = r'C:\GeoNet\ZonasTimeCalls\data.gdb\NearTable'
fld_id_zone = "IN_FID"
fld_id_pnt = "NEAR_FID"
fld_dist = "NEAR_DIST"
# create a list of zones
zones = list(set([r[0] for r in arcpy.da.SearchCursor(tbl, (fld_id_zone))]))
print zones
# create list of id pnt for each zone, store in nested list, sorted by near dist
dct_zone_value = {}
for zone in zones:
flds = (fld_id_pnt, fld_dist)
where = "{0} = {1}".format(fld_id_zone, zone)
lst = [[r[0], r[1]] for r in arcpy.da.SearchCursor(in_table=tbl,
field_names=flds,
where_clause=where,
sql_clause=(None, 'ORDER BY {0} ASC'.format(fld_dist)))]
dct_zone_value[zone] = lst
dct_scores = {zone: 0 for zone in zones}
## for zone, lst in dct_zone_value.items():
## print zone, len(lst)
dct_pnt_zone = {}
print dct_scores
items = len(dct_zone_value[zones[0]])
print items
for i in range(items):
if i % 1000 == 0:
print "Processing:", i
for zone in zones:
lst = dct_zone_value[zone][i]
id_pnt = lst[0]
if not id_pnt in dct_pnt_zone:
# dist = lst[1]
val = dct_pnts[id_pnt]
score = dct_scores[zone]
if score + val < zone_score:
dct_pnt_zone[id_pnt] = zone
dct_scores[zone] += val
# print dct_pnt_zone
print dct_scores
# add field to points:
arcpy.AddField_management(fc_pnt, fld_zone, "SHORT")
flds = ('OID@', fld_zone)
cnt = 0
with arcpy.da.UpdateCursor(fc_pnt, flds) as curs:
for row in curs:
cnt += 1
if cnt % 1000 == 0:
print "Updating point:", cnt
oid = row[0]
if oid in dct_pnt_zone:
row[1] = dct_pnt_zone[oid]
curs.updateRow(row)
if __name__ == '__main__':
main()... and this was the result in the console window:
10
10236
[1, 2, 3, 4]
{1: 0, 2: 0, 3: 0, 4: 0}
15991
Processing: 0
Processing: 1000
Processing: 2000
...
Processing: 14000
Processing: 15000
{1: 10202, 2: 10235, 3: 10235, 4: 10235}
Updating point: 1000
Updating point: 2000
...
Updating point: 14000
Updating point: 15000So the scores of each zone (sum of the point values) are quite similar. And when you draw the points you will see this:
As you can see the red zone extends to the south and that is not very useful. Perhaps a score should be included that weighs the distance to get a better result.
More later...