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

Py... blog

80 posts

The aggravation of aggregation ...

 

               Raster with some classes                                                 Aggregation using maximum

aggregation_demo.pngaggregation_demo2.png

No Spatial Analyst??? A raster is just a big array... hmmmm. I wonder.

What is available at all license levels?

Not on the list?  There are examples of how to get other geometries to array format elsewhere. But everyone knows about old school ascii.  Splurge on disk space... ascii files are easy to edit.  There are many options available to get data into array format.

 

Code section
ASCII example...

Demo code...

This simple sample code should give you some ideas.  In this version, edge considerations are not accounted for, so should that be an issue, you can pad the array with an appropriate collar.

"""
Script:  aggregate_demo.py
Modified: 2016-01-30
Author:  Dan.Patterson AT  carleton.ca
Purpose:  To demonstrate aggregation of raster data without the
spatial analyst extension.  A sample raster is created and methods
to convert an array to  a raster and vice versa are shown.
Notes:
- RasterToNumPyArray(in_raster, {lower_left_corner},
                        {ncols}, {nrows}, {nodata_to_value})
  arr = arcpy.RasterToNumPyArray(rast)
- NumPyArrayToRaster(in_array, {lower_left_corner}, {x_cell_size},
                    {y_cell_size}, {value_to_nodata})
  rast = arcpy.NumPyArrayToRaster(a, arcpy.Point(300000,5025000),
                                  10,10,-9999)
  rast.save(r"F:\Demos\raster_ops\test_agg")
    # esri grid, or add tif, jpg etc
"""
import numpy as np
from numpy.lib.stride_tricks import as_strided
np.set_printoptions(edgeitems=3,linewidth=80,precision=2,
                    suppress=True,threshold=50)
from textwrap import dedent

import arcpy
arcpy.env.overwriteOutput = True
def block_a(a, block=(3,3)):
    """
Provide a 2D block view of a 2D array. No error checking made.
    Columns and rows outside of the block are truncated.
    """
    a = np.ascontiguousarray(a)
    r, c = block
    shape = (a.shape[0]/r, a.shape[1]/c) + block
    strides = (r*a.strides[0], c*a.strides[1]) + a.strides
    b_a = as_strided(a, shape=shape, strides=strides)
    return b_a
def agg_demo(n):
    """
Run the demo with a preset array shape and content.
       See the header.
    """
    a = np.random.random_integers(0,high=5,size=n*n).reshape((n,n))
    rast = arcpy.NumPyArrayToRaster(a,x_cell_size=10)
    agg_rast = arcpy.sa.Aggregate(rast,2,"MAXIMUM")
    agg_arr = arcpy.RasterToNumPyArray(agg_rast)
    # --- a_s is the strided array, a_agg_max is the strided array max
    a_s  = block_a(a, block=(2,2))
    a_agg_max = a_s.max(axis=(2,3))
    # ---
    frmt = "\nInput array..\n{}\n\n" \
          "Arcpy.sa aggregate..\n{}\n\n" \
          "Numpy aggregate..\n{}\n\n" \
          "All close? {}"
    yup = np.allclose(agg_arr,a_agg_max)
    print(dedent(frmt).format(a, agg_arr, a_agg_max, yup))
    return a, agg_arr, a_s, a_agg_max

if __name__=="__main__":
    """ Returns the input array, it's aggregation raster from
    arcpy.sa, the raster representation of the raster and the
    block representation and the aggregation array.
    """

    n=10
    a, agg_arr, a_s, a_agg_max = agg_demo(n)

 

This is a sample ascii file so that you can

see the numeric inputs and structure of

the ascii header and its data.

ncols         10
nrows         10
xllcorner     300000
yllcorner     5025000
cellsize      10
NODATA_value  -9999
0 0 2 5 -1 3 5 2 5 0
-1 -1 2 1 1 -1 -1 4 0 4
4 5 5 5 1 4 1 1 2 2
1 1 3 1 2 0 5 -1 2 3
1 1 5 2 4 5 4 2 5 -1
2 1 0 -1 2 3 2 1 5 1
3 1 5 5 0 3 -1 3 -1 2
1 -1 0 5 2 5 2 1 -1 2
5 4 4 5 2 3 0 5 4 0
1 1 3 5 4 -1 4 5 1 5

 

Does it work with big rasters?

Does it support other summary statistics?

Of course, I covered the details of raster statistics in an earlier post.

 

That's all for now...

I have been working with arrays and lists a lot lately and I wanted to view them in column mode rather than in row mode to make documentation easier to follow and to make the output more intuitative.  The demo relies on the numpy module, but that is no issue since everyone has it with their ArcMap and ArcGIS Pro installation.

You can alter where an array gets parsed by changing the values in this line...

 

np.set_printoptions(edgeitems=2,linewidth=80,precision=2,suppress=True,threshold=8)

 

I altered the [ and ] characters common in lists and arrays, just to show you could do it.  I also added an indentation option.

If you like the array output better than the list output, you can simply make an array from a list using ... new_array np.array(input_list) ...

 

"""
Script:    array_print_demo.py
Author:   
Dan.Patterson@carleton.ca

Modified:  2016-01-16
Purpose:   Demonstration of formatting arrays/lists in alternate formats
Functions:
  -  frmt_arr(arr, indents=0)
  -  indent_arr(arr, indents=1)  called by above
"""

import numpy as np
np.set_printoptions(edgeitems=2,linewidth=80,precision=2,suppress=True,threshold=8)

def frmt_arr(arr, indents=0):
    """
    Format arrays or lists with dimensions >=3 printed by rows rather
    than the usual column expression.  Indentation can be employed as well as
    the number of indent levels.  See string.expandtabs to alter the tab level.
    """

    if isinstance(arr, list):
        arr = np.asarray(arr)  # make sure inputs are array
    result=[]; temp = []
    arr_dim = arr.ndim
    if arr_dim < 3:
        if (arr_dim == 1) and (len(arr[0]) > 1):
            arr = arr.reshape((arr.shape[-1],1))
        if indents:
            arr = indent_arr(arr, indents)
        return arr
    elif arr_dim == 3:
        temp.append(arr)
    elif arr_dim > 3:
        temp = [sub for sub in arr]   
    for arr in temp:
        for x in zip(*arr):
            result.append("   ".join(map(str,x))) # use tabs \t, not space
        if arr_dim > 3:
            result.append("----")
        out = "\n".join(result)
    if indents:
        out = indent_arr(out, indents)
        #out = (str(out).replace("[","|")).replace("]","|")
        #tabs = "    "*indents # "\t"
        # out = tabs + out.replace("\n","\n" + tabs)
    return out
def indent_arr(arr, indents=1):
    """
    Add an indent to a str or repr version of an array.
    The number of indents is determined by the indents option.
    """
       
    out = (str(arr).replace("[","|")).replace("]","|")
    tabs = "    "*indents # "\t"
    out = tabs + out.replace("\n","\n" + tabs)
    return out
if __name__=='__main__':
    """largely a demonstration of slicing by array dimensions
    """

    # syntax  frmt_arr(arr, indents=1)
    a = np.random.randint(1,1000,16)
    a1 = a.tolist()
    e = a.reshape((4,2,2)); e1 = d.tolist()

 

Your preference

>>> print(a1)
[[135, 944], [196, 335], [761, 521], [529, 687], [803, 393], [254, 797], [610, 605], [328, 516]]
>>> # or ...
>>> print(frmt_arr(a1,indents=1))
    ||135 944|
     |196 335|
     ...,
     |610 605|
     |328 516||

A 3D list or array
>>> print(e)
[[[135 944]
  [196 335]]
[[761 521]
  [529 687]]
[[803 393]
  [254 797]]
[[610 605]
  [328 516]]]
>>> # or ...
>>> print(frmt_arr(e,indents =2))
        |135 944|   |761 521|   |803 393|   |610 605|
        |196 335|   |529 687|   |254 797|   |328 516|

Date: 2015-10-08   Modified: 2017-04-19 new ***

Included:

  • Add X Y coordinates to table by distance or percent along poly* features  ***
  • Distance to a specific point
  • Cumulative distance
  • Inter-point distance
  • Azimuth to a point
  • Angle between consecutive points
  • Convert Azumith to compass bearing   ***
  • Convert Degrees decimal minutes to decimal degrees ***

 

Purpose:

Many posts on GeoNet come from people looking to write a script or produce a tool which can be handled in a much simpler fashion.  This is first in a series of posts that will address how to use the field calculator to its fullest.  The use, and perhaps the limitations to using, the field calculator should be guided by the following considerations:

  1. you need to perform some rudimentary task for which there is no existing tool (ie.  it is so simple to do...what is the point of having it builtin)
  2. your datasets are relatively small ...
    • of course this is a subjective classification...  so if processing takes more than a minute... you have left the realm of small
  3. you don't have to perform the task that often...
    • basically you are trying to fix up an oversight is the initial data requirements, or the requirements have changed since project inception
    • you got the data elsewhere and you are trying to get it to meet project standards
    • you want to explore
  4. you have no clue how to script but want to dabble as a prelude to full script development.
  5. whatever else I left out ....

 

Notes:

  • Make sure you have projected data and a double field or nothing will make sense.  I care not about performing great circle calculations or determining geodetic areas.  The field calculator...IMHO... is not the place to do that.
  • Ensure that the sequence of points is indeed correct.  If you have been editing or fiddling around with the data, make sure that nothing is amiss with the data.  What you see, is often not what is.  You can get erroneous results from any environment
  • Sorting will not change the outcome, the results will be in feature order.
    • If you want them in another order, then use tools to do that.
    • I will be paralleling geometric calculations using NumPy and Arcpy is a separate series which removes this issue.
  • Selections can be used but they will be in feature order.
    • sorting and selecting can get it to look like you want it ... but it doesn't make it so... want and is,  are two separate concepts like project and define projection (but I digress)
  • VB anything is not covered.
    • It is time to move on to the next realm in the sequence and evolution of languages.
    • in ArcGIS Pro, VB isn't even an option, so get used to it
  • code is written verbosely, for the most part
    • I will try to use the simplest of programming concepts...  simplicity trumps coolness and speed.
    • parallel, optional and/or optimized solutions will be discussed elsewhere.
  • the math module is already imported

 

Useage:

For all functions, do the following:

  • select the Python parser ... unfortunately, it can't be selected by default
  • toggle on, Show code block
  • Pre-logic Script Code:
    • paste the code in this block
  • Expression box
    • paste the function call in the expression box
    • ensure that there is preceding space before the expression...it will generate an error otherwise
    • ie  dist_between(!Shape!)     ... call the dist_between function ... the shape field is required
    • field names are surrounded by exclamation marks ( ! ) ... the shapefile and file geodatabase standard

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

Add X or Y coordinate to table field by distance or percentage

For use in tables, to retrieve the X or Y coordinates in an appropriate field.  See the header for requirements.


Pre-logic Script Code:

def pnt_along(shape, value=0.0, use_fraction=False, XorY="X"): 
    """Position X or Y coordinate, x/y meters or decimal fraction along a line.
    :Requires:
    :--------
    : shape field: python parser use !Shape!
    : value: (distance or decimal fraction, 0-1)
    : use_fraction: (True/False)
    : XorY: specify X or Y coordinates
    :
    :Returns: the specified coordinate (X or Y) meters or % along line or boundary
    :-------
    :
    :Useage: pnt_along(!Shape!, 100, False, "X") # X, 100 m from start point
    :
    """
    XorY = XorY.upper()
    if use_fraction and (value > 1.0):
        value = value/100.0
    if shape.type.lower() == "polygon":
        shape = shape.boundary()
    pnt = shape.positionAlongLine(value,use_fraction)
    if XorY == 'X':
        return pnt.centroid.X
    else:
        return pnt.centroid.Y

expression =

pnt_along(!Shape!, 100, False, "X")  # X, 100 m from start point

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

Distance to a specific point

Calculate the distance to a specific point within a dataset.   For example, you can determine the distance from a point cloud's center to other points.  Features not in the file can be obtained by other means.  Useful to get summary information on inter-point distances as a prelude to a full clustering or KD-means study.

For example, the potential workflow to get the distance of every point to the point clouds center might include the following steps:

  • add X and Y style fields to contain the coordinates
  • use Field Statistics to get the mean center (use a pen to record them... not everything is automagic)
  • add a Dist_to field to house the calculations
  • use the script below

 

Pre-logic Script Code:

"""  dist_to(shape, from_x, from_y)
input:      shape field, origin x,y
returns:    distance to the specified point
expression: dist_to(!Shape!, x, y)
"""

def dist_to(shape, from_x, from_y):
    x = shape.centroid.X
    y = shape.centroid.Y
    distance = math.sqrt((x - from_x)**2 + (y - from_y)**2)
    return distance

expression =

dist_to(!Shape!, x, y)

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

Cumulative Distance

To determine the distance from a point in a data table to every other point in sequence.  Essentially a sequential perimeter calculation.

 

Pre-logic Script Code:

""" input shape field: returns cumulative distance between points
dist_cumu(!Shape!) #enter into the expression box"""
x0 = 0.0
y0 = 0.0
distance = 0.0
def dist_cumu(shape):
    global x0
    global y0
    global distance
    x = shape.centroid.X
    y = shape.centroid.Y
    if x0 == 0.0 and y0 == 0.0:
        x0 = x
        y0 = y
    distance += math.sqrt((x - x0)**2 + (y - y0)**2)
    x0 = x
    y0 = y
    return distance

expression =

dist_cum(!Shape!)

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

Inter-point distance

Determine distances between sequential point pairs (not cumulative).

 

Pre-logic Script Code:

""" dist_between(shape)
input: shape field
returns: distance between successive points
expression: dist_between(!Shape!)
"""

x0 = 0.0
y0 = 0.0
def dist_between(shape):
    global x0
    global y0
    x = shape.centroid.X
    y = shape.centroid.Y
    if x0 == 0.0 and y0 == 0.0:
        x0 = x
        y0 = y
    distance = math.sqrt((x - x0)**2 + (y - y0)**2)
    x0 = x
    y0 = y
    return distance

expression =

dist_between(!Shape!)

 

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

Azimuth to a specific point

Determine the azimuth to a specific point, for example, a point cloud's center.

 

Pre-logic Script Code:

"""  azimuth_to(shape, from_x, from_y)
input:      shape field, from_x, from_y
returns:    angle between 0 and <360 between a specified point and others
expression: azimuth_to(!Shape!, from_x, from_y)
"""

def azimuth_to(shape, from_x, from_y):
    radian = math.atan((shape.centroid.X - from_x)/(shape.centroid.Y - from_y))
    degrees = math.degrees(radian)
    if degrees < 0:
        return degrees + 360.0
    else:
        return degrees

expression =

azimuth_to(!Shape!,from_x, from_y)

azimuth_1.png

 

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

Angle between successive points

Determine the angle between points, for example, angle changes between waypoints.

 

Pre-logic Script Code:

"""  angle_between(shape)
input:      shape field
returns:    angle between successive points,
            NE +ve 0 to 90, NW +ve 90 to 180,
            SE -ve <0 to -90, SW -ve <-90 to -180
expression: angle_between(!Shape!)
"""

x0 = 0.0
y0 = 0.0
angle = 0.0
def angle_between(shape):
    global x0
    global y0
    x = shape.centroid.X
    y = shape.centroid.Y
    if x0 == 0.0 and y0 == 0.0:
        x0 = x
        y0 = y
        return 0.0
    radian = math.atan2((shape.centroid.Y - y0),(shape.centroid.X - x0))
    angle = math.degrees(radian)
    x0 = x
    y0 = y
    return angle

expression =

angle_between(!Shape!)

 

-----------------------------------------------------------------------------------------------------------------------------------------
Line direction or Azimuth to Compass Bearing

 

 

Can be used to determine the direction/orientation between two points which may or may not be on a polyline.  Provide the origin and destination points.  The origin may be the 0,0 origin or the beginning of a polyline or a polyline segment.

def line_dir(orig, dest, fromNorth=False):
    """Direction of a line given 2 points
    : orig, dest - two points representing the start and end of a line.
    : fromNorth - True or False gives angle relative to x-axis)
    :
    """

    orig = np.asarray(orig)
    dest = np.asarray(dest)
    dx, dy = dest - orig
    ang = np.degrees(np.arctan2(dy, dx))
    if fromNorth:
        ang = np.mod((450.0 - ang), 360.)
    return ang

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

Convert Azimuth to Compass Bearing

 

Once the Azimuth to a particular point is determined, this can be converted to a compass direction, centered in 22.5 degree bands.  The type of compass can be altered to suit... see script header.

 

Pre-logic Script Code:

import numpy as np
global a
global c
c = np.array(['N', 'NNE', 'NE', 'ENE', 'E', 'ESE', 'SE', 'SSE',
              'S', 'SSW', 'SW', 'WSW', 'W', 'WNW', 'NW', 'NNW', 'N'])
a = np.arange(11.25, 360., 22.5)
def compass(angle):
    """Return the compass direction based on supplied angle.
    :Requires:
    :--------
    : angle - angle(s) in degrees, no check made for other formats.
    : - a single value, list or np.ndarray can be used as input.
    : - angles are assumed to be from compass north, alter to suit.
    :
    :Returns: The compass direction.
    :-------
    :
    :Notes:
    :-----
    : Compass ranges can be altered to suit the desired format.
    : See various wiki's for other options. This incarnation uses 22.5
    : degree ranges with the compass centered on the range.
    : ie. N between 348.75 and 11.25 degrees, range equals 22.5)
    :
    :----------------------------------------------------------------------
    """

    if isinstance(angle, (float, int, list, np.ndarray)):
        angle = np.atleast_1d(angle)
    comp_dir = c[np.digitize(angle, a)]
    if len(comp_dir) == 1:
        comp_dir[0]
    return comp_dir

expression =

compass(!Azimuth_Field_Name!)  # python parser, field name enclosed in quotes

 

:--------------------------------------------------------------------------------------------------------------------------------------------------------------

Degrees decimal minutes to decimal degrees

def ddm_ddd(a, sep=" "):
    """ convert decimal minute string to decimal degrees
    : a - degree, decimal minute string
    : sep - usually a space, but check
    """

    d, m = [abs(float(i)) for i in a.split(sep)]
    sign = [-1, 1][d < 0]
    dd = sign*(d + m/60.)
    return dd

 

For example   ddm_ddd(!YourStringField!, sep=" ") will convert a string/text field into a double in your new field

This post seemed interesting for several reasons:   Iterate and Add Values from many fields

  • it is generally the type of thing that we know how to do in a spreadsheet easily
  • it is a different thought process that one uses to translate into code
  • I got to experiment with block based functions that weren't related to raster or image data

 

Reflecting before problem solving

As a preamble, let's examine the year in different ways.

Shaping up your year in various ways....

the conventional year, by day

 

[[  1,   2,   3, ..., 363, 364, 365]] 

the year by average month

      [[  1,   2,   3, ...,  28,  29,  30], 

       [ 31,  32,  33, ...,  58,  59,  60],

       [ 61,  62,  63, ...,  88,  89,  90],

       ...,   ... snip ....

       [271, 272, 273, ..., 298, 299, 300],

       [301, 302, 303, ..., 328, 329, 330],

       [331, 332, 333, ..., 358, 359, 360]]

the year by weeks

      [[  1,   2,   3, ...,   5,   6,   7], 

       [  8,   9,  10, ...,  12,  13,  14],

       [ 15,  16,  17, ...,  19,  20,  21],

       ...,  ... snip ....

       [344, 345, 346, ..., 348, 349, 350],

       [351, 352, 353, ..., 355, 356, 357],

       [358, 359, 360, ..., 362, 363, 364]]

how about the average lunar cycle?

      [[  1,   2,   3, ...,  26,  27,  28],

       [ 29,  30,  31, ...,  54,  55,  56],

       [ 57,  58,  59, ...,  82,  83,  84],

       ...,

       [281, 282, 283, ..., 306, 307, 308],

       [309, 310, 311, ..., 334, 335, 336],

       [337, 338, 339, ..., 362, 363, 364]])

 

Your year can be shaped in different ways... but notice that the year is never always equal to 365 (not talking about leap years here).

The year is divisible by chunks... the fiddly-bits get truncated.  Notice the different year lengths in the above (365, 364, 360, 364).

Did you ever wonder why your birthday is on a different sequentially changing day every year?

 

On to the problem

So here is some data.  It is constructed so that you can do the 'math' yourself.  The human brain is pretty good at picking out the patterns...it is just less adept at formulating a plan to implement it.  So here goes..

 

The task is to take a set of data and determine some properties on a block by block basis.  In this example, the maximum needs to be determined for blocks of 4 values for each block in a row.  Coincidently ... (aka, to simplify the problem) this yields 5 blocks of 4 in each row.  This type of block is denoted as a 1x4 block... one row by 4 columns, if you wish.

 

1.  Input data.... subdivide the data into 1x4 blocks

[[  1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20]

[ 21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40]

[ 41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60]

[ 61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80]

[ 81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100]]

 

2.  Imagine that the data are split into 5 blocks of 1x4 values by 5 rows which means that you would have 25 1x4 blocks of values.   Here is an example [1, 2, 3,4].  Cleverly we can detect that the maximum as being 4 which most of you know how to derive in some programming language using some interface, whether it be the field calculator or a tool in arctoolbox.

 

3.  Using your minds eye, or a spreadsheet, determine the maximum for each block of 4, within each row, for all rows.

The maximum of each 1x4 block in each row, would yield.

[[  4   8  12  16  20]

[ 24  28  32  36  40]

[ 44  48  52  56  60]

[ 64  68  72  76  80]

[ 84  88  92  96 100]]

 

4.  The final stage of the cited problem was to determine the maximum of each of the 1x4 blocks within each row.  That is, the maximum of the values produced in step 2.

The results on a row basis is as follows:

 

[ 20  40  60  80 100]

 

: ---- Example 2 ----
5.  Let's do the sum of each 1x4 block and determine the row max

     Again, I won't show the intermediate step since it is visually confusing unless you are interested in it.

 

6. sums of each 1x4 block in each row

[[[[ 10  26  42  58  74]

  [[ [ 90 106 122 138 154]

  [[ [170 186 202 218 234]

  [[ [250 266 282 298 314]

  [[ [330 346 362 378 394]]

 

7. maximum of (6) by row        [ 74 154 234 314 394]

 

: ---- Example 3 ----

8.  Do the sum but with a masked array. The mask locations should

:    not be included in calculations

[[1 2 3 4 5 -- 7 8 9 10 11 -- 13 14 15 16 17 -- 19 20]

[[ [21 22 23 -- 25 26 27 28 29 -- 31 32 33 34 35 -- 37 38 39 40]

[[[41 -- 43 44 45 46 47 -- 49 50 51 52 53 -- 55 56 57 58 59 --]

[[ [61 62 63 64 65 -- 67 68 69 70 71 -- 73 74 75 76 77 -- 79 80]

[[[81 82 83 -- 85 86 87 88 89 -- 91 92 93 94 95 -- 97 98 99 100]]

 

9. now with masked values

 

10. sum of the 1x 4 blocks accounting for the mask

[[10 20 30 58 56]

[[ [66 106 92 102 154]

[[ [128 138 202 164 174]

[[ [250 200 210 298 236]

[[ [246 346 272 282 394]]

 

11. maximum of 10 by row      [58 154 202 298 394]

# -------------------------------------------------------------------------------------------------------------------------

Now, this can obviously be extended to look at longer data lists.  The following example looks at how to determine the maximum monthly value and maximum sum over a 5 year period.  To simplify the example and to facilitate mental math, a year has 360 days, and 1 month has 30 days.

---- the trick ----

a = int_range((5,360),1,1)  # rows, columns, start, step
b = block_reshape(a,block=(1,30))

 

:---- Example 4 ----

1.  Input data.... subdivide the data into 1x30 blocks... each 360 day year, by 30 days per month

[[   1    2    3         ...,  358  359  360]

[ 361  362  363    ...,  718  719  720]

[ 721  722  723    ..., 1078 1079 1080]

[1081 1082 1083 ..., 1438 1439 1440]

[1441 1442 1443 ..., 1798 1799 1800]]

 

2.  Now imagine what that would look like if subdivided

3.  maximum of each 1x30 block in each row

[[  30   60   90  120  150  180  210  240  270  300  330  360]

[ 390  420  450  480  510  540  570  600  630  660  690  720]

[ 750  780  810  840  870  900  930  960  990 1020 1050 1080]

[1110 1140 1170 1200 1230 1260 1290 1320 1350 1380 1410 1440]

[1470 1500 1530 1560 1590 1620 1650 1680 1710 1740 1770 1800]]

 

4.  maximum of (3) by row

[ 360  720 1080 1440 1800]

 

: ---- Example 5 ----

5.  Let's do the sum of each 1x30 block and determine the row max

6. sums of each 1x30 block in each row

[[  465  1365  2265  3165  4065  4965  5865  6765  7665  8565  9465 10365]

[11265 12165 13065 13965 14865 15765 16665 17565 18465 19365 20265 21165]

[22065 22965 23865 24765 25665 26565 27465 28365 29265 30165 31065 31965]

[32865 33765 34665 35565 36465 37365 38265 39165 40065 40965 41865 42765]

[43665 44565 45465 46365 47265 48165 49065 49965 50865 51765 52665 53565]]

7. maximum of (6) by row

[10365 21165 31965 42765 53565]

 

Summary:

Still not convinced?  Well what is the cumulative sum of 1 to 30 inclusive:

[  1,   3,   6,  10,  15,  21,  28,  36,  45,  55,  66,  78,  91, 105, 120, 136,
    153, 171, 190, 210, 231, 253, 276, 300, 325, 351, 378, 406, 435, 465 ]

How about the next 30 days:

[  31,   63,   96,  130,  165,  201,  238,  276,  315,  355,  396,  438,  481,  525,  570,

        616,  663,  711,  760,  810,  861,  913,  966, 1020, 1075, 1131, 1188, 1246, 1305, 1365]

Now notice the first 2 entries in (6) above... 465 and 1365

 

---- What is the magic??? -----

>>> a.shape   # 5 years, 360 days in a year

(5, 360)

>>> b.shape   # 5 years, 12 months, 30 days per 1 month

(5, 12, 1, 30)

>>> c.shape   # reshape by month

(5, 12)

>>> d.shape   # number of years

(5,)

>>> c1.shape  # reshape the 5 years by month and sum each month, take the max

(5, 12)

>>> d1.shape  # sums for the 5 years, take the yearly max

(5,)

>>>

 

And finally in code

 

let's make up 5 years of data where it rains 1mm per day, day in-day out
>>> # here we go
>>> x = np.ones((5, 360), dtype='int32')
>>> x = np.ones((5, 360), dtype='int32')
>>> x.shape
(5, 360)
>>> x.min(), x.max()
(1, 1)
>>> # so far so good
>>>
>>> x1 = block_reshape(x, block=(1, 30)) # produce the monthly blocks
>>> x2 = ((x1.T).max(axis=0))            # start reorganizing
>>> x2 = (x2.T).squeeze()                # more reorganizing
>>> x3 = x2.max(axis=1)                  # get the maximum per year
>>> #
>>> x4 =((x1.T).sum(axis=0))             # do the sums
>>> x4 = (x4.T).squeeze()                # finish reorganizing
>>> x5 = x4.max(axis=1)                  # and finish up with the maximums
>>>
>>> x   # the daily precipitation, for 360 days for 5 years
array([[1, 1, 1, ..., 1, 1, 1],
       [1, 1, 1, ..., 1, 1, 1],
       [1, 1, 1, ..., 1, 1, 1],
       [1, 1, 1, ..., 1, 1, 1],
       [1, 1, 1, ..., 1, 1, 1]])
>>> x2  # The monthly maximum for 5 years (12 months for 5 years)
array([[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]])
>>> x3  # the yearly summary for the maximums for 5 years
array([1, 1, 1, 1, 1])
>>> x4  # the monthly sums  12 months for 5 years
array([[30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30],
       [30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30],
       [30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30],
       [30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30],
       [30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30]])
>>> x5  # the annual montly maximum for each year
array([30, 30, 30, 30, 30])
>>>

 

That's all for now

I had made reference before of using numpy and arcpy together to perform particular tasks.  Sometimes the combination of the two enables the user to do things they wouldn't otherwise be able to do.  This ditty presents a very simple way to sort a dataset and produce a permantly sorted version while maintaining the original feature identification numbers.  So to begin with, I will remind/introduce the reader to the arcpy data access module's functions as well as comparable ones in numpy.

 

On the arcmap side, the user should become familiar with:

  • arcpy.Describe
  • arcpy.Sort
  • arcpy.da.FeatureClassToNumPyArray and
  • arcpy.da.NumPyArrayToFeatureClass. 

From the  numpy side:

  • np.setprintoptions
  • np.arange and
  • np.sort

The latter performs the same work as arcpy.Sort, but it doesn't have any license restrictions as the Sort tool does.  It also allows you to produce a sequential list of numbers like the oft used Autoincrement field calculator function.

 

The script will perform the following tasks:

  • import the required modules,
  • specify the input file,
  • the spatial reference of the input,
  • the fields used to sort on, in the order of priority,
  • a field to place the rank values, and
  • specify the output file.

 

One can add the field to place the ranks into ahead of time or add a line of code to do it within the script (see arcpy.AddField_management).

 

# imports ......

import numpy as np
import arcpy

np.set_printoptions(edgeitems=4,linewidth=75,precision=2,
                    suppress=True,threshold=10)

arcpy.env.overwriteOutput = True

# inputs  ......

input_shp = r'c:\!Scripts\Shapefiles\RandomPnts.shp'
all_fields = [fld.name for fld in arcpy.ListFields(input_shp)]
SR = arcpy.Describe(input_shp).spatialReference

sort_fields = ['X','Y']
order_field = 'Sort_id'

output_shp = r'c:\!Scripts\Shapefiles\SortedPnts.

# outputs ......

arr = arcpy.da.FeatureClassToNumPyArray(input_shp, "*", spatial_reference=SR)
arr_sort = np.sort(arr, order=sort_fields)

arr_sort['Sort_id'] = np.arange(arr_sort.size)

arcpy.da.NumPyArrayToFeatureClass(arr_sort, output_shp, ['Shape'], SR)


# summary .....

print('\nInput file: {}\nField list: {}'.format(input_shp,all_fields))

print('Sort by: {},  ranks to: {}'.format(sort_fields,orderfield))

 

result

 

Input file: c:\!Scripts\Shapefiles\RandomPnts.shp

Field list: [u'FID', u'Shape', u'Id', u'X', u'Y', u'Sort_id']

Sort by: ['X', 'Y'],  ranks to: Sort_id

 

Sort_np_arc.png 

Give it a try with your own files.

It is really a pain that certain highly used functions are only available an advanced license level.  This is an alternate to the options of using Excel to produce a pivot table from ArcMap tabular data.

 

Flipping rows and columns in data generally works smoothly when the table contains one data type, whether it be integer, float or text.  Problems arise when you add stuff to Excel is that it allows you do so without regard to the underlying data.  So, columns get mixed data types and rows do as well. Uniformity by column is the rule.

 

In NumPy, each column has a particular data type.  The data type controls the operations that can be performed on it.  Numeric fields can have all the number type operations used...similarly for string/text fields.  It is possible to cast a field as an "object" type allowing for mixed type entries.  The nice thing about this type, is that you can't really do anything with it unless it is recast into a more useful form...but it does serve as a conduit to other programs or just for presentation purposes.

 

In the following example, the line

a = # your array goes here

can be derived using

a = arcpy.FeatureClasstoNumPyArray(....)  FeatureClassToNumPyArray

 

The nature of np.unique change in version 1.9 to get the number of unique classes as well.  So if you are using ArcGIS Pro, then you can use the newer version if desired by simply changing line 04 below.

 

a_u, idx, counts = np.unique(a_s, return_inverse=True, unique_counts=True)

 

Array conversion to summary table or pivot tableInput and output

Well... who needs an advanced license or excel ...

Assume we have an array of the format shown in the Input section.  We can determine the counts or sums of unique values in a field, using the following.

  • sort the array on a field,
  • get unique values in that field,
  • sum using the values in another field as weights
  • rotate if desired
    import numpy as np
    a = # your array goes here
    a_s = a[np.argsort(a, order="Class")]
    a_u, idx = np.unique(a_s["Class"], return_inverse=True)
    bin_cnt = np.bincount(idx,weights=a_s['Count'])
    ans = np.array((a_u, bin_cnt), dtype='object')
    print("a_u\n{}\nidx {}\nanswer\n{}".format(a_u, idx, ans))
    rot90 = np.rot90(ans, k=1) 
    and_flipud = np.flipud(rot90) #np.flipud(np.rot90(a,k=1))))
    frmt = ("pivot table... rotate 90, flip up/down\n{}" 
    print(frmt.format(and_flipud))

 

The trick is to set the data type to 'object'. You just use FeatureClassToNumPyArray or TableToNumPyArray and their inverses to get to/from array format.  Ergo....pivot table should NOT be just for an advanced license

For all-ish combos, you can just add the desired lines to the above

 

for i in range(4):
    print("\nrotated {}\n{}".format(90*i, np.rot90(a, k=i)))
for i in range(4):
    f = "\nrotate by {} and flip up/down\n{}"
    print(f.format(90*i, np.flipud(np.rot90(a, k=i))))
for i in range(5):
    f = "\nrotate by {} and flip left/right\n{}"
    print(f.format(90*i, np.fliplr(np.rot90(a, k=i))))

Input table with the following fields

'ID', 'X', 'Y', 'Class', 'Count'

 

>>> input array...
[[(0, 6.0, 0.0, 'a', 10)]
[(1, 7.0, 9.0, 'c', 1)]
[(2, 8.0, 6.0, 'b', 2)]
[(3, 3.0, 2.0, 'a', 5)]
[(4, 6.0, 0.0, 'a', 5)]
[(5, 2.0, 5.0, 'b', 2)]
[(6, 3.0, 2.0, 'a', 10)]
[(7, 8.0, 6.0, 'b', 2)]
[(8, 7.0, 9.0, 'c', 1)]
[(9, 6.0, 0.0, 'a', 10)]]

>>> # the results

a_u  # unique values
['a' 'b' 'c']
idx [0 0 0 0 0 1 1 1 2 2]

answer # the sums
[['a' 'b' 'c']
[40.0 6.0 2.0]]

pivot table... rotate 90, flip up/down
[['a' 40.0]
['b' 6.0]
['c' 2.0]]


This post was inspired by the GeoNet thread... https://community.esri.com/thread/116880

 

I'm a student and I need a python script that i can use for ArcMap

 

I usually suggest that my students use Modelbuilder to build workflows, export to a python script, then modify the script for general use with the existing, or other, data sets.  I personally don't use Modelbuilder, but I have used one of two methods to generate the needed workflow .... Method 1 will be presented in this post...method 2 follows.

 

Method 1

 

Do it once...get the script...modify and reuse


Because of the imbedded images...please open the *.pdf file to view the complete discussion...

 

Regards

Dan

This is part 2 which follows up on my previous blog post.  In this example, the assignment restrictions have changed and one must now develop a script from what they have read about Python and the tools that are used in everyday ArcMap workflows.

 

Details are given in the attached pdf of the same name.


Regards
Dan

 

Homework... make this into a script tool for use in arctoolbox

This blog is inspired by this thread https://community.esri.com/docs/DOC-2436#comment-9031 by Steve Offermann (Esri).  He suggested a very simple way to extend the capabilities of tool results and how to parse arguments for them.  I recommended the use of the Results window outputs in a previous blog.  Hats off to Steve.

 

I am only going to scratch the surface by presenting a fairly simple script...which could easily be turned into a tool.

In this example, a simple shapefile of hexagons, (presented in another blog post) was processed to yield:

 

  • an extent file, giving the bounds of each of the input hexagons,
  • the hexagon corners as points and sent to a shapefile, and,
  • the centroids of each hexagon was treated in a similar fashion

 

The whole trick is to parse your processes down into parameters that can be shared amongst tools.  In this case, tools that can be categorized as:

 

  • one parameter tools like:  AddXY_management and CopyFeatures_management
  • two parameter tools like:
    • FeatureEnvelopeToPolygon_management,
    • FeatureToPoint_management and
    • FeatureVerticesToPoints_management

 

This can then be amended by, or supplemented with, information on the input/output shape geometries.  I demonstrate this by calculating the X,Y coordinates for the point files.

 

So you are saying ... I don't do that stuff ... well remember, I don't do that webby stuff either.   Everyone has a different workflow and if my students are reading this, just think how you could batch project a bunch of files whilst simultaneously renaming them etc etc.  The imagination is only limited by its owner...

 

First the output....

 

Hexagon_outputs.png

 

And now the script....

"""
Script:   run_tools_demo.py
Author:   Dan.Patterson@carleton.ca

Source:   Stefan Offermann
Thread:   https://community.esri.com/docs/DOC-2436

Purpose:  Results window on steroids
  - take a polygon shapefile, determine its envelop,
  - convert the feature to centroids,
  - convert to feature points
  - calculate X,Y for all of the above
  - then make a back of everything
Requires:
  - a source file
  - an output folder
  - a list of tools to run
"""
import os
import arcpy
arcpy.env.overwriteOutput = True
in_FC = "c:/!BlogPosts/Runtools_Demo/Shapefiles/pointy_hex.shp"
path,in_File = os.path.split(in_FC)
path += "/"
backup = "c:/temp/shapefiles/"    # some output folder
# file endings
end = ["_env","_fp","_vert"]      # envelop, feature to point, feature vertices
# two argument tools
two_arg = ["FeatureEnvelopeToPolygon_management",
           "FeatureToPoint_management",
           "FeatureVerticesToPoints_management"
          ]
# one argument tools
one_arg =["AddXY_management", "CopyFeatures_management"]
#
outputs = [in_FC.replace(".shp", end[i]+".shp") for i in range(len(end))]
backups =  [outputs[i].replace(path, backup) for i in range(len(end))]
#
polygons = []
points = []
for i in range(len(two_arg)):                  # run the two argument tools
    args = [in_FC, outputs[i]]                 # select the output file
    result = getattr(arcpy, two_arg[i])(*args) # run the tool...and pray
    frmt = '\Processing tool: {}\n  Input: {}\n  Output: {}'
    print(frmt.format(tools[i], args[0], args[1]))
#
for i in range(len(outputs)):
    args = [outputs[i]]
    print(outputs[i], arcpy.Describe(outputs[i]).shapeType)
    if arcpy.Describe(outputs[i]).shapeType == 'Point':
        result = getattr(arcpy, one_arg[0])(*args) # calculate XY
    result_bak = getattr(arcpy, one_arg[1])(result, backups[i]) # backup
    print('Calculate XY for: {}'.format(result))
    print('Create Backups: {}\n  Output: {}'.format(result,result_bak))

Enjoy and experiment with your workflows.

UPDATE:   take the poll first before you read on How do you write Python path strings?

 

I am sure everyone is sick of hearing ... check your filenames and paths and make sure there is no X or Y.  Well, this is going to be a work in progress which demonstrates where things go wrong while maintaining the identity of the guilty.

Think about it
>>> import arcpy
>>> aoi = "f:\test\a"
>>> arcpy.env.workspace = aoi
>>> print(arcpy.env.workspace)
f: est
>>>

 

>>> print(os.path.abspath(arcpy.env.workspace))
F:\ est
>>> print(os.path.exists(arcpy.env.workspace))
False
>>> print(arcpy.Exists(arcpy.env.workspace))
False
>>>
>>> print("{!r:}".format(arcpy.env.workspace))
'f:\test\x07'
>>>

 

>>> os.listdir(aoi)
Traceback (most recent call last):
  File "<interactive input>", line 1, in <module>
OSError: [WinError 123] The filename, directory name,
or volume label syntax is incorrect: 'f:\test\x07'
>>>

 

>>> arcpy.ListWorkspaces("*","Folder")
>>>
>>> "!r:{}".format(arcpy.ListWorkspaces("*","Folder"))
'!r:None'
>>>

 

 

Examples... Rules broken and potential fixes

Total garbage... as well as way too long.  Time to buy an extra drive.

>>> x ="c:\somepath\aSubfolder\very_long\no_good\nix\this"
>>> print(x)                  # str notation
c:\somepath Subfolder ery_long
o_good
ix his
>>> print("{!r:}".format(x))  # repr notation
'c:\\somepath\x07Subfolder\x0bery_long\no_good\nix\this'
>>>
  • No r in front of the path.
  • \a \b \n \t \v are all escape characters... check the result
  • Notice the difference between plain str and repr notation

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

Solution 1... raw format

>>> x = r"c:\somepath\aSubfolder\very_long\no_good\nix\this"

>>> print(x)                  # str notation
c:\somepath\aSubfolder\very_long\no_good\nix\this

>>> print("{!r:}".format(x))  # repr notation
'c:\\somepath\\aSubfolder\\very_long\\no_good\\nix\\this'
>>>
  • Use raw formatting, the little r in front goes a long way.

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

Solution 2... double backslashes

>>> x ="c:\\somepath\\aSubfolder\\very_long\\no_good\\nix\\this"
>>> print(x)                  # str notation
c:\somepath\aSubfolder\very_long\no_good\nix\this

>>> print("{!r:}".format(x))  # repr notation
'c:\\somepath\\aSubfolder\\very_long\\no_good\\nix\\this'
>>>
  • Yes! I cleverly used raw formatting and everything should be fine but notice the difference between str and repr.

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

Solution 3... forward slashes

>>> x ="c:/somepath/aSubfolder/very_long/no_good/nix/this"
>>> print(x)                  # str notation
c:/somepath/aSubfolder/very_long/no_good/nix/this
>>> print("{!r:}".format(x))  # repr notation
'c:/somepath/aSubfolder/very_long/no_good/nix/this'
>>>

 

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

Solution 4... os.path functions

There are very useful functions and properties in os.path.  The reader is recommended to examine the contents after importing the os module (ie dir(os.path)  and help(os.path)

 

>>> x = r"F:\Writing_Projects\Before_I_Forget\Scripts\timeit_examples.py"
>>> base_name = os.path.basename(x)
>>> dir_name = os.path.dirname(x)
>>> os.path.split(joined)  # see splitdrive, splitext, splitunc
('F:\\Writing_Projects\\Before_I_Forget\\Scripts', 'timeit_examples.py')
>>> joined = os.path.join(dir_name,base_name)
>>> joined
'F:\\Writing_Projects\\Before_I_Forget\\Scripts\\timeit_examples.py'
>>>
>>> os.path.exists(joined)
True
>>> os.path.isdir(dir_name)
True
>>> os.path.isdir(joined)
False
>>>

ad nauseum

 

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

Gotcha's

Fixes often suggest the following ... what can go wrong, if you failed to check.

(1)

>>> x = "c:\somepath\aSubfolder\very_long\no_good\nix\this"
>>> new_folder = x.replace("\\","/")
>>> print(x)                  # str notation
c:\somepath Subfolder ery_long
o_good
ix his
>>> print("{!r:}".format(x))  # repr notation
'c:\\somepath\x07Subfolder\x0bery_long\no_good\nix\this'
>>>

 

(2)

>>> x = r"c:\new_project\aSubfolder\"
  File "<string>", line 1
    x = r"c:\new_project\aSubfolder\"
                                    ^
SyntaxError: EOL while scanning string literal

 

(3)

>>> x = "c:\new_project\New_Data"
>>> y = "new_grid"
>>> out = x + "\\" + y
>>> print(out)
c:
ew_project\New_Data\new_grid

 

(4)

>>> x = r"c:\new_project\New_Data"
>>> z = "\new_grid"
>>> out = x + z
>>> print(out)
c:\new_project\New_Data
ew_grid

 

(5)  This isn't going to happen again!

>>> x = r"c:\new_project\New_Data"
>>> z = r"\new_grid"
>>> out = x + y
>>> print(out)
c:\new_project\New_Datanew_grid

 

(6)  Last try

>>> x = r"c:\new_project\New_Data"
>>> z = r"new_grid"
>>> please = x + "\\" + z
>>> print(please)
c:\new_project\New_Data\new_grid

 

Well this isn't good!   Lesson?  Get it right the first time. Remember the next time someone says...

Have you checked your file paths...?????   Remember these examples.

 

Curtis pointed out this helpful link...I will include it here as well

Paths explained: Absolute, relative, UNC, and URL—Help | ArcGIS for Desktop

 

That's all for now.

I will deal with spaces in filenames in an update.  I am not even to go to UNC paths.

Code formatting tips

 

Updated - 2017/06/27  Added another reference and some editing.

 

This topic has been covered by others as well...

 

We all agree the Geonet code editor is horrible... but it has been updated.

Here are some other tips.

 

To begin... introduction or review

  • don't try to post code while you are responding to a thread in your inbox
  • access the More button, from the main thread's title... to do this:
    • click on the main thread's title
    • now you should see it... you may proceed

 

Step 1... select ... More ... then ... Syntax highlighter
  • Go to the question and press Reply ...
  • Select the Advanced editor if needed (or ...),  then select

If you can't see it, you didn't select it

 

 More...Syntax highlighter ,

 

Your code can now be pasted in and highlighted with the language of your choice .........

Your code should be highlighted something like this ............

--- Python -----------------------------

import numpy as np
a = np.arange(5)
print("Sample results:... {}".format(a))

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

Now the above example gives you the language syntax highlighting that you are familiar with..

Alternatives include just using the HTML/XML option

-----HTML/XML ---------------------

# just use the HTML/XML option.. syntax colors will be removed
import numpy as np
a = np.arange(5)
print("simple format style {}".format(a))
simple format style [0 1 2 3 4]

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

 

NOTE:   you can only edit code within this box and not from outside!

 

 

 

Script editing tipscont'd

HTML editing tips:....

  • You can get into the html editor for fine-tuning, but it can get ugly for the uninitiated.
  • Comments get deleted ... this not a full editor under your control
  • If you have lots of text to enter, do it first then enter and format code
  • If editing refresh is slow, use the HTML editor or you will have retired before it completes.
  • The editor seems to edit each character as you type and it gets painfully slower as the post gets bigger.
  • You can improve comments and code output by using tables like this and in the example below.

 

Here is a simple script with code and results published in columns (2 columns * 1 row).  If the contents are wider than the screen, the scroll-bar will be located at the end of the document rather than attached to each table (except for iThingys, then just use swipe).

 

Sample script using a plain format... 1920x1080px screen sizeResult 2
>>> import numpy as np
>>> a = np.arange(5)
>>> print("Sample results:... {}".format(a))
>>> # really long comment 30 |------40 | -----50 | -----60 | -----70 | ---- 80|
Sample results:... [0 1 2 3 4]
>>> # really long comment 30 |------40 | -----50 | -----60 | -----70 | ---- 80|

 

Leave space after a table so you can continue editing after the initial code insertion.

It is often hard to select the whitespace before or after a table and you may need to go to the html editor < > just above the More toggle

 

Larger script sample...

Before code tip:  try to keep your line length <70 characters

# -*- coding: UTF-8 -*-
"""
:Script:   demo.py
:Author:   Dan.Patterson@carleton.ca
:Modified: 2016-08-14
:Purpose:  none
:Functions:  help(<function name>) for help
:----------------------------
: _demo  -  This function ...
:Notes:
:References
:
"""

#---- imports, formats, constants ----

import sys
import numpy as np
from textwrap import dedent

ft = {'bool':lambda x: repr(x.astype('int32')),
      'float': '{: 0.3f}'.format}
np.set_printoptions(edgeitems=10, linewidth=80, precision=2, suppress=True,
                    threshold=100, formatter=ft)
script = sys.argv[0]

#---- functions ----

def _demo():
    """  
    :Requires:
    :--------
    :
    :Returns:
    :-------
    :
    """

    return None

#----------------------
if __name__ == "__main__":
    """Main section...   """
    #print("Script... {}".format(script))
    _demo()

 

Some space for editing after should be left since positioning the cursor is difficult after the fact.

 

Output options

 

  • You can paste text and graphics with a table column.
  • You can format a column to a maximum pixel size.

 

Sample output with a graph

Option 0: 1000 points
[[ 2. 2.]
[ 3. 3.]] extents
.... snip
Time results: 1.280e-05 s, for 1000 repeats

point_in_polygon.png

point_in_polygon.png

 

So there has been some improvement. 

Again...

You just have to remember that to edit code...

you have to go back to the syntax highlighter.

You can't edit directly on the page.

Reclassifying raster data can be a bit of a challenge, particularly if there are nodata values in the raster.  This is a simple example of how to perform classifications using a sample array. (background 6,000+ rasters the follow-up  )

 

An array will be used since it is simple to bring in raster data to numpy using arcpy's:

  •   RasterToNumPyArray
    •   RasterToNumPyArray (in_raster, {lower_left_corner}, {ncols}, {nrows}, {nodata_to_value})

and sending the result back out using:

  • NumPyArrayToRaster
    • NumPyArrayToRaster (in_array, {lower_left_corner}, {x_cell_size}, {y_cell_size}, {value_to_nodata})

On with the demo...

 

Raster with full data

old                                new

[ 0  1  2  3  4  5  6  7  8  9]    [1 1 1 1 1 2 2 2 2 2]

[10 11 12 13 14 15 16 17 18 19]    [3 3 3 3 3 4 4 4 4 4]

[20 21 22 23 24 25 26 27 28 29]    [5 5 5 5 5 6 6 6 6 6]

[30 31 32 33 34 35 36 37 38 39]    [7 7 7 7 7 7 7 7 7 7]

[40 41 42 43 44 45 46 47 48 49]    [7 7 7 7 7 7 7 7 7 7]

[50 51 52 53 54 55 56 57 58 59]    [7 7 7 7 7 7 7 7 7 7]

# basic structure
a = np.arange(60).reshape((6, 10))
a_rc = np.zeros_like(a)
bins = [0, 5, 10, 15, 20, 25, 30, 60, 100]
new_bins = [1, 2, 3, 4, 5, 6, 7, 8]
new_classes = zip(bins[:-1], bins[1:], new_bins)
for rc in new_classes:
    q1 = (a >= rc[0])
    q2 = (a < rc[1])
    z = np.where(q1 & q2, rc[2], 0)
    a_rc = a_rc + z
return a_rc
# result returned

 

Lines 2, 3, 4 and 5 describe the array/raster, the classes that are to be used in reclassifying the raster and the new classes to assign to each class.  Line 5 simply zips the bins and new_bins into a new_classes arrangement which will subsequently be used to query the array, locate the appropriate values and perform the assignment (lines 6-10 )

 

Line 3 is simply the array that the results will be placed.  The np.zeros_like function essentially creates an array with the same structure and data type as the input array.  There are other options that could be used to create containment or result arrays, but reclassification is going to be a simple addition process...

 

  • locate the old classes
  • reclass those cells to a new value
  • add the results to the containment raster/array

 

Simple but effective... just ensure that your new classes are inclusive by adding one class value outside the possible range of the data.

 

Line 10 contains the np.where statement which cleverly allows you to put in a query and assign an output value where the condition is met and where it is not met.  You could be foolish and try to build the big wonking query that handles everything in one line... but you would soon forget when you revisit the resultant the next day.  So to alleviate this possibility, the little tiny for loop does the reclassification one grouping at a time and adds the resultant to the destination array.  When the process is complete, the final array is returned.

 

Now on to the arrays/rasters that have nodata values.  The assignment of nodata values is handled by RasterToNumPyArray so you should be aware of what is assigned to it.

 

Raster with nodata values

old                                  new

[--  1  2  3  4  5  6 --  8  9]      [--  1  1  1  1  2  2 --  2  2]

[10 11 12 13 -- 15 16 17 18 19]      [ 3  3  3  3 --  4  4  4  4  4]

[20 -- 22 23 24 25 26 27 -- 29]      [ 5 --  5  5  5  6  6  6 --  6]

[30 31 32 33 34 -- 36 37 38 39]      [ 7  7  7  7  7 --  7  7  7  7]

[40 41 -- 43 44 45 46 47 48 --]      [ 7  7 --  7  7  7  7  7  7 --]

[50 51 52 53 54 55 -- 57 58 59]]     [ 7  7  7  7  7  7 --  7  7  7]

Make a mask (aka ... nodata values) where the numbers are divisible by 7 and the remainder is 0.

Perform the reclassification using the previous conditions.

 

# mask the values
a_mask = np.ma.masked_where(a%7==0, a)
a_mask.set_fill_value(-1)
# set the nodata value

 

The attached sample script prints out the test with the following information:

 

Input array ... type ndarray

...snip...

Reclassification using

:  from [0, 5, 10, 15, 20, 25, 30, 60, 100]

:  to   [1, 2, 3, 4, 5, 6, 7, 8]

:  mask is False value is None

Reclassed array

...snip...

 

Input array ... type MaskedArray

...snip

Reclassification using

:  from [0, 5, 10, 15, 20, 25, 30, 60, 100]

:  to   [1, 2, 3, 4, 5, 6, 7, 8]

:  mask is True value is -1

Reclassed array

...snip....

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

That's about all for now.  Check the documentation on masked arrays and their functions.  Most functions and properties that apply to ndarrays also apply to masked arrays... it's like learning a new language just by altering the pronounciation of what you already know.

Anaconda ... python ... spyder ... conda .... arcgis pro

Ok... change is inevitable... at least we still have IDEs, having harkened from the card reader days.

So this is a running visual of installing Spyder so I can use with ArcGIS PRO and for general python 3.4/5/6 work

I have used pythonwin and pyscripter for some time.

Some people are charmed by pycharm and so people can't wait to be IDLE. ... but for now, I will document Spyder.

I will add to this as I find more..

 

NOTE See the Comments section for other sample images of Spyder

 

Updates:

 2017-10-31 Anaconda 5.0 is out

see the changelog for details

2017-02-18 Anaconda 4.3 is out  

see the full changelog for details...  some highlights...

 

  2016-10-20 Conda and ArcGIS Pro | ArcPy Café 

 

Things just got a whole load easier....

 

Current distribution supports up to and including python 3.6, and a nice new Spyder plus many more upgraded packages.  I am using this for future-proofing my work.  Arc* will eventually be there so you may as well test while you can.

 

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

I will leave the material below as a legacy record since much of it still applies

The original link on ArcGIS Pro and the changes to managing python environments can be found here

.... Python and ArcGIS Pro 1.3 : Conda

 

Related python information can also be found in

.....   The ...py... links

        Coming to Python.... preparation and anticipation

        Python 3.5 comes to iThings

 

Additions and modificationsDocumentation

-

-  2016-07-15  importing xlsxwriter

-  2016-07-15  initial post

  1. Anaconda | Continuum Analytics: Documentation
  2. Get started — Conda documentation
  3. Anaconda package list
  4. Excel plug-ins for Anaconda |
  5. Spyder Documentation

 

:------------------------------------------------------------------------------------------------------------------------:

State the obvious..... Install ArcGIS PRO 

Just follow the instructions.  Don't try and monkey-patch an old machine that barely runs ArcMap.  Start fresh.

 

1.  Install your Python IDE

I have never used Arc*'s built in IDE.  I am not sure why they include it, except for the fact that they can control its installation and no one needs to worry about problems associated with a separate  IDE.  I installed spyder, Jupyter Qt-console and PythonWin.  Spyder and Jupyter are included in the Anaconda distribution that esri provides.  If you don't use one now, then install spyder.

 

2.  Setting your default editor

 If you want to use another one, go to the Project pane, the Geoprocessing Options and do some setup.  Spyder is located in a somewhat cryptic folder path, which I have show in navigation mode and in step 2. as a visual with cutouts for the fiddle bits.  In reality, the portion of the path C:\Program Files\ArcGIS\Pro is the default installation path and the one most likely to be offered during installation and for multi-user machines.  I personally install mine where I know what it is in the folder.... C:\ArcPro.

 

spyder_01.png

 

:------------------------------------------------------------------------------------------------------------------------:

3. The file path to locate the executable

I assume a default path installation in the following.  Everything prior to the ...\bin folder is specific to your machine install path.

spyder_02.png

:------------------------------------------------------------------------------------------------------------------------:

4. The site-package folder

 

What is included without by default from an ArcGIS Pro installation... this is not a complete list of available packages... the list of those, is given above in the table.  The packages come by python version.  We are currently using python 3.4.x in ArcGIS PRO and ArcMap

spyder_03.png

 

:------------------------------------------------------------------------------------------------------------------------:

5. The Spyder folder contents

 

What is in the spyder folder... scientific_startup does some standard imports for what-it-says-type work

 

spyder_04.png

:------------------------------------------------------------------------------------------------------------------------:

6. The pkgs folder

 

A record of the packages that were installed.

 

spyder_05.png

 

:------------------------------------------------------------------------------------------------------------------------:

7.  Importing packages... xlsxwriter demo

 

Here is the workflow I used to import the xlsxwriter module for use in python and arcmap (presumably).

Here is the workflow I used to import the xlsxwriter module for use in python and arcmap (presumably).

From the start button (windows 10, bottom left) navigate to the ArcGIS folder via All Apps find the Python Command Prompt and right-click on it and Run as Administrator

xlsxwriter1.png

Do the conda install xlswriter entry as suggested in the originating post.

xlsxwriter2.png

Hit Enter and away you go. The magic happens and it should be installed.

xlsxwriter3.png

 

At this stage, I went back to Spyder and from the IPython console I tested... looks good

 

xlsxwriter4.png

:------------------------------------------------------------------------------------------------------------------------:

 

More as I find it...

Sliding/Moving windows ...

 

This is the companion to block functions introduced earlier.  I will keep it simple.  The stats functions for rasters with and without nodata values still apply to this type of treatment.  The only difference is how the sub-arrays are generated.  I won't get into the details of how they are done here, but the stride function uses some pretty tricky handling of the input array structure to parse out the moving windows.

 

Remember, block functions move one window size across and down, per movement cycle.  Sliding functions move one cell across and down and sample the cell based on the window size...hence, there are many more windows to sample in sliding function than in block functions.  I will just show the code snippets and you can play on your own.  I have posted sample timing results for comparative purposes for both implementations.

 

 

Code section and comments
  • The stride_tricks does the work.  You can examine the code in the ...site-packages, numpy lib folder.
  • functools provides decorator capabilities.  I tend to import it by default since I use decorators a lot for timing and documentation functions. 
  • The textwrap module provide several functions, with dedent being useful for output documentation. 
  • The numpy set_printoptions function allows you to control the appearance of output.  The version here is quite simplified...refer to the documentation for more information.

 

Stride and block functions provide the means to sample the input array according to the type of data, the size and layout.  You must read the documentation and experiment in order to fully understand.  It really messes with your head and makes the Rubic's cube look simple. The decorator and timer functions I have talked about before.  Refer to my previous blogs for more information.

imports

import numpy as np
from numpy.lib.stride_tricks import as_strided
from functools import wraps
from textwrap import dedent
np.set_printoptions(edgeitems=3,linewidth=80,precision=2,
                    suppress=True,threshold=100)

 

slide (stride) and block functions

#---- functions ----
def _check(a, r_c):
    """Performs the array checks necessary for stride and block.
    : a   - Array or list.
    : r_c - tuple/list/array of rows x cols. 
    :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
    a = np.atleast_2d(a)
    shp = a.shape
    r, c = r_c = ( min(r, a.shape[0]), min(c, shp[1]) )
    a = np.ascontiguousarray(a)
    return a, shp, r, c, tuple(r_c)
   
def stride(a, r_c=(3, 3)):
    """Provide a 2D sliding/moving view of an array. 
    :  There is no edge correction for outputs.
    :
    :Requires
    :--------
    : 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, shp, 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


def block(a, r_c=(3, 3)):
    """See _check and/or stride for documentation.  This function
    :  moves in increments of the block size, rather than sliding
    :  by one row and column
    :
    """

    a, shp, r, c, r_c = _check(a, r_c)
    shape = (a.shape[0]/r, a.shape[1]/c) + r_c
    strides = (r*a.strides[0], c*a.strides[1]) + a.strides
    a_b = as_strided(a, shape=shape, strides=strides).squeeze()
    return a_b

 

timing decorator, decorated stat function and time test

def delta_time(func):
    """simple timing decorator function"""
    import time
    @wraps(func)
    def wrapper(*args, **kwargs):
        print("Timing function for... {}".format(func.__name__))
        t0 = time.perf_counter()        # start time
        result = func(*args, **kwargs)  # ... run the function ...
        t1 = time.perf_counter()        # end time
        print("Results for... {}".format(func.__name__))
        print("  time taken ...{:12.9e} sec.".format(t1-t0))
        #print("\n {}".format(result))  # print within wrapper
        return result                   # return result
    return wrapperdelta_time
def array_mean(a):
    """change the func"""
    a_mean = a.mean(axis=(2, 3))
    return None
def time_test():
    """time test for block and sliding raster"""
    N = [100, 500, 1000, 2000, 3000, 4000, 5000]
    for n in N:
        a = np.arange(n*n).reshape((n, n))
        #a0 = stride(a, block=(3, 3))  # uncomment this or below
        a0 = block(a, block=(3, 3))
        print("\nArray size... {0}x{0}".format(n))
        a_stat = array_mean(a0)       # time function
    return None

 

main section

if __name__=="__main__":
    """comparison between block and slide view"""
    prn=False
    time_test()

Description: This code does...timing  examples for sliding and block mean

 

Timing results

   sliding mean                block mean

array size   total time    total time (seconds)

1000x1000   3.85e-02      4.65...e-03

2000x2000   1.51e-01      1.75...e-02 

3000x3000   3.41e-01      3.86...e-02 

4000x4000   6.37e-01      7.20...e-02 

5000x5000   9.46e-01      1.04...e-01

 

The slowest is obviously the sliding function since it does 9 samples for ever 1 that block does since I used a 3x3 window...I think...it is the easy stuff that confuses.  Enjoy and study the code  and its origins should you be interested.  I should note that the whole area of sliding implementation has been quite the area of investigation.  The implementation that I presented here is one that I chose to use because I liked it.  It may not be the fastest or the bestest... but I don't care...if it takes less than a sip of coffee, it is good enough for me.

 

That's all for now....

Yes ... for a whole $7 Canadian, you too can do cutting edge programming on your iThing ... iPad 2 or above, the iPhone maybe even the iWatch. See ...  Pythonista for iOS

 

Well, if you don't have an iPhone, iPad or iWhatever to program on, you can keep up with conventional python innovations and update here:

 

Coming to Python.... preparation and anticipation

 

Oh... and os.walk alone, is so yesterday... os.scandir is the new kid providing speedups... but failing python 3.5, you can get it from github thankfully:

 

GitHub - benhoyt/scandir: Better directory iterator and faster os.walk(), now in the Python 3.5 stdlib

 

or just wait until 3.5 comes along.

 

Now ... the iStuff

I have mentioned the program Pythonista (available at the iStore) as a great program for doing programming.  I have been using the 2.7 version for about 2 years now and I am impressed with its IDE.  The new version runs both python 3.5 (yes... 3.5) and python 2.7.  It contains many, many modules by default, including numpy and MatPlotLib to name just a few.  You can ship your cool scripts off to github, OneDrive, all those cloud places.  You can even produce pdfs amongst other things.  A downside... some bizarre Apple policy  (no clue what it means) makes it just a tad hard to load a script directly... (but we all know how to copy and paste content don't we, or consult a 12 year old).

 

So, have a gander at your site-packages folder in your python distribution and just add more.  Do you like line numbering? code completion? indent/dedent? cool IDE themes (like the Dark Theme) and whole load of other stuff, then you might be interested.   SciPy and Pandas isn't included yet, so you will have to do your Voronoi Diagrams and Delauney Triangulations the old fashioned way through pure python and/or numpy or MatPlotlib.  Ooops, I forgot, most of you don't have it since you aren't using ArcGIS Pro yet, so you won't miss it (you will have to install Pro to get python 3.4.x and the SciPy stack with other goodies).

 

But how cool... standing at the bus stop... whip out the iPad and run a quick triangulation in Matplotlib... how to impress potential employers who might be watching.

 

image.png

I can't wait until esri gets something that allows me to do analysis on the iThings so I can work with arcpy, numpy and all my favorite libs and pys.  Imagine being able to graph data in the field as it is being collected...

image.png

 

 

Two screen grabs just to give you a quick view of other stuff.  I will update when I play with the new version more.

 

image.png

 

 

Sample script, using a partial split screen for output on an iPad.  Imagine how cool this would look on your iPhone

 

image.png

 

Alternate theme, for the moody programmer...

image.png

 

The python 2.7 and 3.5 distributions have their own access points with their own site_packages folder and even a shared one.  Kind of like a version of condas... maybe iCondas...

image.jpeg

 

Later

... Dan