Skip navigation
All People > Dan_Patterson > Py... blog
1 2 3 Previous Next

Py... blog

85 posts

I rarely get a great idea.  This is another one of those cases.  But I do recognize a great idea when I see it.  Sadly I usually forget it, then rediscover it many months or years later.  I spotted a fantastic idea today and I had to drop everything (between answering questions) and jot it down.

 

The premise that twigged my interest was about how to 'combine' values of elevation, slope and aspect on Stack Overflow.  I did a wowwser!  Terrain derivatives on Stack!  I was enthralled.  The author.... unutbu ... of mysterious persona, put together a short missive using built-in numpy tools.  The code is short and sweet, works with floating point or integer rasters (aka, arrays) and is fast (at least in my testing so far).

 

Before I ruin the rest of my day with some droll corner case I hadn't thought of, I will jot this down now for posterity.

 

The large print

  • This doesn't require a Spatial Analyst License
  • A complete implementation (coming soon) uses arcpy's RasterToNumPyArray and NumPyArrayToRaster to communicate between the thin veil separating Arc* from the array world
  • This is in the initial stages of development of a full-fledged tool, so hold you ...but!... comments

 

Start with 3 small rasters/arrays.  This way, you can see how the system works.

 

Array 'a'Array 'b'
array([[0, 0, 0, 4, 4, 4, 1, 1, 1],
       [0, 0, 0, 4, 4, 4, 1, 1, 1],
       [0, 0, 0, 4, 4, 4, 1, 1, 1],
       [2, 2, 2, 1, 1, 1, 2, 2, 2],
       [2, 2, 2, 1, 1, 1, 2, 2, 2],
       [2, 2, 2, 1, 1, 1, 2, 2, 2],
       [1, 1, 1, 4, 4, 4, 0, 0, 0],
       [1, 1, 1, 4, 4, 4, 0, 0, 0],
       [1, 1, 1, 4, 4, 4, 0, 0, 0]])
array([[0, 0, 0, 1, 1, 1, 2, 2, 2],
       [0, 0, 0, 1, 1, 1, 2, 2, 2],
       [0, 0, 0, 1, 1, 1, 2, 2, 2],
       [3, 3, 3, 4, 4, 4, 5, 5, 5],
       [3, 3, 3, 4, 4, 4, 5, 5, 5],
       [3, 3, 3, 4, 4, 4, 5, 5, 5],
       [0, 0, 0, 1, 1, 1, 2, 2, 2],
       [0, 0, 0, 1, 1, 1, 2, 2, 2],
       [0, 0, 0, 1, 1, 1, 2, 2, 2]])

 

Array 'c'unique combinations
array([[0, 0, 0, 3, 3, 3, 0, 0, 0],
       [0, 0, 0, 3, 3, 3, 0, 0, 0],
       [0, 0, 0, 3, 3, 3, 0, 0, 0],
       [1, 1, 1, 4, 4, 4, 1, 1, 1],
       [1, 1, 1, 4, 4, 4, 1, 1, 1],
       [1, 1, 1, 4, 4, 4, 1, 1, 1],
       [2, 2, 2, 5, 5, 5, 2, 2, 2],
       [2, 2, 2, 5, 5, 5, 2, 2, 2],
       [2, 2, 2, 5, 5, 5, 2, 2, 2]])
array([[0, 0, 0, 6, 6, 6, 1, 1, 1],
       [0, 0, 0, 6, 6, 6, 1, 1, 1],
       [0, 0, 0, 6, 6, 6, 1, 1, 1],
       [2, 2, 2, 7, 7, 7, 3, 3, 3],
       [2, 2, 2, 7, 7, 7, 3, 3, 3],
       [2, 2, 2, 7, 7, 7, 3, 3, 3],
       [4, 4, 4, 8, 8, 8, 5, 5, 5],
       [4, 4, 4, 8, 8, 8, 5, 5, 5],
       [4, 4, 4, 8, 8, 8, 5, 5, 5]],

 

So here it is.

def combine_(*arrs):
    """Combine arrays to produce a unique classification scheme
    :
    :References:
    :-----------
    : https://stackoverflow.com/questions/48035246/
    :       intersect-multiple-2d-np-arrays-for-determining-zones
    : original def find_labels(*arrs):
    """

    indices = [np.unique(arr, return_inverse=True)[1] for arr in arrs]
    M = np.array([item.max()+1 for item in indices])
    M = np.r_[1, M[:-1]]
    strides = M.cumprod()
    indices = np.stack(indices, axis=-1)
    vals = (indices * strides).sum(axis=-1)
    uniqs, labels = np.unique(vals, return_inverse=True)
    labels = labels.reshape(arrs[0].shape)
    return labels, uniqs

 

Now stop reading here if you aren't interested in the details, but have a quasi-firm grasp on the process.  Read on if you must, but suffice to say unubtu was on to something.

 

-----------------------------  details .... details .... details... read at your peril  ------------------------------------------------------------

In essence, what the code does, is find the unique values of each array shoved into the function.  In our case, the list of unique values is short (see line 3 below).  Not too onerous a combination slate.

Line 10 can be modified to show the unique values in each array... useful if they aren't sequential.  This is show below in line 1.

induni = [np.unique(arr, return_inverse=True)[0] for arr in arrs]  # arrays a, b, c

[array([0, 1, 2, 4]), array([0, 1, 2, 3, 4, 5]), array([0, 1, 2, 3, 4, 5])]

 

The nifty bit, is how unutbu constructs the 'indices' in two steps, using the maximum value in each array, determines the unique combination, applies it to the arrays to get a unique value ('val') for each combination and then reassigns a class ('label' to it).

 

This process, not documented in the code, is expanded upon here.  In lines 1-4 below, an array, combos, is constructed from the values in the array plus the 'val' derived from the above process.

idx = [np.unique(arr, return_inverse=True)[1] for arr in arrs]
vals = (indices * strides).sum(axis=-1)
dt = [('a', '<i8'), ('b', '<i8'), ('c', '<i8'), ('vals', '<i8')]
combos = np.array(list(zip(*idx)), dtype= dt)

 

If the 'combos' array is sorted, the unique values can be extracted from it.  The values in the 'vals' column is used to assign a 'label' to each unique combination.

srt_order = np.argsort(combos, order=['a','b', 'c'])  # sort the combos

srted = combos[srt_order]  # ---- extract them in sorted order

np.unique(srted)  # ---- determine the unique combinations

array([(0, 0, 0,   0),
       (0, 2, 2,  56),
       (1, 0, 2,  49),
       (1, 2, 0,   9),
       (1, 4, 4, 113),
       (2, 3, 1,  38),
       (2, 5, 1,  46),
       (3, 1, 3,  79),
       (3, 1, 5, 127)],
      dtype=[('a', '<i8'), ('b', '<i8'),
             ('c', '<i8'), ('vals', '<i8')])

 

Hope it is clear as mud, and I have some experimenting to do, but I just thought I would throw it out there.

Yup... adding points to poly* features

Already offered... Densify at all license levels?  NO!

My missive .... Densification by Factor  YES!

 

 

There is also a Geodesic Densify, but I don't do Geographic in any event, so use esri tools

Why did I reinvent the wheel?  Because the densification basis is to simply provide a sequence of locations, for X and Y then subdivide the intervening location by some factor (hence the ... by Factor)

 

Here is an example....

a = np.array([[ 20.,  20.],
              [ 20.,  30.],
              [ 30.,  30.],
              [ 30.,  20.],
              [ 20.,  20.]])

 

Then you run some code

def _densify_2D(a, fact=2):
    """Densify a 2D array using np.interp.
    """

    # Y = a changed all the y's to a
    a = np.squeeze(a)
    n_fact = len(a) * fact
    b = np.arange(0, n_fact, fact)
    b_new = np.arange(n_fact - 1)
    c0 = np.interp(b_new, b, a[:, 0])
    c1 = np.interp(b_new, b, a[:, 1])
    n = c0.shape[0]
    c = np.zeros((n, 2))
    c[:, 0] = c0
    c[:, 1] = c1
    return c

And you get this.

a_dens = _densify_2D(a, fact=2)

a_dens

array([[ 20.,  20.],
       [ 20.,  25.],
       [ 20.,  30.],
       [ 25.,  30.],
       [ 30.,  30.],
       [ 30.,  25.],
       [ 30.,  20.],
       [ 25.,  20.],
       [ 20.,  20.]])
a.shape  # (49874, 2)

%timeit _densify_2D(a, fact=2)

5.25 ms ± 323 µs per loop
     (mean ± std. dev. of 7 runs,
      100 loops each)

 

 

Where it really shines, is for long features that just need to be densified without any particular requirement for a set spacing.  If is fairly speedy as well.

Hope someone finds this useful, if you don't have the Standard or Advanced license in particular.

This is the current state of the Point and PointGeometry .... not the best....

pt2 = arcpy.Point()

pt2
<Point (0.0, 0.0, #, #)>        # ---- make a point

pg2 = arcpy.PointGeometry(pt2)  # ---- make a point geometry

pg2.WKT                         # ---- why does it have valid X, Y values????
'POINT (0 0)'

So there is no such thing as a 'null' point (but a 'nominal' point as one of our esteemed participants noted). 

'None' isn't good enough.  You can create an 'empty' geometry recognized by other representations with a little trickery.

import arcpy

n = np.nan

a = np.array([n, n, n, n])

pt = arcpy.Point()

pt.X, pt.Y, pt.Z, pt.M = a

pg = arcpy.PointGeometry(pt)

pg
<PointGeometry object at 0x200c6ed6ac8[0x200bf9a39e0]>

pg.WKT
'POINT EMPTY'

pt
<Point (nan, nan, nan, nan)>

 

And when you look at the PointGeometry, you decide which is best\.

# ---- by using NaN (not a number) for point properties ----
pt        # <Point (nan, nan, nan, nan)>

pg = arcpy.PointGeometry(pt)

pg.JSON   # '{"x":"NaN","y":"NaN","spatialReference":{"wkid":null}}'

pg.WKT    #  'POINT EMPTY'

# ---- The current state of affairs ----

pt2       # <Point (0.0, 0.0, #, #)>

pg2 = arcpy.PointGeometry(pt2)

pg2.JSON  # '{"x":0,"y":0,"spatialReference":{"wkid":null}}'

pg2.WKT   # 'POINT (0 0)'

 

I should point out that you can use math.nan inplace of np.nan

m = math.nan
b = np.array([m, m, m, m])
pt3 = arcpy.Point()
pt3.X, pt3.Y, pt3.Z, pt3.M = b
pt3  # ---- <Point (nan, nan, nan, nan)>

pg3 = arcpy.PointGeometry(pt3)

pg3.JSON  # ---- '{"x":"NaN","y":"NaN","spatialReference":{"wkid":null}}'

pg3.WKT   # ---- 'POINT EMPTY'

Natural terrain is so fickle.  Wouldn't it be nice to have a sink-less terrain so that your flowdirection, flowaccumulation and stream generation would go as nature planned.

 

If you have a need to test algorithms, then natural surfaces aren't the best, but generating artificial terrain can be a pain.  They tend to have repetitive patterns and lack the finesse that real surface has.  

 

This is a brief look at one of the many algorithms that have been employed in early gaming software and largely ignored by the 'derivative' people.  The algorithm is pretty simple.... provide random seeds to the corners of a square array/raster, then divide it into triangles and squares in a predictable fashion calculating the averages of the intervening space, assigning the central location, the value of the pattern points.  Now, to mix it up, one can add some randomness to the calculated values to 'roughen' the surface to so speak.

 

The following is the sequence of the pattern generation... the code is on my GitHub account and the Code Sharing site for those that want to explore on their own.  I repurposed an implementation adding a few twists, like leaving out the random jitters and the ability to replicate in both rows and columns should you want to generate a surface with some repetition.

 

More importantly, this is for fun and is a good way to learn how to use numpy arrays to work with array data that forms images or rasters in the GIS sense of the term.

 

Input template array with 4 corners defined...  The empty space is represented
by 'nan' aka, 'not a number', you could use zero, but 'not a number' is a
raster null, without having to get into the issues of masks and the like.

Step 1... throw in 4 corners

[[ 0.250  nan  nan  nan  nan  nan  nan  nan  0.750]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ 0.500  nan  nan  nan  nan  nan  nan  nan  1.000]]

Step 2... calculate the central location...

[[ 0.250  nan  nan  nan  nan  nan  nan  nan  0.750]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ nan  nan  nan  nan  0.625  nan  nan  nan  nan]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ 0.500  nan  nan  nan  nan  nan  nan  nan  1.000]]

Step 3...
Split into the diamond structure and calculate the compass points

[[ 0.250  nan  nan  nan  0.542  nan  nan  nan  0.750]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ 0.458  nan  nan  nan  0.625  nan  nan  nan  0.792]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ 0.500  nan  nan  nan  0.708  nan  nan  nan  1.000]]

Step 4....
Calculate the central values for the squares

[[ 0.250  nan  nan  nan  0.542  nan  nan  nan  0.750]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ nan  nan  0.469  nan  nan  nan  0.677  nan  nan]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ 0.458  nan  nan  nan  0.625  nan  nan  nan  0.792]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ nan  nan  0.573  nan  nan  nan  0.781  nan  nan]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ 0.500  nan  nan  nan  0.708  nan  nan  nan  1.000]]

Step ....
Repeat......
[[ 0.250  nan  0.420  nan  0.542  nan  0.656  nan  0.750]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ 0.392  nan  0.469  nan  0.578  nan  0.677  nan  0.740]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ 0.458  nan  0.531  nan  0.625  nan  0.719  nan  0.792]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ 0.510  nan  0.573  nan  0.672  nan  0.781  nan  0.858]
[ nan  nan  nan  nan  nan  nan  nan  nan  nan]
[ 0.500  nan  0.594  nan  0.708  nan  0.830  nan  1.000]]

Repeat
[[ 0.250  nan  0.420  nan  0.542  nan  0.656  nan  0.750]
[ nan  0.383  nan  0.502  nan  0.613  nan  0.706  nan]
[ 0.392  nan  0.469  nan  0.578  nan  0.677  nan  0.740]
[ nan  0.463  nan  0.551  nan  0.650  nan  0.732  nan]
[ 0.458  nan  0.531  nan  0.625  nan  0.719  nan  0.792]
[ nan  0.518  nan  0.600  nan  0.699  nan  0.787  nan]
[ 0.510  nan  0.573  nan  0.672  nan  0.781  nan  0.858]
[ nan  0.544  nan  0.637  nan  0.748  nan  0.867  nan]
[ 0.500  nan  0.594  nan  0.708  nan  0.830  nan  1.000]]


Step.... finally... all the squares and triangles are filled in

[[ 0.250  0.351  0.420  0.488  0.542  0.604  0.656  0.704  0.750]
[ 0.342  0.383  0.443  0.502  0.559  0.613  0.663  0.706  0.732]
[ 0.392  0.427  0.469  0.525  0.578  0.630  0.677  0.714  0.740]
[ 0.438  0.463  0.503  0.551  0.601  0.650  0.694  0.732  0.754]
[ 0.458  0.493  0.531  0.577  0.625  0.673  0.719  0.757  0.792]
[ 0.496  0.518  0.556  0.600  0.649  0.699  0.747  0.787  0.812]
[ 0.510  0.536  0.573  0.620  0.672  0.725  0.781  0.823  0.858]
[ 0.518  0.544  0.587  0.637  0.691  0.748  0.807  0.867  0.908]
[ 0.500  0.546  0.594  0.646  0.708  0.762  0.830  0.899  1.000]]

 

Of course, the image that follows uses a much larger array (1025x1025).  The generated array, was sink-less, and flow properties and a stream network was generated from the result.  Pretty smooth eh?

 

Next up... terrain generation from polygons and other vector geometry

Filter

Started this with rolling stats. 

 

Rolling Statistics for grid data... an alternative

 

Once you have a filter... kernel whatever you want to call it, all you need is a moving window of some form and shape.  I won't go into those details since my previous work covers this.  I will just report a few simple ones, and show you some of the possible kernels that you can create.

 

# ---- my input array ----

Input array.......
-shape (1, 10, 10), ndim 3
.... 1 1 1 1 1 2 2 2 2 2 
.... 1 1 1 1 1 2 2 2 2 2 
.... 1 1 1 1 1 2 2 2 2 2 
.... 1 1 1 1 1 2 2 2 2 2 
.... 1 1 1 1 1 2 2 2 2 2 
.... 3 3 3 3 3 4 4 4 4 4 
.... 3 3 3 3 3 4 4 4 4 4 
.... 3 3 3 3 3 4 4 4 4 4 
.... 3 3 3 3 3 4 4 4 4 4 
.... 3 3 3 3 3 4 4 4 4 4

From the above, you can see that I kept it simple... 4 blocks of values, the perfect study area.

 

Let's see if we can pick up some lines that divide the raster/array

Sobel horizonal filter...
-shape (1, 3, 3), ndim 3

  . -1 -2 -1 
  .  0  0  0 
  .  1  2  1 


Sobel horizontal result.......
-shape (1, 10, 10), ndim 3, masked array, fill value 999999

.... -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 
.... -1  0  0  0  0  0  0  0  0 -1 
.... -1  0  0  0  0  0  0  0  0 -1 
.... -1  0  0  0  0  0  0  0  0 -1 
.... -1  8  8  8  8  8  8  8  8 -1 
.... -1  8  8  8  8  8  8  8  8 -1 
.... -1  0  0  0  0  0  0  0  0 -1 
.... -1  0  0  0  0  0  0  0  0 -1 
.... -1  0  0  0  0  0  0  0  0 -1 
.... -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 

See those 8's! nailed the horizontal divide.

 

vertical filter...

  . -1  0  1 
  . -2  0  2 
  . -1  0  1 

Sobel vertical result.......
-shape (1, 10, 10), ndim 3, masked array, fill value 999999
.... -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 
.... -1  0  0  0  4  4  0  0  0 -1 
.... -1  0  0  0  4  4  0  0  0 -1 
.... -1  0  0  0  4  4  0  0  0 -1 
.... -1  0  0  0  4  4  0  0  0 -1 
.... -1  0  0  0  4  4  0  0  0 -1 
.... -1  0  0  0  4  4  0  0  0 -1 
.... -1  0  0  0  4  4  0  0  0 -1 
.... -1  0  0  0  4  4  0  0  0 -1 
.... -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 

 

So easy... vertical nailed as well.

 

One more...

Array...
-shape (1, 3, 3), ndim 3
  . -1 -1  0 
  . -1  0  1 
  .  0  1  1 

Emboss vertical result.......

-shape (1, 10, 10), ndim 3, masked array, fill value 999999
.... -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 
.... -1  0  0  0  2  2  0  0  0 -1 
.... -1  0  0  0  2  2  0  0  0 -1 
.... -1  0  0  0  2  2  0  0  0 -1 
.... -1  4  4  4  6  6  4  4  4 -1 
.... -1  4  4  4  6  6  4  4  4 -1 
.... -1  0  0  0  2  2  0  0  0 -1 
.... -1  0  0  0  2  2  0  0  0 -1 
.... -1  0  0  0  2  2  0  0  0 -1 
.... -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 

 

How about some filters that I have found.

    n = np.nan
    all_f = [1, 1, 1, 1, 1, 1, 1, 1, 1]
    no_cnt = [1, 1, 1, 1, n, 1, 1, 1, 1]
    cross_f = [1, 0, 1, 0, 0, 0, 1, 0, 1]
    plus_f = [0, 1, 0, 1, 0, 1, 0, 1, 0]
    grad_e = [1, 0, -1, 2, 0, -2, 1, 0, -1]
    grad_n = [-1, -2, -1, 0, 0, 0, 1, 2, 1]
    grad_ne = [0, -1, -2, 1, 0, -1, 2, 1, 0]
    grad_nw = [2, -1, 0, -1, 0, 1, 0, 1, 2]
    grad_s = [1, 2, 1, 0, 0, 0, -1, -2, -1]
    grad_w = [-1, 0, 1, -2, 0, 2, -1, 0, 1]
    lap_33 = [0, -1, 0, -1, 4, -1, 0, -1, 0]
    line_h = [-1, -1, -1, 2, 2, 2, -1, -1, -1]
    line_ld = [2, -1, -1, -1, 2, -1, -1, -1, 2]
    line_rd = [-1, -1, 2, -1, 2, -1, 2, -1, -1]
    line_v = [-1, 0, -1, -1, 2, -1, -1, 2, -1]
    high = [-0.7, -1.0, -0.7, -1.0, 6.8, -1.0, -0.7, -1.0, -0.7]
    sob_hor = [-1, -2, -1, 0, 0, 0, 1, 2, 1]   # sobel y
    sob_vert = [-1, 0, 1, -2, 0, 2, -1, 0, 1]  # sobel x
    emboss = [-1, -1, 0, -1, 0, 1, 0, 1, 1]
    sharp1 = [0., -0.25, 0., -0.25, 2.0, -0.25, 0., -0.25, 0.]
    sharp2 = [-0.25, -0.25, -0.25, -0.25, 3.0, -0.25, -0.25, -0.25, -0.25]
    sharp3 = [-1, -1, -1, -1, 9, -1, -1, -1, -1]
    lowpass = [1, 2, 1, 2, 4, 2, 1, 2, 1]

    # ---- assemble the dictionary ----

    d = {1: all_f, 2: no_cnt, 3: cross_f, 4: plus_f, 5: grad_e,
         6: grad_n, 7: grad_ne, 8: grad_nw, 9: grad_s, 10: grad_w,
         11: lap_33, 12: line_h, 13: line_ld, 14: line_rd, 15: line_v,
         16: high, 17: sob_hor, 18: sob_vert, 19: emboss, 20: sharp1,
         21: sharp2, 23: sharp3, 24: lowpass}

    filter_ = np.array(d[mode]).reshape(3, 3)

You call the filter using the mode number (ie d[23] is sharp3.

import numpy as np
from numpy.lib.stride_tricks import as_strided

def stride(a, win=(3, 3), stepby=(1, 1)):
    """Provide a 2D sliding/moving view of an array.
    :  There is no edge correction for outputs. Use _pad first.
    """

    err = """Array shape, window and/or step size error.
    Use win=(3,3) with stepby=(1,1) for 2D array
    or win=(1,3,3) with stepby=(1,1,1) for 3D
    ----    a.ndim != len(win) != len(stepby) ----
    """

    assert (a.ndim == len(win)) and (len(win) == len(stepby)), err
    shape = np.array(a.shape)  # array shape (r, c) or (d, r, c)
    win_shp = np.array(win)    # window      (3, 3) or (1, 3, 3)
    ss = np.array(stepby)      # step by     (1, 1) or (1, 1, 1)
    newshape = tuple(((shape - win_shp) // ss) + 1)
    newshape += tuple(win_shp)
    newstrides = tuple(np.array(a.strides) * ss) + a.strides
    a_s = as_strided(a, shape=newshape, strides=newstrides).squeeze()
    return a_s


def a_filter(a, mode=1, ignore_nodata=True, nodata=None):
    """put the filter stuff below this line"""
    #
    # ---- stride the input array ----
    a_strided = stride(a)
    if ignore_nodata:
        c = (np.sum(a_strided * filter_, axis=(2, 3)))
    else:
        c = (np.nansum(a_strided * filter_, axis=(2, 3)))
    pad_ = nodata
    if nodata is None:
        if c.dtype.name in ('int', 'int32', 'int64'):
            pad_ = min([0, -1, a.min()-1])
        else:
            pad_ = min([0.0, -1.0, a.min()-1])
    c = np.lib.pad(c, (1, 1), "constant", constant_values=(pad_, pad_))
    m = np.where(c == pad_, 1, 0)
    c = np.ma.array(c, mask=m, fill_value=None)
    return filter_, c

# ----------------------------------------------------------------------
# __main__ .... code section
if __name__ == "__main__":
    """Optionally...
    : - print the script source name.
    : - run the _demo
    """

#    you can put your stuff here... use RasterToNumPyArray to get an array

 

So that is pretty well some other stuff.

I should point out that the same process can be used for statistical calculations, terrain derivatives (slope, aspect, hillshade, curvature etc ) and I have reported on these before.

 

Now I should point out that none of the above requires the Spatial Analyst extension and the RasterToNumPyArray and NumPyArrayToRaster are provided with arcpy.da 's module.  Lots of possibilities, lots of capability is built-in.

 

More later.

Rolling stats...  The 3x3 moving/sliding/rolling window.. This kernel type will focus on returning statistics like your FocalStatistics functions... but without needing the Spatial Analyst extension...

 

They are used ....

And there are other uses....

 

First the results.

a      # ---- A grid/raster/array ----

array([[5, 1, 2, 0, 5, 4, 7, 3, 2, 2],
       [2, 5, 9, 9, 6, 9, 2, 1, 3, 0],
       [7, 2, 7, 4, 4, 4, 3, 7, 7, 7],
       [3, 9, 6, 8, 0, 1, 8, 5, 0, 7],
       [8, 4, 6, 8, 1, 9, 6, 9, 5, 5],
       [3, 9, 9, 5, 3, 1, 5, 5, 2, 3],
       [0, 5, 2, 0, 9, 4, 7, 9, 9, 0],
       [2, 5, 5, 4, 3, 7, 7, 8, 3, 5],
       [5, 2, 4, 5, 3, 6, 5, 2, 7, 2],
       [4, 6, 4, 9, 8, 6, 8, 6, 2, 0]])


a_s = stride(a, r_c=(3,3))  # ---- create an array of rolling windows 3x3 ----

# ---- Now calculate all the statistics --------------------------------------

rolling_stats(a_s, no_null=False, prn=True)
Min...
[[1 0 0 0 2 1 1 0]
[2 2 0 0 0 1 0 0]
[2 2 0 0 0 1 0 0]
[3 4 0 0 0 1 0 0]
[0 0 0 0 1 1 2 0]
[0 0 0 0 1 1 2 0]
[0 0 0 0 3 2 2 0]
[2 2 3 3 3 2 2 0]]
Max...
[[9 9 9 9 9 9 7 7]
[9 9 9 9 9 9 8 7]
[9 9 8 9 9 9 9 9]
[9 9 9 9 9 9 9 9]
[9 9 9 9 9 9 9 9]
[9 9 9 9 9 9 9 9]
[5 5 9 9 9 9 9 9]
[6 9 9 9 8 8 8 8]]
Mean...
[[ 4.44  4.33  5.11  5.    4.89  4.44  3.89  3.56]
[ 5.56  6.56  5.89  5.    4.11  4.44  4.    4.11]
[ 5.78  6.    4.89  4.33  4.    5.78  5.56  5.78]
[ 6.33  7.11  5.11  4.    3.78  5.44  5.    4.56]
[ 5.11  5.33  4.78  4.44  5.    6.11  6.33  5.22]
[ 4.44  4.89  4.44  4.    5.11  5.89  6.11  4.89]
[ 3.33  3.56  3.89  4.56  5.67  6.11  6.33  5.  ]
[ 4.11  4.89  5.    5.67  5.89  6.11  5.33  3.89]]
Med...
[[ 5.  4.  5.  4.  4.  4.  3.  3.]
[ 6.  7.  6.  4.  4.  4.  3.  5.]
[ 6.  6.  6.  4.  4.  6.  6.  7.]
[ 6.  8.  6.  3.  3.  5.  5.  5.]
[ 5.  5.  5.  4.  5.  6.  6.  5.]
[ 5.  5.  4.  4.  5.  7.  7.  5.]
[ 4.  4.  4.  4.  6.  7.  7.  5.]
[ 4.  5.  4.  6.  6.  6.  6.  3.]]
Sum...
[[40 39 46 45 44 40 35 32]
[50 59 53 45 37 40 36 37]
[52 54 44 39 36 52 50 52]
[57 64 46 36 34 49 45 41]
[46 48 43 40 45 55 57 47]
[40 44 40 36 46 53 55 44]
[30 32 35 41 51 55 57 45]
[37 44 45 51 53 55 48 35]]
Std...
[[ 2.67  3.2   2.85  2.62  2.02  2.5   2.28  2.59]
[ 2.59  2.36  2.73  3.09  2.88  2.83  2.71  2.96]
[ 2.2   2.16  2.73  3.16  2.98  2.62  2.59  2.39]
[ 2.4   1.79  3.    3.37  3.15  2.83  2.58  2.5 ]
[ 3.    2.91  3.26  3.34  2.87  2.56  2.26  3.08]
[ 2.91  2.73  2.83  2.62  2.42  2.28  2.38  3.03]
[ 1.76  1.71  2.33  2.45  1.94  2.02  2.36  3.2 ]
[ 1.29  1.79  2.    2.    1.79  1.73  2.31  2.56]]
Var...
[[  7.14  10.22   8.1    6.89   4.1    6.25   5.21   6.69]
[  6.69   5.58   7.43   9.56   8.32   8.02   7.33   8.77]
[  4.84   4.67   7.43  10.     8.89   6.84   6.69   5.73]
[  5.78   3.21   8.99  11.33   9.95   8.02   6.67   6.25]
[  8.99   8.44  10.62  11.14   8.22   6.54   5.11   9.51]
[  8.47   7.43   8.02   6.89   5.88   5.21   5.65   9.21]
[  3.11   2.91   5.43   6.02   3.78   4.1    5.56  10.22]
[  1.65   3.21   4.     4.     3.21   2.99   5.33   6.54]]
Range

# ---- If you have nodata values, they can be accommodated ------------------

 

Now the code functions

 

import numpy as np
from numpy.lib.stride_tricks import as_strided

def _check(a, r_c, subok=False):
    """Performs the array checks necessary for stride and block.
    : a   - Array or list.
    : r_c - tuple/list/array of rows x cols.
    : subok - from numpy 1.12 added, keep for now
    :Returns:
    :------
    :Attempts will be made to ...
    :  produce a shape at least (1*c).  For a scalar, the
    :  minimum shape will be (1*r) for 1D array or (1*c) for 2D
    :  array if r<c.  Be aware
    """

    if isinstance(r_c, (int, float)):
        r_c = (1, int(r_c))
    r, c = r_c
    if a.ndim == 1:
        a = np.atleast_2d(a)
    r, c = r_c = (min(r, a.shape[0]), min(c, a.shape[1]))
    a = np.array(a, copy=False, subok=subok)
    return a, r, c, tuple(r_c)


def _pad(a, nan_edge=False):
    """Pad a sliding array to allow for stats"""
    if nan_edge:
        a = np.pad(a, pad_width=(1, 2), mode="constant",
                   constant_values=(np.NaN, np.NaN))
    else:
        a = np.pad(a, pad_width=(1, 1), mode="reflect")
    return a


def stride(a, r_c=(3, 3)):
    """Provide a 2D sliding/moving view of an array.
    :  There is no edge correction for outputs.
    :
    :Requires:
    :--------
    : _check(a, r_c) ... Runs the checks on the inputs.
    : a - array or list, usually a 2D array.  Assumes rows is >=1,
    :     it is corrected as is the number of columns.
    : r_c - tuple/list/array of rows x cols.  Attempts  to
    :     produce a shape at least (1*c).  For a scalar, the
    :     minimum shape will be (1*r) for 1D array or 2D
    :     array if r<c.  Be aware
    """

    a, r, c, r_c = _check(a, r_c)
    shape = (a.shape[0] - r + 1, a.shape[1] - c + 1) + r_c
    strides = a.strides * 2
    a_s = (as_strided(a, shape=shape, strides=strides)).squeeze()
    return a_s


# ----------------------------------------------------------------------
# (18) rolling stats .... code section
def rolling_stats(a, no_null=True, prn=True):
    """Statistics on the last two dimensions of an array.
    :Requires
    :--------
    : a - 2D array  **Note, use 'stride' above to obtain rolling stats
    : no_null - boolean, whether to use masked values (nan) or not.
    : prn - boolean, to print the results or return the values.
    :
    :Returns
    :-------
    : The results return an array of 4 dimensions representing the original
    : array size and block size
    : eg.  original = 6x6 array   block = 3x3
    :      breaking the array into 4 chunks
    """

    a = np.asarray(a)
    a = np.atleast_2d(a)
    ax = None
    if a.ndim > 1:
        ax = tuple(np.arange(len(a.shape))[-2:])
    if no_null:
        a_min = a.min(axis=ax)
        a_max = a.max(axis=ax)
        a_mean = a.mean(axis=ax)
        a_med = np.median(a, axis=ax)
        a_sum = a.sum(axis=ax)
        a_std = a.std(axis=ax)
        a_var = a.var(axis=ax)
        a_ptp = a_max - a_min
    else:
        a_min = np.nanmin(a, axis=(ax))
        a_max = np.nanmax(a, axis=(ax))
        a_mean = np.nanmean(a, axis=(ax))
        a_med = np.nanmedian(a, axis=(ax))
        a_sum = np.nansum(a, axis=(ax))
        a_std = np.nanstd(a, axis=(ax))
        a_var = np.nanvar(a, axis=(ax))
        a_ptp = a_max - a_min
    if prn:
        s = ['Min', 'Max', 'Mean', 'Med', 'Sum', 'Std', 'Var', 'Range']
        frmt = "...\n{}\n".join([i for i in s])
        v = [a_min, a_max, a_mean, a_med, a_sum, a_std, a_var, a_ptp]
        args = [indent(str(i), '... ') for i in v]
        print(frmt.format(*args))
    else:
        return a_min, a_max, a_mean, a_med, a_sum, a_std, a_var, a_ptp

Curses at cursors... pretty well a regular occurrence on this site.  People love them and love to hate them.

They try to nest them, and nest them within 'for' loops and 'with' statements with calls that sound poetic ( ... for row in rows, with while for arc thou ... )

 

There are old cursors and new cursors (the da-less and the da cursors).  Cursors appear in other guises such as the new, and cleverly named, 'arcgis' module (digression # 1 ... really? something else with arcgis in it! Who is in charge of branding).

 

Perhaps cursors are cloaked, in other arcpy and  data access module methods (ie. blank-to-NumPyArray and NumPyArray-to-blank).  Who knows for sure since much is locked in arcgisscripting.pyd.

 

Sadly, we deal in a work of mixed data types.  Our tables contain columns of attributes, organized sequentially by rows.  Sometimes the row order has meaning, sometimes not.   Each column contains one data type in a well-formed data structure.  This is why spreadsheets are purely evil for trying to create and maintain data structure, order and form (you can put anything anywhere).

 

Harkening back to the fields of mathematics, arrays are assemblages of data, in 1, 2, 3 or more dimensions.  If an array of any dimension has a uniform data type, then life is easier from a structural and usage perspective (this is one reason why Remote Sensing is easier than GIS ... bring on the mail).  We need to maintain an index which ties our geometry to our attributes so what goes where, and where is what, doesn't get mixed up (digression # 2... I am sure this isn't what the branders meant by The Science of Where but we can only hope)

 

Nerd Stuff 

Enough with the boring stuff... bring on the code.

Some of this has been looked from a slightly different perspective in 

 

Get to the Points... arcpy, numpy, pandas

We need some data to work with so... how about a square.

    in_fc = r"C:\Your_spaceless_path\Your.gdb\square"

The 'Describe' object

 

The 'describe' object does just that: describes an object, in this case a FeatureClass.
desc = arcpy.da.Describe(in_fc)
In the newer arcpy.da module, the values can be accessed from a dictionary.  A quick way to get the sorted dictionary keys is to use a list comprehension.  If you want the values, then you can obtain them in a similar fashion.
sk = sorted([k for k in desc.keys()])  # sorted keys
kv = [(k, desc[k]) for k in sk]  # key/value pairs
kv = "\n".join(["{!s:<20} {}".format(k, desc[k]) for k in sk])
Some useful keys associated with featureclasses are extracted as follows:
With appropriate snips in the full list

[..., 'MExtent', 'OIDFieldName', 'ZExtent', 'aliasName', 'areaFieldName',
  'baseName', ... 'catalogPath', ... 'dataElementType', 'dataType',
  'datasetType', ... 'extent', 'featureType', 'fields', 'file', ... 'hasM',
  'hasOID', 'hasSpatialIndex', 'hasZ', 'indexes', ... 'lengthFieldName',
  ... 'name', 'path', ... 'shapeFieldName', 'shapeType', 'spatialReference',
  ...]

The Cursor object
A cursor gives access to the geometry and attributes of a featureclass.  It is always recommended to use the spatial reference when creating the cursor.  Historically (fixed?), if it was omitted, geometry discrepancies would arise. This information can easily be obtained from the describe object from the previous section.
SR = desc['spatialReference']  # Get the search cursor object.
flds = "*"
args = [in_fc, flds, None, SR, True, (None, None)]
cur = arcpy.da.SearchCursor(*args)
See what it reveals...
dir(cur)

['__class__', '__delattr__', '__dir__', '__doc__', '__enter__', '__eq__',
'__esri_toolinfo__', '__exit__', '__format__', '__ge__', '__getattribute__',
'__getitem__', '__gt__', '__hash__', '__init__', '__iter__', '__le__',
'__lt__', '__ne__', '__new__', '__next__', '__reduce__', '__reduce_ex__',
'__repr__', '__setattr__', '__sizeof__', '__str__', '__subclasshook__',
'_as_narray', '_dtype', 'fields', 'next', 'reset']
 
Individual properties are:
cur.__class__
<class 'da.SearchCursor'>

cur.__class__.__mro__
(<class 'da.SearchCursor'>, <class 'object'>)
The search cursor inherits from 'object' and the difference in their properties and methods offered by the da.SearchCursor can be determined as follows:
s0 = set(dir(mro[0]))
s1 = set(dir(mro[1]))

sorted(list(set.difference(s0, s1)))

 ['__enter__', '__esri_toolinfo__', '__exit__', '__getitem__', '__iter__',
  '__next__', '_as_narray', '_dtype', 'fields', 'next', 'reset']
 
The cursor offers the means to process its objects.
cur.__esri_toolinfo__
 ['FeatureLayer|Table|TableView|Dataset|FeatureDataset::::', 'String::*::',
  'Python::None::', 'CoordinateSystem::::']
 
This handy one returns a numpy structured/recarray
type(cur._as_narray())
 <class 'numpy.ndarray'>
How about information on the attributes of in_fc (more about this later)
cur._as_narray().__array_interface__

 {'version': 3,
  'strides': None,
  'shape': (0,),
  'typestr': '|V36',
  'descr': [('OBJECTID', '<i4'),
            ('Shape', '<f8', (2,)),
            ('Shape_Length', '<f8'),
            ('Shape_Area', '<f8')],
  'data': (2044504703824, False)}
 
We got a glimpse of the field names and data types from the __array_interface__, but this information can be accessed directly as well.
cur.fields
 ('OBJECTID', 'Shape', 'Shape_Length', 'Shape_Area')

cur._dtype
dtype([('OBJECTID', '<i4'),
       ('Shape', '<f8', (2,)),
       ('Shape_Length', '<f8'),
       ('Shape_Area', '<f8')])    
Now, the gotcha's.  We created our search cursor at the beginning and each record was cycled through until it reached the end.  If we attempt to get its properties we may in for a surprise, so we need to 'reset' the cursor back to the start.
cur._as_narray()  # try to get its properties, all we get is the dtype
array([],
    dtype=[('OBJECTID', '<i4'), ('Shape', '<f8', (2,)),
     ('Shape_Length', '<f8'), ('Shape_Area', '<f8')])
Once the cursor is reset, the array values for the square are revealed with the appropriate data type.
cur.reset()
cur._as_narray()  # reset to the beginning  

array([(1, [342000.0, 5022000.0], 4000.0, 1000000.0),
       (1, [342000.0, 5023000.0], 4000.0, 1000000.0),
       (1, [343000.0, 5023000.0], 4000.0, 1000000.0),
       (1, [343000.0, 5022000.0], 4000.0, 1000000.0),
       (1, [342000.0, 5022000.0], 4000.0, 1000000.0)],
    dtype=[('OBJECTID', '<i4'), ('Shape', '<f8', (2,)),
           ('Shape_Length', '<f8'), ('Shape_Area', '<f8')])
    
There is no automatic reset, so be careful.  You can print the objects in the array in a couple of ways.
cur.reset()
for row in cur:
   print(("{} "*len(row)).format(*row))  # print individual elements

  1 (342000.0, 5022000.0) 4000.0 1000000.0
 1 (342000.0, 5023000.0) 4000.0 1000000.0
 1 (343000.0, 5023000.0) 4000.0 1000000.0
 1 (343000.0, 5022000.0) 4000.0 1000000.0
 1 (342000.0, 5022000.0) 4000.0 1000000.0
Resetting the cursor, and print again.
cur.reset()
for row in cur:
    print(row)  # print the whole row as a tuple

   (1, (342000.0, 5022000.0), 4000.0, 1000000.0)
   (1, (342000.0, 5023000.0), 4000.0, 1000000.0)
   (1, (343000.0, 5023000.0), 4000.0, 1000000.0)
   (1, (343000.0, 5022000.0), 4000.0, 1000000.0)
   (1, (342000.0, 5022000.0), 4000.0, 1000000.0)
Of course since generator-like objects can be converted to a list, that can be done as an alternative, particularly if you have the memory and wish to deal with list objects instead.
cur.reset()
list(cur)
 [(1, (342000.0, 5022000.0), 4000.0, 1000000.0),
  (1, (342000.0, 5023000.0), 4000.0, 1000000.0),
  (1, (343000.0, 5023000.0), 4000.0, 1000000.0),
  (1, (343000.0, 5022000.0), 4000.0, 1000000.0),
  (1, (342000.0, 5022000.0), 4000.0, 1000000.0)]
 
So if you know the data type of the components of the cursor, you can go to the ndarray in an indirect fashion.
cur.reset()
dt = cur._dtype
c_lst = list(cur)

np.asarray(c_lst, dtype=dt)

array([(1, [342000.0, 5022000.0], 4000.0, 1000000.0),
       (1, [342000.0, 5023000.0], 4000.0, 1000000.0),
       (1, [343000.0, 5023000.0], 4000.0, 1000000.0),
       (1, [343000.0, 5022000.0], 4000.0, 1000000.0),
       (1, [342000.0, 5022000.0], 4000.0, 1000000.0)],
    dtype=[('OBJECTID', '<i4'), ('Shape', '<f8', (2,)),
           ('Shape_Length', '<f8'), ('Shape_Area', '<f8')])
The ndarray can be viewed as a record array.  Since the data type and structure remain the same, a 'view' of the array as a record array (recarray).  Record arrays allow the user to slice the array using conventional array slicing or by object dot notation.  
a = a.view(np.recarray)

a.Shape == a['Shape']  # check to see if slicing equals dot notation

array([[ True,  True],
      [ True,  True],
       [ True,  True],
       [ True,  True],
       [ True,  True]], dtype=bool)
Or more simply...
np.all(a.Shape == a['Shape'])
True       

a.Shape  # or a['Shape']

array([[  342000.,  5022000.],
       [  342000.,  5023000.],
       [  343000.,  5023000.],
       [  343000.,  5022000.],
       [  342000.,  5022000.]])
 
You can calculate the properties of the objects simply, but in the case of polygons, the duplicate start/end point should be reduced to a singleton.  In the examples, the object's shape is obtained, then the desired property is derived on a column basis.
pnts = a.Shape[:-1]       # get the unique points

cent = pnts.mean(axis=0)  # return the mean by column

cent array([  342500.,  5022500.])
 
With some fancy work, and calling one of my previously defined array functions in the 'arraytools' module, you can do things like determine interpoint distances.
import arraytools as art

art.e_dist(cent, pnts)

array([ 707.11,  707.11,  707.11,  707.11])
Which is correct given the square polygon's shape.
Another example, to demonstrate array functions.  In the case of polygons, it is useful to have the first and last point (ie duplicates) retained to ensure closure of the polygon. 
poly = a.Shape

art.e_leng(poly)  # method to return polygon perimeter/length, total, then by segment

 (4000.0, [array([[ 1000.,  1000.,  1000.,  1000.]])])

art.e_area(poly)

 1000000.0
 
------------------------------------------------------------------------------ 
Working with cursors
--------------------
Cursors can access the columns in the tabular data in a variety of ways.  One of the easiest to follow is to simply refer to the columns by the order in which they appear when they were retrieved.  This is fine if one writes scripts in the same way.  In the example that follows, the list of fields to be used with the cursor operations is defined as:
in_flds = ['OID@', 'SHAPE@X', 'SHAPE@Y', 'Int_fld', 'Float_fld', 'String_fld']

 

When using the above notation, the position of the fields is used to reference their values.   So you may see code that uses ' for row in cursor '  with row[0] being the feature object id (OID@) and row[3] being the value from an integer field (Int_fld).  If you are like me, anything beyond 2, means you are finger counting remembering the python counting is zero-based.  I now prefer to spend the extra time assigning variable names rather than using positional notation.  You can see this in lines 12-13 below.

in_fc = r'C:\Folder\path_to\A_Geodatabase.gdb\FeatureClass   # or Table

desc = arcpy.Describe(in_fc)
SR = desc.spatialReference
in_flds = ['OID@', 'SHAPE@X', 'SHAPE@Y', 'Int_fld', 'Float_fld', 'String_fld']
where_clause = None
spatial_reference = SR
explode_to_points = True
sql_clause = (None, None)

results = []
with arcpy.da.SearchCursor(in_tbl, in_flds) as cursor:
    for id, x, y, i_val, f_val, s_val in cursor:
        if id > 10:
            do stuff
        else:
            do other stuff
        results.append(... put the stuff here ...)
return results

-----------------------------------------------------------------------------
arcgisscripting
---------------
arcgisscripting can be located in your ArcGIS Pro distribution once everything is imported (arcgisscripting.__file__).  It is located in the installation path (substitute C:\\ArcPro for your Pro path, everything else is the same)
'C:\\ArcPro\\bin\\Python\\envs\\arcgispro-py3\\lib\\site-packages\\arcgisscripting.pyd'
Now importing arcpy, also imports parts of arcgisscripting and it also imports the geoprocessor from
C:\ArcPro\Resources\ArcPy\arcpy\geoprocessing\__init__
which imports _base.py which uses the Geoprocessor class as 'gp'

dir(arcgisscripting)
['ExecuteAbort', 'ExecuteError', 'ExecuteWarning', 'NumPyArrayToRaster', 'Raster', 'RasterToNumPyArray', '__cleanup__', '__doc__', '__file__', '__loader__', '__name__', '__package__', '__spec__', '_addTimeInterval', '_analyzeForSD', '_chart', '_convertWebMapToMapDocument', '_createGISServerConnectionFile', '_createGeocodeSDDraft', '_createMapSDDraft', '_createimageservicesddraft', '_getImageEXIFProperties', '_getUTMFromLocation', '_hasLocalFunctionRasterImplementation', '_listDateTimeStringFormats', '_listStyleItems', '_listTimeZones', '_mapping', '_ss', '_wrapLocalFunctionRaster', '_wrapToolRaster', 'arcgis', 'create', 'da', 'getmytoolboxespath', 'getsystemtoolboxespath', 'getsystemtoolboxespaths', 'na']

dir(arcgisscripting.da)
['Describe', 'Domain', 'Editor', 'ExtendTable', 'FeatureClassToNumPyArray', 'InsertCursor', 'ListDomains', 'ListFieldConflictFilters', 'ListReplicas', 'ListSubtypes', 'ListVersions', 'NumPyArrayToFeatureClass', 'NumPyArrayToTable', 'Replica', 'SearchCursor', 'TableToNumPyArray', 'UpdateCursor', 'Version', 'Walk', '__doc__', '__loader__', '__name__', '__package__', '__spec__', '_internal_eq', '_internal_sd', '_internal_vb'
References
----------

 

More later...

The background story of syntax highlighting

and other tales from behind the scenes...

: ---- before table ---- :

A small snippet of code

---- before code ----

# comment line place using html...
import sys
import numpy as np
#
a = [1, 2, 3, 4]
print(("{} "*len(a)).format(*a))

---- what it really looks like ----

..... Now I had to break up the lines into bits in order for people to read it... so bear with me.

.....  Do remember, someone wrote these rules and the syntax stuff, so don't blame the

...... messenger


<pre class="language-python line-numbers">

<code>

<span class="token comment"># comment line place using html...</span><br />

<span class="token keyword">import</span>

 sys<br />

<span class="token keyword">import</span>

 numpy

<span class="token keyword">as</span>

 np<br />

<span class="token comment">#</span><br />

a

<span class="token operator">=</span>

<span class="token punctuation">[</span>

<span class="token number">1</span>

<span class="token punctuation">,</span>

<span class="token number">2</span>

<span class="token punctuation">,</span>

<span class="token number">3</span>

<span class="token punctuation">,</span>

<span class="token number">4</span>

<span class="token punctuation">]</span><br />

<span class="token keyword">print</span>

<span class="token punctuation">(</span>

<span class="token punctuation">(</span>

<span class="token string">"{} "</span>

<span class="token operator">*</span>

len

<span class="token punctuation">(</span>

a

<span class="token punctuation">)</span>

<span class="token punctuation">)</span>

<span class="token punctuation">.</span>

format

<span class="token punctuation">(</span>

<span class="token operator">*</span>

a

<span class="token punctuation">)</span>

<span class="token punctuation">)</span>

&zwj;&zwj;&zwj;&zwj;&zwj;&zwj;&zwj;&zwj;&zwj;&zwj;&zwj;&zwj;&zwj;&zwj;&zwj;&zwj;&zwj;&zwj;

&zwj;&zwj;&zwj;&zwj;

<span class="line-numbers-rows">

<span>&zwj;</span>

<span>&zwj;</span>

<span>&zwj;</span>

<span>&zwj;</span>

<span>&zwj;</span>

<span>&zwj;</span>

</span>

</code>

</pre>

 

OooooK  so a lot of stuff that makes editing them nearly impossible.  It must be hard to identify punctuation, i guess... hence all the bloat (paid by the letter?)

 

Now this is plain text formatting..... no fancy colors... not too bad

import numpy as np
#
a = [1, 2, 3, 4]
print(("{} "*len(a)).format(*a))

 

Here is what it looks like as plain text....

<pre class="language-none line-numbers">

<code>import numpy as np<br />

#<br />

a = [1, 2, 3, 4]<br />

print(("{} "*len(a)).format(*a))

<span class="line-numbers-rows"><span>&zwj;</span><span>&zwj;</span><span>&zwj;</span><span>&zwj;</span></span></code></pre>

 

You can almost read it ... much better... there will be less need for **ml programmers... a small sacrifice, but I am sure they can be redeployed.  

Sadly, though there is still all that line number stuff... so I am going to try a direct copy into the html section ( see the <> button up top) so see if I can dump it and reduce the bloat

 

import numpy as np
#
a = [1, 2, 3, 4]
print(("{} "*len(a)).format(*a))

 

So the price you pay for syntax highlighting and line numbering is buried in the code.  So, if you don't start syntax highlighting first then paste your code when you are there... it is way easier just to start again or not post your code snippet for perusal at all.

All done...

Look

 

Code as verbose as possible.  Based on a generic function that I use.

Rules

  1. Do not have Pro or Map open
  2. Run your code
  3. Open the project when you are done.
  4. Review rules 1-3
  5. You can sort using multiple fields or a single field.  You can reverse sort by just slicing and reversing at the same time... arr_sort = arr_sort[::-1] way quicker than trying to get all the options on a dialog.  This can be done on a column or multiple columns basis.
  6. Still too confusing? .... attach the code to a new script in a toolbox... it doesn't need any parameters, but you have to substitute the paths.  ExtendTable requires a close tie to Pro even if your external python IDE is the editing IDE default when you set up Pro.

ExtendTable will do a permanent join if the script is run from a tool within ArcToolbox, making it a nice alternative to the standard join procedure.  

 

Numpy has functionality within the recfunctions module located at https://github.com/numpy/numpy/blob/master/numpy/lib/recfunctions.py 

These would include

__all__ = [
'append_fields', 'drop_fields', 'find_duplicates',
'get_fieldstructure', 'join_by', 'merge_arrays',
'rec_append_fields', 'rec_drop_fields', 'rec_join',
'recursive_fill_fields', 'rename_fields', 'stack_arrays',
]

So you have options to do your work within or without Arc* depending on what is the easiest in the workflow

 

Less than 1000 words too.

Let's get to the point... in this case, the points of a square.

 

We can convert it to a structured array using nothing but arcpy methods, a mix of arcpy and numpy and more recently, and extension into Pandas.

 

Here is what they look like using various functions.

-----------------------------------------------------------------------------------------------------------------------------------------------------------

1.  Rudimentary steps

in_fc = r"C:\Git_Dan\a_Data\testdata.gdb\square"
flds = ['SHAPE@X', 'SHAPE@Y']
cur = arcpy.da.SearchCursor(in_fc, flds, None, None, True, (None, None))
a = cur._as_narray()

# ----- the reveal ----
array([(342000.0, 5022000.0), (342000.0, 5023000.0), (343000.0, 5023000.0),
       (343000.0, 5022000.0),  (342000.0, 5022000.0)],
      dtype=[('SHAPE@X', '<f8'), ('SHAPE@Y', '<f8')])

 

We can also do it directly using another arcpy.da interface and it yields the same results perhaps there is a direct link between the cursor's _as_narray and FeatureClassToNumPyArray, details are buried within a *.pyd.  Nonetheless, the journey reaches the same destination.

flds = ['SHAPE@X', 'SHAPE@Y']
arcpy.da.FeatureClassToNumPyArray(in_fc, field_names=flds, explode_to_points=True)

array([(342000.0, 5022000.0), (342000.0, 5023000.0), (343000.0, 5023000.0),
       (343000.0, 5022000.0), (342000.00000000186, 5022000.0)],
      dtype=[('SHAPE@X', '<f8'), ('SHAPE@Y', '<f8')])

 

In its absolute simplest form, the searchcursor simply needs a list of fields and a flag to convert a polygon to points.  The _as_narray method handles conversion to a numpy recarray with a specified dtype (2 floating point numbers) and a shape ( (5,) ) indicating that there are 5 pairs of numbers forming the polygon

 

-----------------------------------------------------------------------------------------------------------------------------------------------------------

2.  To an ndarray with just the coordinates

The _xy yields a simpler form...

_xy(in_fc)

    array([[  342000.,  5022000.],
           [  342000.,  5023000.],
           [  343000.,  5023000.],
           [  343000.,  5022000.],
           [  342000.,  5022000.]])
... because it uses arcpy's searchcursor and the _as_narray method to speedily get the geometry (lines 5-8).  At this point, the array's data type is converted (as a 'view') to 64-bit floating point numbers and reshaped to 5 rows of 2 points each.
def _xy(in_fc):
    """Convert featureclass geometry (in_fc) to a simple 2D point array.
    :  See _xyID if you need id values.
    """
    flds = ['SHAPE@X', 'SHAPE@Y']
    args = [in_fc, flds, None, None, True, (None, None)]
    cur = arcpy.da.SearchCursor(*args)
    a = cur._as_narray()
    N = len(a)
    a = a.view(dtype='float64')
    a = a.reshape(N, 2)
    del cur
    return a

 

Working with either form of the returned array has its advantages and disadvantages.  I won't cover those here, but it is worthy to note, that I regularly use both forms of array.

 

-----------------------------------------------------------------------------------------------------------------------------------------------------------

3.  Throw in the IDs

Now, we can go back to working with a structured array if we need to maintain different data types in our columns.  This can be accomplished using _xyID which returns the same point values with the addition of the identification number of the polygon that it came with (ie. polygon with objectid 1).  Rather than the cumbersome datatype field descriptions that we got with _xy, I decided to provide simpler names by specifying a dtype whose names were different.  This can be done as long as you don't try to alter the underlying data type (more on this later).

 

def _xyID(in_fc, to_pnts=True):
    """Convert featureclass geometry (in_fc) to a simple 2D structured array
    :  with ID, X, Y values. Optionally convert to points, otherwise centroid.
    """
    flds = ['OID@', 'SHAPE@X', 'SHAPE@Y']
    args = [in_fc, flds, None, None, to_pnts, (None, None)]
    cur = arcpy.da.SearchCursor(*args)
    a = cur._as_narray()
    a.dtype = [('IDs', '<i4'), ('Xs', '<f8'), ('Ys', '<f8')]
    del cur
    return a

# ---- yields ----
array([(1, 342000.0, 5022000.0), (1, 342000.0, 5023000.0), (1, 343000.0, 5023000.0),
       (1, 343000.0, 5022000.0), (1, 342000.0, 5022000.0)],
      dtype=[('IDs', '<i4'), ('Xs', '<f8'), ('Ys', '<f8')])

 

 

-----------------------------------------------------------------------------------------------------------------------------------------------------------

4.  How about some attributes

Yes, but what about the attributes.  Lets get fancy and use _ndarray which will pull in all the data from the table as well.

 

def _ndarray(in_fc, to_pnts=True, flds=None, SR=None):
    """Convert featureclass geometry (in_fc) to a structured ndarray including
    :  options to select fields and specify a spatial reference.
    :
    :Requires:
    :--------
    : in_fc - input featureclass
    : to_pnts - True, convert the shape to points. False, centroid returned.
    : flds - '*' for all, others: 'Shape',  ['SHAPE@X', 'SHAPE@Y'], or specify
    """
    if flds is None:
        flds = "*"
    if SR is None:
        desc = arcpy.da.Describe(in_fc)
        SR = desc['spatialReference']
    args = [in_fc, flds, None, SR, to_pnts, (None, None)]
    cur = arcpy.da.SearchCursor(*args)
    a = cur._as_narray()
    del cur
    return a

array([ (1, [342000.0, 5022000.0], 342000.0, 5022000.0, 343000.0, 5023000.0, 4000.0, 1000000.0),
       (1, [342000.0, 5023000.0], 342000.0, 5022000.0, 343000.0, 5023000.0, 4000.0, 1000000.0),
       (1, [343000.0, 5023000.0], 342000.0, 5022000.0, 343000.0, 5023000.0, 4000.0, 1000000.0),
       (1, [343000.0, 5022000.0], 342000.0, 5022000.0, 343000.0, 5023000.0, 4000.0, 1000000.0),
       (1, [342000.0, 5022000.0], 342000.0, 5022000.0, 343000.0, 5023000.0, 4000.0, 1000000.0)],
      dtype=[('OBJECTID', '<i4'), ('Shape', '<f8', (2,)), ('EXT_MIN_X', '<f8'),
             ('EXT_MIN_Y', '<f8'), ('EXT_MAX_X', '<f8'), ('EXT_MAX_Y', '<f8'),
             ('Shape_Length', '<f8'), ('Shape_Area', '<f8')])

 

I suppose that you cleverly noticed that the other attributes are replicated during the process of converting the polygon to points.  This can also be accomplished with arcpy.da.FeatureClassToNumPyArray simply by using '*' to use all fields rather than prescribed ones. 

 

-----------------------------------------------------------------------------------------------------------------------------------------------------------

5.  Got to have both, but not together

Working with an array which has attributes and point geometry AND the attributes are replicated is a tad of a pain, so it is time to split the attributes from the geometry.

 

def _two_arrays(in_fc, both=True, split=True):
    """Send to a numpy structured/array and split it into a geometry array
    :  and an attribute array.  They can be joined back later if needed.
    """
    a = _xyID(in_fc, to_pnts=True)  # just the geometry
    shp_fld, oid_fld, SR, shp_type = fc_info(in_fc)
    dt_a = [('IDs', '<i4'), ('Xs', '<f8'), ('Ys', '<f8')]  # option 1
    dt_b = [('IDs', '<i4'), ('Xc', '<f8'), ('Yc', '<f8')]
    a.dtype = dt_a
    b = None
    if split:
        ids = np.unique(a['IDs'])
        w = np.where(np.diff(a['IDs']))[0] + 1
        a = np.split(a, w)
        a = np.array([[ids[i], a[i][['Xs', 'Ys']]] for i in range(len(ids))])
    if both:
        b = _ndarray(in_fc, to_pnts=False)
        dt_b.extend(b.dtype.descr[2:])
        b.dtype = dt_b
    return a, b



(array([[1, array([(342000.0, 5022000.0), (342000.0, 5023000.0),
          (343000.0, 5023000.0), (343000.0, 5022000.0),
          (342000.0, 5022000.0)],
        dtype=[('Xs', '<f8'), ('Ys', '<f8')])]], dtype=object),
array([ (1, 342500.0, 5022500.0, 342000.0, 5022000.0, 343000.0, 5023000.0,
            4000.0, 1000000.0)],
        dtype=[('IDs', '<i4'), ('Xc', '<f8'), ('Yc', '<f8'), ('EXT_MIN_X', '<f8'),
                  ('EXT_MIN_Y', '<f8'), ('EXT_MAX_X', '<f8'), ('EXT_MAX_Y', '<f8'),
                  ('Shape_Length', '<f8'), ('Shape_Area', '<f8')]))

This incarnation of getting at the data produces an 'object array' which is an array that can contain records of different size and/or composition.  Why?  Simple!  If you are converting polygons to point objects, then it is highly unlikely that you will have the same number of points per polygon.  You have two options then... collapse and rearrange the structure of the outputs so that all the coordinates are arranged in their respective columns or keep the inherent structure the data.  This is exemplified by the outputs returned by FeatureClassToNumPyArray and the above function.  Lines 24-27 represent the object array which contains an index number then an array of point coordinates with a specified dtype.  Lines 28-32 is the attribute array for that polygon.  As a nice side bonus, the point coordinates have been reduced to the polygon centroid value and not returned as individual points with repetitive attributes as noted previously.  Two for the price of one.

 

-----------------------------------------------------------------------------------------------------------------------------------------------------------

6.  What's this new thing

 

Then there is the new SpatialDataFrame and the tie in to Pandas rather than numpy.  It offers some advantages in terms of analysis but in some other areas not so much.  To get the simple polygon out to points you can do the following.

 

from arcgis import SpatialDataFrame as SDF

a = SDF.from_featureclass(in_fc)

p = a.SHAPE.tolist()[0]

p['rings']
[[[342000, 5022000], [342000, 5023000], [343000, 5023000], [343000, 5022000], [342000, 5022000]]]

 

There are many options available to get to the root geometry of your geometry and the ability to work with arcpy, numpy and now Pandas opens a variety of opportunities for improving your analytical arsenal.

 

Next blog will focus on the analysis side and how the options can be used to your advantage.

You have 12 months of data in some raster form... You want some statistical parameter... There are areas of nodata, the extents are all the same... and ... you have a gazillion of these to do.

Sounds like you have a 'cube'... the original 'space-time' cube' .  You can pull space from a time slice... You can slice time through space.  At every location on a 'grid' you have Z as a sequence over time. 

 

Here is the code.   I will use ascii files, but they don't have to be, you just need to prep your files before you use them.

 

Source data originally from .... here .... thanks Xander Bakker.

 

import os
import numpy as np
from textwrap import dedent, indent
import arcpy

arcpy.overwriteOutput = True

ft = {'bool': lambda x: repr(x.astype('int32')),
      'float': '{: 0.3f}'.format}
np.set_printoptions(edgeitems=3, linewidth=80, precision=2, suppress=True,
                    threshold=80, formatter=ft)
np.ma.masked_print_option.set_display('-')  # change to a single -

# ---- Processing temporal ascii files ----
# Header information
# ncols 720
# nrows 360
# xllcorner -180
# yllcorner -90
# cellsize 0.5
# NODATA_Value -9999
# ---- Begin the process ----
#
cols = 720
rows = 360
ll_corner =  arcpy.Point(-180., -90.0)  # to position the bottom left corner
dx_dy = 0.5
nodata = '-9999'
#
# ---- create the basic workspace parameters, create masked arrays ----
#
out_file = r'c:\Data\ascii_samples\avg_yr.tif'
folder = r'C:\Data\ascii_samples'
arcpy.env.workspace = folder
ascii_files = arcpy.ListFiles("*.asc")
a_s = [folder + '\{}'.format(i) for i in ascii_files]
arrays = []
for arr in a_s:
    a = np.mafromtxt(arr, dtype='int', comments='#',
                     delimiter=' ', skip_header=6,
                     missing_values=nodata, usemask=True)
    arrays.append(a)
#
# ---- A sample calculation from the inputs.... calculate the mean ----
#
N = len(arrays)                    # number of months... in this case
arr_shp = (N,) + arrays[0].shape   # we want a (month, col, row) array
msk = arrays[0].mask               # clone the mask... assume they are the same
e = np.zeros(arr_shp, 'int')       # one way is create a zeros array and fill
for i in range(len(arrays)):
    e[i] = arrays[i]
a = np.ma.array(e, mask=e*msk[np.newaxis, :, :])  # the empty array is ready
#
# ---- your need here... ie. Calculate a mean ----
#
avg = a.mean(axis=0)  # calculate the average down through the months
#
# ---- send it out to be viewed in ArcMap or ArcGIS Pro ----
#
value_to_nodata = int(avg.get_fill_value())
out = avg.view(np.ndarray, fill_value=value_to_nodata)
g = arcpy.NumPyArrayToRaster(out, ll_corner, dx_dy, dx_dy)
g.save(out_file)

So the basic process is simple.  I have coded this verbosely and used input parameters read manually from the ascii header since it is the quickest way.... and you would probably know what they are from the start.

 

So in this example... 12 months of some variable, averaged accounting for the nodata cells.  Do the map stuff, define its projection... project it, give it some symbology and move on.

I will leave that for those that make maps.

Modify to suit... maybe I will get this into a toolbox someday

 

 

NOTE.....

 

Now in the linked example, there was a need to simply convert those to rasters from the input format.  In that case you would simply consolidate the salient portions of the script as follows and create the output rasters within the masked array creation loop

......
for arr in a_s:
    a = np.mafromtxt(arr, dtype='int', comments='#',
                     delimiter=' ', skip_header=6,
                     missing_values=nodata, usemask=True)
    value_to_nodata = int(a.get_fill_value())
    out = a.view(np.ndarray, fill_value=value_to_nodata)
    r = arcpy.NumPyArrayToRaster(out, ll_corner, dx_dy, dx_dy)
    out_file = arr.replace(".asc", ".tif")
    r.save(out_file)
    del r, a

.....

 

So for either doing statistical calculations for temporal data, or for format conversion... there are options available where arcpy and numpy play nice.

Geometry...

That is what I am normally interested in.  The baggage (aka, attributes) tag along for the ride and I normally find it easier to keep the baggage separate until I am done with the geometry. For those following along, see my previous post.

 

Let us compare some of the ways that we can pull the geometry out of a featureclass.  The following demonstrations can be followed in your own workflow for testing your specific cases.

 

Imports first ___________________________________________________________________________________

import sys
import numpy as np
import arcpy
import arcgis

_______________________________________________________________________________________________

I prefer to import modules in the order of less polluting first to ensure namespace is grabbed by what I want in order of importance in case there is any conflict during import.

 

Searchcursor __________________________________________________________________________________

# ---- pick a featureclass ----
fc0 = r'drive:\your_spaceless_folder_path\your_geodatabase.gdb\polygon'

# ---- using arcpy.Describe ----
desc = arcpy.Describe(fc0)
shp = desc.shapeFieldName    # ---- using dot notation ----

# ---- using arcpy.da.Describe ----
desc = arcpy.da.Describe(fc0)
shp = desc['shapeFieldName'] # ---- note... now a dictionary ----

# ---- geometry, as X, Y ----
flds = [shp + "@"]
shps = [row[0] for row in arcpy.da.SearchCursor(fc0, flds)]

_______________________________________________________________________________________________

 

Which of course you can begin to use to get basic properties.  In this example,

Geometry properties  ____________________________________________________________________________

for s in shps:
    print(s.type, s.partCount, s.pointCount, s.length, s.area)
polygon 1 5 40.0 100.0
polygon 1 10 64.0 64.0
polygon 2 14 99.41640786499875 182.0

_______________________________________________________________________________________________

 

Then you can get to work and convert to a NumPy array quickly and simply.  Since I know this is a polygon featureclass, it only takes a couple of lines to perform the task.

 

Get to the point(s) ______________________________________________________________________________

pnts = []
for shp in shps:
    for j in range(shp.partCount):
        pt = shp.getPart(j)
        p_list = [(pnt.X, pnt.Y) for pnt in pt if pnt]
        pnts.extend(p_list)
dt = [('Xs', '<f8'), ('Ys', '<f8')]
a = np.asarray(pnts, dtype=dt)

_______________________________________________________________________________________________

The basic difference between the array forms is how you want to work with them.

In the example above array 'a' has a specified data type (dtype).  The fields/columns of the array can be accessed by name (Xs and Ys).  Since this array is a structured array, the access would be a['Xs'] and a['Ys'].  If I converted this to a record array, one could use object.field notation.

 

The coordinates are nothing more than the same type of number, so the field names can be dispensed with altogether.  In this case, the sentient being is responsible for knowing what they are working with.  Both forms of the same data are shown below

_______________________________________________________________________________________________

a # a structured array
array([(300020.0, 5000000.0), (300010.0, 5000000.0), (300010.0, 5000010.0), ...,
          (300002.0, 5000002.0), (300008.0, 5000002.0), (300005.0, 5000008.0)],
      dtype=[('Xs', '<f8'), code="" ('ys',="" '<f8')])

a_nd = np.asarray(pnts)  # as an np.ndarray
array([[ 300020.00,  5000000.00],
       [ 300010.00,  5000000.00],
       [ 300010.00,  5000010.00],
       ...,
       [ 300002.00,  5000002.00],
       [ 300008.00,  5000002.00],
       [ 300005.00,  5000008.00]])
a_nd.dtype
dtype('float64')

_______________________________________________________________________________________________

 

The following demo function can be used on your data to examine some of the options and explore some of the properties and methods available to each.  Just don't forget the imports at the start of the post

 

The demo  _____________________________________________________________________________________

def _demo(in_fc):
    """Run the demo and return some objects
    : create a SpatialDataFrame class object from a featureclass
    : create a record array from a_sdf
    : create a numpy.ndarray using the da module
    """

    SR = arcpy.Describe(in_fc).SpatialReference

    sdf = arcgis.SpatialDataFrame
    a_sdf = sdf.from_featureclass(in_fc,
                                  sql_clause=None,
                                  where_clause=None,
                                  sr=SR)
    a_rec = sdf.to_records(a_sdf)  # record array
    a_np = arcpy.da.FeatureClassToNumPyArray(in_fc,
                                             field_names="*",
                                             spatial_reference=SR,
                                             explode_to_points=True)
    a_np2 = fc_g(in_fc)
    return sdf, a_sdf, a_rec, a_np, a_np2

 

Now... behind all the magic of the above options, a searchcursor is behind the acquisition of all of the options shown above.  The options do, however, provide access to methods and properties which are unique to their data class.  In many instances these are shared.

 

Here is what the 'look like' .  In the case of the SpatialDataFrame, a_sdf, the same when converted to a record array and the conventional da.FeatureClassToNumPyArray... all fields were included.  in the last option, a_np2, just the geometry is returned as demonstrated in the above examples, with the addition of handling the geometry parts and point when the polygon is 'exploded' to it's constituent points.

 

a_sdf
Out[40]:
   OBJECTID  Id Text_fld                                              SHAPE
0         1   1     None  {'rings': [[[300020, 5000000], [300010, 500000...
1         2   2        C  {'rings': [[[300010, 5000020], [300010, 500001...
2         3   3        A  {'rings': [[[300010, 5000010], [300010, 500002...

a_rec
Out[41]:
rec.array([ (0, 1, 1, None, {"rings": [[[300020, 5000000], [300010, 5000000], [300010, 5000010], [300020, 5000010], [300020, 5000000]]], "spatialReference": {"latestWkid": 2951, "wkid": 2146}}),
(1, 2, 2, 'C', {"rings": [[[300010, 5000020], [300010, 5000010], [300000, 5000010], [300000, 5000020], [300010, 5000020]], [[300002, 5000018], [300002, 5000012], [300008, 5000012], [300008, 5000018], [300002, 5000018]]], "spatialReference": {"latestWkid": 2951, "wkid": 2146}}),
(2, 3, 3, 'A', {"rings": [[[300010, 5000010], [300010, 5000020], [300020, 5000020], [300020, 5000010], [300010, 5000010]], [[300010, 5000010], [300010, 5000000], [300000, 5000000], [300000, 5000010], [300010, 5000010]], [[300005, 5000008], [300002, 5000002], [300008, 5000002], [300005, 5000008]]], "spatialReference": {"latestWkid": 2951, "wkid": 2146}})],
          dtype=[('index', '<i8'), ('OBJECTID', '<i8'), ('Id', '<i8'), ('Text_fld', 'O'), ('SHAPE', 'O')])

a_np
Out[42]:
array([(1, [300020.0, 5000000.0], 1, 'None', 40.0, 100.0),
       (1, [300010.0, 5000000.0], 1, 'None', 40.0, 100.0),
       (1, [300010.0, 5000010.0], 1, 'None', 40.0, 100.0), ...,
       (3, [300002.0, 5000002.0], 3, 'A', 99.41640786499875, 182.0),
       (3, [300008.0, 5000002.0], 3, 'A', 99.41640786499875, 182.0),
       (3, [300005.0, 5000008.0], 3, 'A', 99.41640786499875, 182.0)],
      dtype=[('OBJECTID', '<i4'), ('Shape', '<f8', (2,)), ('Id', '<i4'), ('Text_fld', '<U255'), ('Shape_Length', '<f8'), ('Shape_Area', '<f8')])

a_np2
Out[43]:
array([(1, 0, 300020.0, 5000000.0), (1, 0, 300010.0, 5000000.0),
       (1, 0, 300010.0, 5000010.0), ..., (3, 1, 300002.0, 5000002.0),
       (3, 1, 300008.0, 5000002.0), (3, 1, 300005.0, 5000008.0)],
      dtype=[('ID_num', '<i4'), ('Part_num', '<i4'), ('Xs', '<f8'), ('Ys', '<f8')])

 

When the SpatialDataFrame is converted to a numpy record array, the geometry field (Shape) has a dtype of 'O', which is an object array.  This will be the 'normal' case since it is unlikely that all polygons will contain the same number of parts and points per part.  To truly be a numpy array, the 'shape' of the array needs to be consistent.  

 

The latter two representations, (a_np and a_np2) deal with this but converting the polygons to points.  These points can be converted back to polygons after they are used in some process such as moving, calculating parameters, reprojecting.

 

Next installation... working with geometry in its various representations.

ArcGIS API for Python  ... version 1.2.0

Another cleverly named product to provide more clarity to the other named kindred... ArcGIS for Desktop, ArcGIS Pro, Arc... etcetera.  The link to the help is here.  The ability to work with Jupyter notebooks, NumPy, SciPy and Arcpy is touted and welcomed (...and there is something about web mapping and stuff as well).

 

Where stuff is

Locating ArcGIS in your installation path, depends on how you installed it... for a single user (aka no sharesies) or for anyone.  This describes the single user situation.

To begin, import the module and check its __path__ and __file__ property.  My installation path is in an ArcPro folder that I created... yours will differ, but beyond the base folder, everything should be the same.

 

Basic import______________________________________________________________________

In [1]: import arcgis
In [2]: arcgis.__path__
Out[2]: ['C:\\ArcPro\\bin\\Python\\envs\\arcgispro-py3\\lib\\site-packages\\arcgis']

In [3]: arcgis.__file__

Out[3]: 'C:\\ArcPro\\bin\\Python\\envs\\arcgispro-py3\\lib\\site-packages\\arcgis\\__init__.py'

_________________________________________________________________________________

 

The __init__.py file is read during the import and a whole load of imports are done to cover access to pretty well all available modules contained in the ArcGIS path.

 

If you only want to work with the geometry or Pandas related functionality, you can import the SpatialDataFrame class directly.

 

SpatialDataFrame___________________________________________________________________

In [1]: from arcgis import SpatialDataFrame as sdf

# ---- or ----

In [1]: from arcgis.features._data.geodataset import SpatialDataFrame as sdf

In [2]: sdf.__module__  # to confirm the path

Out[2]: 'arcgis.features._data.geodataset.geodataframe'

_________________________________________________________________________________

 

The SpatialDataFrame class is pretty well all that is in geodataframe.py script.

Another useful call within that class is one to from_featureclass and to_featureclass which can be obtained by a variety of other imports or by a direct call to the io module's fileops.py

 

Featureclass access________________________________________________________________

In [3]: from arcgis.features._data.geodataset.io import from_featureclass, to_featureclass

# ---- or ----

In [4] from arcgis.features._data.geodataset import io

_________________________________________________________________________________

 

If you prefer the module.function method of calling, the option in ln [4] can be used to convert a featureclass to a SpatialDataFrame.

 

The Reveal________________________________________________________________________

In [5] fc0 = r'drive:\your_space_free_path\your_geodatabase.gdb\polygon'

In [6] a = io.from_featureclass(fc0)

In [7] print("\nSpatialDataFrame 'a' \n{}".format(a))  # ---- the sdf contents ----

 

SpatialDataFrame 'a'
   OBJECTID  Id Text_fld                                              SHAPE
0         1   1     None  {'rings': [[[300020, 5000000], [300010, 500000...
1         2   2        C  {'rings': [[[300010, 5000020], [300010, 500001...
2         3   3        A  {'rings': [[[300010, 5000010], [300010, 500002...

ln [8] type(a)

Out[8]: <class 'arcgis.features._data.geodataset.geodataframe.spatialdataframe'>

_________________________________________________________________________________

 

Behind the Scenes

Great so far... BUT when you delve deeper into what is going on under the hood, you will find out that the from_featureclass method...

  1. imports the SpatialDataFrame class
  2. uses arcpy.da.SearchCursor to convert the featureclass into lists of
    1. values (non-geometry fields)
    2. geometry (geometry field)
  3. the value fields are then converted to a pandas dataframe (pd.DataFrame)
  4. the dataframe in step 3 is then used with the geometry values to produce a SpatialDataFrame modelled after a Pandas dataframe.

 

Final Comments

So... in conclusion, an arcpy.da searchcursor is used to get all the necessary geometry and attribute data then ultimately to a SpatialDataFrame ... which is like a Pandas dataframe ... but it is geometry enabled.

Sort of like geopandas... (geopandas on GitHub) but relying on arc* functionality and code for the most part.

 

More to Come

There will be more posts as I have time... the next post... Geometry.... part II

Requirements:  You need to have ArcGIS Pro installed and the associated installation of jupyter notebooks as described in previous blogs. 

 

This is a novel way of presenting all things python, arcpy and eventually through the Python API for ArcGIS Pro...

 

Background link   Arcgis Pro 2... Jupyter Notebook Setup Basics

 

The link here...distance.pdf .... should open up the explanatory pdf for the accompanying jupyter notebook.  If not... click the attachment below.

 

Unfortunately, Jive isn't at the stage to allow for this type of interactive content.

 

If you have comments, questions or found a bug, send them to me via email.

In the last blog post... ArcGIS Pro 2... Creating Desktop Shortcuts .. I showed how to create desktop shortcuts to some of the more commonly used applications within esri's distribution of python.  In this post, the Jupyter Notebook will be addressed.

 

Before you start on this venture, make sure that you have read up on notebooks and see if they are for you and your workflow and not something that '... everyone is doing it, and so should I.... '.  They do have their place, but they require maintenance.

 

The first order of business.... 

  • Create a folder location to store your notebooks... it is just way easier to keep them in one location and distribute them from there.  I chose to put them in my local machine's GitHub folder, in a separate folder within. ... C:\GitDan\JupyterNoteBooks ... pretty clever ehhh?
  • Right-click on the file ---- "C:/ArcPro/bin/Python/envs/arcgispro-py3/Scripts/jupyter-notebook-script.py" ---- and select create shortcut  which will simply create a shortcut to that file.
    • Go to the Shortcut tab and edit it the Target line and put ...
    • ---- C:\ArcPro\bin\Python\envs\arcgispro-py3\pythonw.exe ---- in front what is there... this will yield the whole path to python in the Pro distribution, plus the script to run (yes, the paths are that long and cryptic).
    • C:\ArcPro\bin\Python\envs\arcgispro-py3\pythonw.exe "C:/ArcPro/bin/Python/envs/arcgispro-py3/Scripts/jupyter-notebook-script.py"
  • In the Start in: line, put the path to the folder that you are going to house your notebooks.  In my example, this was the folder ---- C:\Git_Dan\JupyterNoteBooks ----
  • Finally, right-click on the shortcut, select Copy, flip over to your Desktop and Paste it there.  
  • yes... I know you can go to the command interface and run the command line from there, but why.  You can also use Anaconda Navigator in other non-ArcGIS Pro environments.  The installation and setup of the application within the Pro environment isn't worth the extra effort.

 

 

Python runs that script, which imports notebook.notebookapp as shown below.  That import does the work of setting up notebooks to work in your target folder.

""" Source code of ---- jupyter-notebook-script.py ---- located in
:  _drive_:\_folder_\bin\Python\envs\arcgispro-py3\Scripts ... where _drive_ and _folder
:  are unique to your machine... for example
:  ---- C:\ArcPro\bin\Python\envs\arcgispro-py3\Scripts ----
"""

if __name__ == '__main__':
    import sys
    import notebook.notebookapp

    sys.exit(notebook.notebookapp.main())

 

The details aren't really relevant... Just click on the shortcut, create your notebook and begin sharing... but before you do that, begin your reading here

 

And from a recent post... perhaps the future might look like this...