topic Re: Is it possible to divide a raster into equal parts based off the sum of its cells? in Public Safety Questions

Is it possible to divide a raster into equal parts based off the sum of its cells?

Anonymous User — Thu, 09 Nov 2017 13:17:23 GMT

I'm working on a project to redraw district boundaries for our police department to help spread the work load more evenly among officers.

I created a raster layer showing overall time spent on calls across the city and was wondering if there is a way to divide the layer into equal value amounts rather than equal areas.

We have 4 districts and I'd like to divide my raster so that the amount of time (sum of the cells) is equally split. The boundaries have been arbitrarily drawn throughout the years as new areas where annexed so completely redrawing the district lines is not a problem.

Any ideas on how to approach this or am I hoping for too much?

Re: Is it possible to divide a raster into equal parts based off the sum of its cells?

DanPatterson_Retired — Thu, 09 Nov 2017 13:30:24 GMT

Totally untried... but what about Cost Allocation which would be different than

Euclidean Allocation .... suggesting that a cost surface might apportion the space differently than purely distance based.

Re: Is it possible to divide a raster into equal parts based off the sum of its cells?

XanderBakker — Thu, 09 Nov 2017 13:55:59 GMT

Interesting thought to use Cost Distance for this purpose. Will probably not result in areas with an equal sum of the value, but might come close and requires the user to define the starting points.

I have done something (remotely) related to this (vector based) where the objective was to select as many parcels as possible based on a minimum distance between them. A similar approach could be used. Would be interesting to take a look. Brandon Jumonville , would you be able to share (a part of) the time raster or other input data to have a look at what is possible?

Re: Is it possible to divide a raster into equal parts based off the sum of its cells?

Anonymous User — Thu, 09 Nov 2017 19:19:21 GMT

Re: Is it possible to divide a raster into equal parts based off the sum of its cells?

Anonymous User — Thu, 09 Nov 2017 19:24:38 GMT

So this is where I'm at right now. I took the total calls for a 5.75 year period, determined the average time officers spent per incident type, applied that value to the calls and created a 250'X250' raster from the points (1.25 million if you were wondering). So now I have these time values ranging from .36 to 4076.57 hours per cell but they aren't evenly divided between the current district boundaries... that's where I'd like to redraw the boundaries to evenly spread, at least reasonably close, the time spent between four districts.

Re: Is it possible to divide a raster into equal parts based off the sum of its cells?

XanderBakker — Thu, 09 Nov 2017 19:25:21 GMT

The image provide some more information on what you have, but is it possible to share some of the actual data?

Re: Is it possible to divide a raster into equal parts based off the sum of its cells?

Anonymous User — Thu, 09 Nov 2017 19:57:58 GMT

Here's a link to a tile service on our portal

Services Directory - Average_CAD_Call_Total_Time(MapServer)

Re: Is it possible to divide a raster into equal parts based off the sum of its cells?

XanderBakker — Sat, 11 Dec 2021 06:10:53 GMT

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: 15000‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

So 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...

Re: Is it possible to divide a raster into equal parts based off the sum of its cells?

XanderBakker — Fri, 10 Nov 2017 22:34:52 GMT

Let's tag Dan_Patterson‌ and see if he knows of some numpy or other python magic that he can toss in.

Re: Is it possible to divide a raster into equal parts based off the sum of its cells?

DanPatterson_Retired — Fri, 10 Nov 2017 22:40:31 GMT

Only to add, that if they don't have the license to use Near or its table equivalent,

generate-near-table-for-points

can be used, with a bit of upfront work particularly for a large . There is no license restriction on it, but some coding familiarity

Re: Is it possible to divide a raster into equal parts based off the sum of its cells?

XanderBakker — Sat, 11 Nov 2017 21:03:14 GMT

Just thought of another alternative to do this calculation by using location - allocation in ArcGIS Online (or using Network Analyst in Desktop). See below the description (with Spanish interface, sorry) of how to do this in ArcGIS Online:

Now before you start, you should publish the 4 locations and the demand point (raster values) to ArcGIS Online. Note that ArcGIS Online is limited to 10.000 demand points to be included in the analysis. I reduced the points filtering out the lowest grid_code values reducing the 16.000 points to a little over 6.000 points.

In ArcGIS Online, load the layers a a new map:

Switch on analysis and select "Choose best facilities:

Next choose "Maximize coverage with capacity":

Then specify the way you want to assign the demand to the facilities (in this case drive time), select the option installations to demand and select the demand layer which are the raster points with the grid_code:

Next you will have to select "Field" and the field "grid_code" in de demand layer, to use the time value for assigning the demand per point. Also indicate that the drive range is unlimited or specify a valid value. Make sure that the layer with the 4 locations is selected at number 4:

Before the next stop, you will have to indicate the maximum demand for each of the 4 locations. This you need to defined based on the statistics of the grid_code field. See below the stats (values ranging from 2 to 10 and the sum is 31235):

The value 31235 is the total demand. You should divide this by 4 (4 locations) and add the maximum value of the range to it as security to not keep any demand point out (±7820). Specify that at number 4 as constant capacity and you are ready to do the analysis:

Always click on the link "show credits" te review the number of credits this operation will consume. It will be the number of demand points divided by 10:

I did not execute the operation, since I didn't want to consume 631.9 credits just to see if this works... So no guarantees at this point that the desired result is actually created.

Will try with a much smaller dataset and check in later to share the result.

Re: Is it possible to divide a raster into equal parts based off the sum of its cells?

XanderBakker — Sat, 11 Nov 2017 21:15:58 GMT

So I did a small test with about 500 demand points, filtering in Desktop the points by "grid_code > 4 AND Mod(OBJECTID, 10) = 0" and published that layer as service:

And did the process. However, for some reason it only processed one candidate and assigned about 132 demand points, which also had a strange logic applied:

I guess it's back to the drawing board...

Re: Is it possible to divide a raster into equal parts based off the sum of its cells?

XanderBakker — Sat, 11 Nov 2017 21:57:59 GMT

Today I ran in to a tool that was mentioned at another thread called "Allocation Tools for ArcGIS" shared by a user called "deelesh": http://www.arcgis.com/home/item.html?id=c5ae481f25604d1bb050d9bd1d6e3c06 which allows for location allocation using an open source Linear Programming (LP) solver called COIN.

So I thought I had a go at it, but it requires a network dataset. Since I don't have one, I thought I could create something similar using the Spanning Tree Tools created by Dan Patterson and which you can download here: https://www.arcgis.com/home/item.html?id=6ce9db93533345e49350d30a07fc913a

I ran the Spanning Tree tool on my demand points (a simplified set) and this was the result:

Then I turned it into a network and used it in the location allocation tool:

And this is the result:

Still not perfect, but we are getting somewhere...

Re: Is it possible to divide a raster into equal parts based off the sum of its cells?

DanPatterson_Retired — Sat, 11 Nov 2017 22:06:47 GMT

Cool! Looks good!

Re: Is it possible to divide a raster into equal parts based off the sum of its cells?

XanderBakker — Sun, 12 Nov 2017 18:20:17 GMT

I changed the code to make the growth of each zone depend on the accumulated value. The zone with the lowest sum value gets to pick the next demand point (raster points). This made the result a bit better, but it is still not what I would accept (left before and right after the code change):

Re: Is it possible to divide a raster into equal parts based off the sum of its cells?

XanderBakker — Sun, 12 Nov 2017 18:34:45 GMT

Since it looked cool as Dan mentioned, I thought that I would do another test, creating the network using Dan's tool but now on all the nearly 16.000 points and create the network. This is a detail of how that looks like:

... I did the analysis on the entire set:

Which resulted in something that actually makes sense:

There may still be some manual adjustments required, but for now it seems to be an acceptable result.

Re: Is it possible to divide a raster into equal parts based off the sum of its cells?

DanPatterson_Retired — Sun, 12 Nov 2017 18:47:03 GMT

Excellent MST has purpose ! A big bookmark on this

Re: Is it possible to divide a raster into equal parts based off the sum of its cells?

XanderBakker — Sun, 12 Nov 2017 20:00:42 GMT

It is a great tool Dan, thanks for sharing and especially in my line of work (utilities) the tool can be used in many different planing and design activities.

Re: Is it possible to divide a raster into equal parts based off the sum of its cells?

DanPatterson_Retired — Sun, 12 Nov 2017 20:45:06 GMT

given me the goal to implement weighted MST over the holidays

Re: Is it possible to divide a raster into equal parts based off the sum of its cells?

Anonymous User — Mon, 13 Nov 2017 12:50:34 GMT

Very cool!