Euclidean distance, allocation and other stuff...

804
3
10-03-2018 10:48 PM
Labels (1)
DanPatterson_Retired
MVP Esteemed Contributor
2 3 804

Your standard fare in raster world

Fill in some holes … 'buffer' a line... allocate space to stuff in space... identify unique regions

Fancy names.  Nibble, Euclidean distance, Euclidean allocation, Regiongroup

---- (1) The task ----

Start with a raster or array.  Fancy labels in the top left, some random-ish color scheme with values noted in the middle.

Now, zero ( 0 ) we will say is nodata.  The other numbers represent some class value. 

Fill in the 'gaps'... aka, nodata, with the value of a cell with data based on the closest distance.  Euclidean, crow flies, for our purposes, but it need not be. 

Go! 

What did you get for cells A04, H01 and H02?  What about cell D07 and H08?

---- (1) The reveal ----

Let's see how we did

A04 - 2   can't be 1 because the diagonal to '1' is 1.414*cell width, so 2 it is

H01 - 2   could have been a 2 or a 3, because both are 3 cells away from cells with values

H02 - 3   no brainer

D07 - 3   could have been 3 or 4, but 3 wins

H08 -3    3 is 2 cells away and 1 and 5 are a greater distance on an angle

---- (3) The distances ----

Pretty boring and obvious... but for completeness.

I don't do maps, but the dominant colors are 0 or root 2 as you can see from the spiffy ArcGIS Pro symbology tab.So no big surprises

 

---- (4) Allocation around linear features ----

Yes, that is possible too, sort of like a variable buffer, trade area, but currently just a simple Euclidean spatial allocation.

---- (5) The references ----

 Euclidean distance

Euclidean Allocation

NumPy/SciPy/  plus Arcpy stuff solution is what I used

import sys
import numpy as np
from scipy import ndimage as nd
from arcpy.geoprocessing import env
from arcpy.arcobjects import Point
from arcgisscripting import NumPyArrayToRaster, RasterToNumPyArray

env.overwriteOutput = True

def euclid_calc(a, mask_vals=0, dist=True, alloc=True):
    """Calculate the euclidean distance and/or allocation

    Parameters:
    a : array
        numpy float or integer array
    mask_vals : number, list or tuple
        If a single number is provided, a `mask` will be created using it.  A
        list or tuple of values can be used to provide multiple value masking.
    dist : boolean
        True, the distance of the closest non-masked value to the masked cell
    alloc : boolean
        True, the value of the closest non-masked value to the masked cell
    """

    if not dist:
        dist = None
    if not alloc:
        alloc = None
    m = np.isin(a, mask_vals).astype('int')
    dist, idx = nd.distance_transform_edt(m, return_distances=True,
                                          return_indices=True)
    alloc = a[tuple(idx)]
    return dist, alloc

# ----------------------------------------------------------------------
# .... final code section
if len(sys.argv) == 1:
    testing = True
    r_in = r"C:\GIS\A_Tools_scripts\Raster_tools\Data\poly.tif"
    a = RasterToNumPyArray(r_in)
    dist, alloc = euclid_calc(a, mask_vals=0, dist=True, alloc=True)
else:
    testing = False
    r_in = sys.argv[1]
    r_out = sys.argv[2]
    LLx = sys.argv[3]
    LLy = sys.argv[4]
    cell_sze = sys.argv[5]
    a = RasterToNumPyArray(r_in)
    dist, alloc = euclid_calc(a, mask_vals=0, dist=True, alloc=True)
    #
    # specify a, dist or alloc below... this could be a tool choice if needed
    #
    r0 = NumPyArrayToRaster(a, Point(LLx, LLy), cell_sze) 
    r0.save(r_out)‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍

Of course, when run from within the amazing Spyder python IDE, you can simply skip saving to output raster if needed, or try various other options from with the scipy slate other than the distance_transform_edt

On to other "special analyst" tools soon.

3 Comments
DouglasLong
New Contributor III

This is quite interesting Dan, I have not written my own Euclidean Allocation but I have my own Euclidean Distance and direction. we will have to chat more as I do a lot of this on a daily basis and would be beneficial to pick each other's brains.

I use a pretty simple method of getting a range of the x and y coordinates and I use numpy broadcasting to get the distance from all the others relative to a source cell. I posted an example of a possible output image. with the source being the little green dot.

DanPatterson_Retired
MVP Esteemed Contributor

Nice!  Sometimes I just use meshgrid to emulate cell locations.

c = np.array([5,5])

shp = 11

m = np.meshgrid(np.arange(shp), np.arange(shp))

xy = np.array(list(zip(m[0].ravel(), m[1].ravel())))

c = np.array([5,5])

e_dist(c, xy).reshape(shp, shp)

array([[7.07, 6.4 , 5.83, 5.39, 5.1 , 5.  , 5.1 , 5.39, 5.83, 6.4 , 7.07],
       [6.4 , 5.66, 5.  , 4.47, 4.12, 4.  , 4.12, 4.47, 5.  , 5.66, 6.4 ],
       [5.83, 5.  , 4.24, 3.61, 3.16, 3.  , 3.16, 3.61, 4.24, 5.  , 5.83],
       [5.39, 4.47, 3.61, 2.83, 2.24, 2.  , 2.24, 2.83, 3.61, 4.47, 5.39],
       [5.1 , 4.12, 3.16, 2.24, 1.41, 1.  , 1.41, 2.24, 3.16, 4.12, 5.1 ],
       [5.  , 4.  , 3.  , 2.  , 1.  , 0.  , 1.  , 2.  , 3.  , 4.  , 5.  ],
       [5.1 , 4.12, 3.16, 2.24, 1.41, 1.  , 1.41, 2.24, 3.16, 4.12, 5.1 ],
       [5.39, 4.47, 3.61, 2.83, 2.24, 2.  , 2.24, 2.83, 3.61, 4.47, 5.39],
       [5.83, 5.  , 4.24, 3.61, 3.16, 3.  , 3.16, 3.61, 4.24, 5.  , 5.83],
       [6.4 , 5.66, 5.  , 4.47, 4.12, 4.  , 4.12, 4.47, 5.  , 5.66, 6.4 ],
       [7.07, 6.4 , 5.83, 5.39, 5.1 , 5.  , 5.1 , 5.39, 5.83, 6.4 , 7.07]])

And a helper function that I have posted before (part of my inprogress, arraytools on GitHub) is below.

import numpy as np  # of course

def e_dist(a, b, metric='euclidean'):
    """Distance calculation for 1D, 2D and 3D points using einsum

    `a`, `b` : array like
        Inputs, list, tuple, array in 1, 2 or 3D form
    `metric` : string
        euclidean ('e', 'eu'...), sqeuclidean ('s', 'sq'...),
    -----------------------------------------------------------------------
    """

    a = np.asarray(a)
    b = np.atleast_2d(b)
    a_dim = a.ndim
    b_dim = b.ndim
    if a_dim == 1:
        a = a.reshape(1, 1, a.shape[0])
    if a_dim >= 2:
        a = a.reshape(np.prod(a.shape[:-1]), 1, a.shape[-1])
    if b_dim > 2:
        b = b.reshape(np.prod(b.shape[:-1]), b.shape[-1])
    diff = a - b
    dist_arr = np.einsum('ijk,ijk->ij', diff, diff)
    if metric[:1] == 'e':
        dist_arr = np.sqrt(dist_arr)
    dist_arr = np.squeeze(dist_arr)
    return dist_arr

at may seem like overkill for small arrays, but I have used it for calculating distance from an origin to multiple destinations in the millions in under a second..  

It scales well also, array a and b can be of any shape, so it emulates SciPy's pdist, cdist and the other one... 

hence the e_dist (an homage to Einstein and einsum )

DanPatterson_Retired
MVP Esteemed Contributor

So if you have some pure raster, please a great idea to share and collaborate

About the Author
Retired Geomatics Instructor at Carleton University. I am a forum MVP and Moderator. Current interests focus on python-based integration in GIS. See... Py... blog, my GeoNet blog...
Labels