nav[aria-label="Primary Navigation"] { padding: 0; & ul { list-style: none; width: 100%; display: flex; flex-direction: row; justify-content: start; align-items: start; gap: 30px; padding: 0; & li { margin: 0; } & ul li { list-style: none; } } }

Numpy Snippets # 3 ... Phish_Nyet ... creating sampling grids using numpy and arcpy

DanPatterson_Retired

Numpy Snippets

Updated: 2016-09-09

The purpose of this post is to show how numpy can play nicely with arcpy to produce geometry and perform tasks like shape translation, rotation and scaling which are not readily available in arcpy's available functions.

To pay homage to ArcMap's ... Fish Net ... I propose a very scaled down version aptly named Phish_Nyet. A fuller version will appear when ArcScript 2.0 opens. This is a demo script. All you do is vary the parameters in the demo() function within the script.

Only the rectangular option will be presented here, however, it is possible to make all kinds of sampling grids, including hexagonal ones as shown above.

The following points apply:

Output is to a shapefile.
The output file must have a projected coordinate system. To create grids in Geographic coordinates, use Fishnet in ArcMap. Why? because invariably people create grids in Geographic coordinates, then project them without densifying the grid blocks. This results in shapes which are in error because the curvature associated with lines of latitude are not accounted for.
A corner point is specified as the origin of the net. One can determine which corner to use by specifying a dX and dY with positive and/or negative values. In this fashion, you can get your corner to be the top or bottom, left or right of the output. If you want the middle to be the origin...do the math and get a corner.
The output is controlled by the cell widths (dX and dY), the number of columns and rows an (i.e. the X and Y directions) and a rotation angle, which is positive for clockwise rotation.

Notes:

- grid_array - does the numpy magic of generating the grid shapes. I have tried to be as verbose as possible. In short, I generate a

seed shape and propagate it. I have intentionally kept rotate and output_polygons as visible function so you can see how they work.

A fuller version will surface as I stated when ArcScripts 2.0 appears. Change the parameters in the demo() function and run it.

<SPAN class="string token">"""
Phish_Nyet.py
Author:  Dan.Patterson@carleton.ca

Purpose: Produce a sampling grid with user defined parameters.
"""</SPAN>
<SPAN class="keyword token">import</SPAN> arcpy
<SPAN class="keyword token">import</SPAN> numpy <SPAN class="keyword token">as</SPAN> np

<SPAN class="keyword token">def</SPAN> <SPAN class="token function">demo</SPAN><SPAN class="punctuation token">(</SPAN><SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">:</SPAN>
    <SPAN class="string token">"""Generate the grid using the following parameter"""</SPAN>
    output_shp <SPAN class="operator token">=</SPAN> r<SPAN class="string token">'C:\temp\Phish_Nyet.shp'</SPAN>
    SR <SPAN class="operator token">=</SPAN>  arcpy<SPAN class="punctuation token">.</SPAN>SpatialReference<SPAN class="punctuation token">(</SPAN><SPAN class="number token">2951</SPAN><SPAN class="punctuation token">)</SPAN> <SPAN class="comment token"># u'NAD_1983_CSRS_MTM_9' YOU NEED ONE!!!</SPAN>
    corner <SPAN class="operator token">=</SPAN> <SPAN class="punctuation token">[</SPAN><SPAN class="number token">340000.0</SPAN><SPAN class="punctuation token">,</SPAN> <SPAN class="number token">5022000.0</SPAN><SPAN class="punctuation token">]</SPAN>     <SPAN class="comment token"># corner of grid </SPAN>
    dX <SPAN class="operator token">=</SPAN> <SPAN class="number token">1000.0</SPAN><SPAN class="punctuation token">;</SPAN>  dY <SPAN class="operator token">=</SPAN> <SPAN class="number token">1000.0</SPAN>          <SPAN class="comment token"># X and Y cell widths</SPAN>
    cols <SPAN class="operator token">=</SPAN> <SPAN class="number token">3</SPAN><SPAN class="punctuation token">;</SPAN>  rows <SPAN class="operator token">=</SPAN> <SPAN class="number token">3</SPAN>                <SPAN class="comment token"># columns/rows...grids in X and Y direction</SPAN>
    angle <SPAN class="operator token">=</SPAN> <SPAN class="number token">0</SPAN>                          <SPAN class="comment token"># rotation angle, clockwise +ve</SPAN>
    <SPAN class="comment token"># create the grid </SPAN>
    pnts <SPAN class="operator token">=</SPAN> grid_array<SPAN class="punctuation token">(</SPAN>corner<SPAN class="punctuation token">,</SPAN>dX<SPAN class="punctuation token">,</SPAN>dY<SPAN class="punctuation token">,</SPAN>cols<SPAN class="punctuation token">,</SPAN>rows<SPAN class="punctuation token">,</SPAN>angle<SPAN class="punctuation token">)</SPAN>
    output_polygons<SPAN class="punctuation token">(</SPAN>output_shp<SPAN class="punctuation token">,</SPAN>SR<SPAN class="punctuation token">,</SPAN>pnts<SPAN class="punctuation token">)</SPAN>
    <SPAN class="keyword token">print</SPAN><SPAN class="punctuation token">(</SPAN><SPAN class="string token">'\nPhish_Nyet has created... {}'</SPAN><SPAN class="punctuation token">.</SPAN>format<SPAN class="punctuation token">(</SPAN>output_shp<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">)</SPAN>

<SPAN class="keyword token">def</SPAN> <SPAN class="token function">grid_array</SPAN><SPAN class="punctuation token">(</SPAN>corner<SPAN class="operator token">=</SPAN><SPAN class="punctuation token">[</SPAN><SPAN class="number token">0</SPAN><SPAN class="punctuation token">,</SPAN><SPAN class="number token">0</SPAN><SPAN class="punctuation token">]</SPAN><SPAN class="punctuation token">,</SPAN>dX<SPAN class="operator token">=</SPAN><SPAN class="number token">1</SPAN><SPAN class="punctuation token">,</SPAN>dY<SPAN class="operator token">=</SPAN><SPAN class="number token">1</SPAN><SPAN class="punctuation token">,</SPAN>cols<SPAN class="operator token">=</SPAN><SPAN class="number token">1</SPAN><SPAN class="punctuation token">,</SPAN>rows<SPAN class="operator token">=</SPAN><SPAN class="number token">1</SPAN><SPAN class="punctuation token">,</SPAN>angle<SPAN class="operator token">=</SPAN><SPAN class="number token">0</SPAN><SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">:</SPAN>
    <SPAN class="string token">"""create the array of pnts to pass on to arcpy using numpy magic"""</SPAN>
    X <SPAN class="operator token">=</SPAN> <SPAN class="punctuation token">[</SPAN><SPAN class="number token">0.0</SPAN><SPAN class="punctuation token">,</SPAN><SPAN class="number token">0.0</SPAN><SPAN class="punctuation token">,</SPAN><SPAN class="number token">1.0</SPAN><SPAN class="punctuation token">,</SPAN><SPAN class="number token">1.0</SPAN><SPAN class="punctuation token">,</SPAN><SPAN class="number token">0.0</SPAN><SPAN class="punctuation token">]</SPAN>                  <SPAN class="comment token"># X,Y values for a unit square</SPAN>
    Y <SPAN class="operator token">=</SPAN> <SPAN class="punctuation token">[</SPAN><SPAN class="number token">0.0</SPAN><SPAN class="punctuation token">,</SPAN><SPAN class="number token">1.0</SPAN><SPAN class="punctuation token">,</SPAN><SPAN class="number token">1.0</SPAN><SPAN class="punctuation token">,</SPAN><SPAN class="number token">0.0</SPAN><SPAN class="punctuation token">,</SPAN><SPAN class="number token">0.0</SPAN><SPAN class="punctuation token">]</SPAN>                  <SPAN class="comment token"># </SPAN>
    seed <SPAN class="operator token">=</SPAN> np<SPAN class="punctuation token">.</SPAN>column_stack<SPAN class="punctuation token">(</SPAN><SPAN class="punctuation token">(</SPAN>X<SPAN class="punctuation token">,</SPAN>Y<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">)</SPAN> <SPAN class="operator token">*</SPAN> <SPAN class="punctuation token">[</SPAN>dX<SPAN class="punctuation token">,</SPAN>dY<SPAN class="punctuation token">]</SPAN>    <SPAN class="comment token"># first array corner values scaled</SPAN>
    u <SPAN class="operator token">=</SPAN> <SPAN class="punctuation token">[</SPAN>seed <SPAN class="operator token">+</SPAN> <SPAN class="punctuation token">[</SPAN>j<SPAN class="operator token">*</SPAN>dX<SPAN class="punctuation token">,</SPAN>i<SPAN class="operator token">*</SPAN>dY<SPAN class="punctuation token">]</SPAN> <SPAN class="keyword token">for</SPAN> i <SPAN class="keyword token">in</SPAN> range<SPAN class="punctuation token">(</SPAN><SPAN class="number token">0</SPAN><SPAN class="punctuation token">,</SPAN>rows<SPAN class="punctuation token">)</SPAN> <SPAN class="keyword token">for</SPAN> j <SPAN class="keyword token">in</SPAN> range<SPAN class="punctuation token">(</SPAN><SPAN class="number token">0</SPAN><SPAN class="punctuation token">,</SPAN>cols<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">]</SPAN>
    pnts <SPAN class="operator token">=</SPAN> np<SPAN class="punctuation token">.</SPAN>array<SPAN class="punctuation token">(</SPAN>u<SPAN class="punctuation token">)</SPAN>                         <SPAN class="comment token"># </SPAN>
    x <SPAN class="operator token">=</SPAN> <SPAN class="punctuation token">[</SPAN>rotate<SPAN class="punctuation token">(</SPAN>p<SPAN class="punctuation token">,</SPAN>angle<SPAN class="punctuation token">)</SPAN> <SPAN class="keyword token">for</SPAN> p <SPAN class="keyword token">in</SPAN> pnts<SPAN class="punctuation token">]</SPAN>        <SPAN class="comment token"># rotate the scaled points </SPAN>
    pnts <SPAN class="operator token">=</SPAN> <SPAN class="punctuation token">[</SPAN> p <SPAN class="operator token">+</SPAN> corner <SPAN class="keyword token">for</SPAN> p <SPAN class="keyword token">in</SPAN> x<SPAN class="punctuation token">]</SPAN>            <SPAN class="comment token"># translate them</SPAN>
    <SPAN class="keyword token">return</SPAN> pnts

<SPAN class="keyword token">def</SPAN> <SPAN class="token function">rotate</SPAN><SPAN class="punctuation token">(</SPAN>pnts<SPAN class="punctuation token">,</SPAN>angle<SPAN class="operator token">=</SPAN><SPAN class="number token">0</SPAN><SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">:</SPAN>
    <SPAN class="string token">"""rotate points about the origin in degrees, (+ve for clockwise) """</SPAN>
    angle <SPAN class="operator token">=</SPAN> np<SPAN class="punctuation token">.</SPAN>deg2rad<SPAN class="punctuation token">(</SPAN>angle<SPAN class="punctuation token">)</SPAN>               <SPAN class="comment token"># convert to radians</SPAN>
    s <SPAN class="operator token">=</SPAN> np<SPAN class="punctuation token">.</SPAN>sin<SPAN class="punctuation token">(</SPAN>angle<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">;</SPAN>  c <SPAN class="operator token">=</SPAN> np<SPAN class="punctuation token">.</SPAN>cos<SPAN class="punctuation token">(</SPAN>angle<SPAN class="punctuation token">)</SPAN>   <SPAN class="comment token"># rotation terms</SPAN>
    aff_matrix <SPAN class="operator token">=</SPAN> np<SPAN class="punctuation token">.</SPAN>array<SPAN class="punctuation token">(</SPAN><SPAN class="punctuation token">[</SPAN><SPAN class="punctuation token">[</SPAN>c<SPAN class="punctuation token">,</SPAN> <SPAN class="operator token">-</SPAN>s<SPAN class="punctuation token">]</SPAN><SPAN class="punctuation token">,</SPAN><SPAN class="punctuation token">[</SPAN>s<SPAN class="punctuation token">,</SPAN> c<SPAN class="punctuation token">]</SPAN><SPAN class="punctuation token">]</SPAN><SPAN class="punctuation token">)</SPAN> <SPAN class="comment token"># rotation matrix</SPAN>
    XY_r <SPAN class="operator token">=</SPAN> np<SPAN class="punctuation token">.</SPAN>dot<SPAN class="punctuation token">(</SPAN>pnts<SPAN class="punctuation token">,</SPAN> aff_matrix<SPAN class="punctuation token">)</SPAN>         <SPAN class="comment token"># numpy magic to rotate pnts</SPAN>
    <SPAN class="keyword token">return</SPAN> XY_r

<SPAN class="keyword token">def</SPAN> <SPAN class="token function">output_polygons</SPAN><SPAN class="punctuation token">(</SPAN>output_shp<SPAN class="punctuation token">,</SPAN>SR<SPAN class="punctuation token">,</SPAN>pnts<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">:</SPAN>
    <SPAN class="string token">"""produce the output polygon shapefile"""</SPAN>
    msg <SPAN class="operator token">=</SPAN> <SPAN class="string token">'\nRead the script header... A projected coordinate system required'</SPAN>
    <SPAN class="keyword token">assert</SPAN> <SPAN class="punctuation token">(</SPAN>SR <SPAN class="operator token">!=</SPAN> None<SPAN class="punctuation token">)</SPAN> <SPAN class="operator token">and</SPAN> <SPAN class="punctuation token">(</SPAN>SR<SPAN class="punctuation token">.</SPAN>type<SPAN class="operator token">==</SPAN><SPAN class="string token">'Projected'</SPAN><SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">,</SPAN> msg
    polygons <SPAN class="operator token">=</SPAN> <SPAN class="punctuation token">[</SPAN><SPAN class="punctuation token">]</SPAN>
    <SPAN class="keyword token">for</SPAN> pnt <SPAN class="keyword token">in</SPAN> pnts<SPAN class="punctuation token">:</SPAN>                           <SPAN class="comment token"># create the polygon geometry</SPAN>
        polygons<SPAN class="punctuation token">.</SPAN>append<SPAN class="punctuation token">(</SPAN>arcpy<SPAN class="punctuation token">.</SPAN>Polygon<SPAN class="punctuation token">(</SPAN>arcpy<SPAN class="punctuation token">.</SPAN>Array<SPAN class="punctuation token">(</SPAN><SPAN class="punctuation token">[</SPAN>arcpy<SPAN class="punctuation token">.</SPAN>Point<SPAN class="punctuation token">(</SPAN><SPAN class="operator token">*</SPAN>xy<SPAN class="punctuation token">)</SPAN> <SPAN class="keyword token">for</SPAN> xy <SPAN class="keyword token">in</SPAN> pnt<SPAN class="punctuation token">]</SPAN><SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">,</SPAN>SR<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">)</SPAN>
    <SPAN class="keyword token">if</SPAN> arcpy<SPAN class="punctuation token">.</SPAN>Exists<SPAN class="punctuation token">(</SPAN>output_shp<SPAN class="punctuation token">)</SPAN><SPAN class="punctuation token">:</SPAN>               <SPAN class="comment token"># overwrite any existing versions</SPAN>
        arcpy<SPAN class="punctuation token">.</SPAN>Delete_management<SPAN class="punctuation token">(</SPAN>output_shp<SPAN class="punctuation token">)</SPAN>
    arcpy<SPAN class="punctuation token">.</SPAN>CopyFeatures_management<SPAN class="punctuation token">(</SPAN>polygons<SPAN class="punctuation token">,</SPAN> output_shp<SPAN class="punctuation token">)</SPAN>

<SPAN class="keyword token">if</SPAN> __name__ <SPAN class="operator token">==</SPAN> <SPAN class="string token">'__main__'</SPAN><SPAN class="punctuation token">:</SPAN>
    <SPAN class="string token">"""Generate the grid using the listed parameters"""</SPAN>
    demo<SPAN class="punctuation token">(</SPAN><SPAN class="punctuation token">)</SPAN>  <SPAN class="comment token"># modify the parameters in demo to run</SPAN>
‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍‍<SPAN class="line-numbers-rows"><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN><SPAN>‍</SPAN></SPAN>

Of course other regular geometric shapes can be generated in a similar fashion, but not all may pack like rectangles and others do.

Finde mehr Beiträge mit dem Tag

Py...Blog

Kommentare

DanPatterson_Retired

deltax and deltay

AbdulMajid

sorry for bothering you so much, i just wanted to ask from where i can control the size of hex either in side length or area?

arcmap window code

DanPatterson_Retired

Projected coordinates only... no support for decimal degree data, use esri's tools instead

AbdulMajid

thankyou for your reply,

i tried using your sampling tool, i want the grid to be starting from a specified lat/lon with specified number of rows and columns. i really tried hard in figuring it out.

DanPatterson_Retired

try Sampling Grids for ArcGIS Pro

AbdulMajid

i am trying to make a hex grid with desired number of rows and columns, your post is helpful in generating a shape file in the end. but i cant figure out how to do the hexagon thing. can you let me know how to do that?

DanPatterson_Retired

In a pure numpy solution, affine has been changed so that the matrix for translation, rotation, scale and skew can be applied to a numpy array of 'points' in one foul swoop which does a polygon all at once rather than by point by point. This vectorization relies on the vectorization properties. I should update the code when I get a chance

BrianWilson

Matrix methods are very efficient because you can roll all the operations into one step: translate, rotate, scale.

If you are doing millions of transforms it will help. If you are only doing a few thousand and it's easier for you to understand, you can just code up the transform one step at a time. I did that to generate a small grid but now that I know about numpy I have to look at it again.

I wrote a function like this to do the transform:

def affine(p, scale, theta, offset):
""" Scale, rotate and translate point """
return arcpy.Point(

(p.X * math.cos(theta) - p.Y * math.sin(theta)) * scale.X + offset.X,
(p.X * math.sin(theta) + p.Y * math.cos(theta)) * scale.Y + offset.Y)

(from GitHub - Geo-CEG/shoveltest: Create grids to facilitate shovel tests for archaeology site work.)

p is an arcpy.Point to be transformed, scale is a point defining how much to scale the point, theta is the rotation angle in radians, and offset says how far to move the point.

In practice, I use a "for" loop to march through the points in a feature and call this function each time to move a point from one place to another.I do the same calculations (sine and cosine) each time but I am only moving about 100 points so I don't care, it takes longer for the "arcpy" module to load than to run my program.

Dan's code relies on numpy to do the looping-- in the "rotate" function above, he builds a transformation matrix and then hands an array of points and the matrix to the "dot" method, and all the work takes place there. The numpy library is written in 'c' so it is probably more than 1000 times faster than what I did. Getting the work done is often more important than finding the most efficient code.

In Dan's example he only rotates a point with the matrix. He could have added code to build a matrix that would also scale and translate the points at the same time. Hmm I guess that will have to go into my rewrite...

Hornbydd

I'm not a mathematician so find the matrices bit a bit bamboozling but this is nifty bit of code and I can see how one could incorporate this into a larger workflow, nice!

DanPatterson_Retired

You're welcome Jeff...I am in the process of organizing 'NumPy Repository' group with more stuff in a centralized location for all things pythonic. Keep your eyes open for it...

JeffWard

I've had fun messing around with this. I have successfully created a grid of hexagons. Thanks for indirectly teaching me about list comprehensions.

Nehmen Sie an der Diskussion teil

Angemeldete Mitglieder können Beiträge verfassen, Updates folgen und mehr. Neu hier? Registrieren Sie ein kostenloses Konto.

Schnellzugriffe

Community Help
Find useful guides, FAQs, and documents to help you navigate and make the most of Esri Community.