Skip navigation
All People > Dan_Patterson > Py... blog > 2020 > January
2020

Buffering

Probably one of the first things you did in a GIS class.

Select all the fast food outlets that are within 1 mile/km of a school. 

To the rescue... a little buffer, a selectorama, intersect, spatial join... whatever.  At some stage you have used buffers for spatial delineation either by creating one as a feature layer, or virtually when you select things within-a-distance-of.

 

So initially I thought I would just show the geometry created by buffering a simple shape or two to show the point densification on the outward corners.  In the old days a buffer around a point was represented by a 36-sided circle, if you like, an n-gon.  Pretty good if the buffer was small, but seriously lame when the buffer radius was very large.  Crrrrankout a buffer on a shapefile to see what I mean.  Many a line of code was written to approximate a true circle.  Now you can.

 

But what about a buffer?  The outward corners have that circular appearance, maintaining the offset radius (aka the buffer size) from the line work.  One would expect teeny-tiny circular arcs to appear in the geometry representation.  At worse, perhaps an n-gon representation of the arc.

 

Not so.  Here is a square green polygon being shown in Spyder with the arcpy's new svg display.  Pretty cool, but I like my numpy array version better (in red)

 

Geometry representationGeometry coordinates

p_0  # ---- arcpy.Polygon

 

 

Pretty standard, lots of extra information in the

geometry, but when you get to it, it is the 

coordinates to the right we are interested in.

 

To make things a bit easier on the eyes, I

subtracted the origin of the view space from the

x, y values.

p_0[0]
<Array
[<Point (300010.0, 5000000.0, #, #)>,
  <Point (300010.0, 5000010.0, #, #)>,
  <Point (300020.0, 5000010.0, #, #)>,
  <Point (300020.0, 5000000.0, #, #)>,
  <Point (300010.0, 5000000.0, #, #)>]>

Shift the coordinates to the origin

 

coords = [(i.X - 300000, i.Y - 5000000)
           for i in p_0[0]]

coords
[(10.0, 0.0),
(10.0, 10.0),
(20.0, 10.0),
(20.0, 0.0),
(10.0, 0.0)]

p0  # ---- a numpy Geo array of the above

p0
Geo([[10., 0.],
     [10., 10.],
     [20., 10.],
     [20., 0.],
     [10., 0.]])

 

In both cases, the coordinates are held in an array of some kind as in the above. 

 

In the example p_0 is the arcpy.Polygon and p_0[0] is the slice of the first object within... an arcpy.Array of arcpy.Point objects.

 

So what do the buffer coordinates look like?

Buffer coordinatesCommentary

arcpy.Polygon buffer coordinates

b_0 = p_0.buffer(1)
b_0[0] # get the array
<Array
[<Point (300020.0449000001, 4999999.001, #, #)>,
<Point (300020.0, 4999999.0, #, #)>,
<Point (300010.0449000001, 4999999.000499999, #, #)>,
<Point (300010.0, 4999999.000600001, #, #)>,
<Point (300009.0, 5000000.0, #, #)>,
<Point (300009.0, 5000010.0, #, #)>,
<Point (300010.0, 5000011.0, #, #)>,
<Point (300020.0, 5000011.0, #, #)>,
<Point (300021.0, 5000010.0, #, #)>,
<Point (300021.0, 5000000.0, #, #)>,
<Point (300020.0449000001, 4999999.001, #, #)>]>

Where are all the points?  I bet you 

can recognize some, but there is no

n-gon representation.

arcpy's __geo_interface__

b_0.__geo_interface__
{'type': 'MultiPolygon',
'coordinates':
[[[(300020.0449000001, 4999999.001),
    (300020.0, 4999999.0),
     (300010.0449000001, 4999999.000499999),
     (300010.0, 4999999.000600001),
     (300009.0, 5000000.0),
     (300009.0, 5000010.0),
     (300010.0, 5000011.0),
     (300020.0, 5000011.0),
     (300021.0, 5000010.0),
     (300021.0, 5000000.0),
     (300020.0449000001, 4999999.001)]]]}

 

Well! No trace of those extra points.

Just a few lame ones with a teeny

offset from the ones making up the

offset lines.

 

What gives?  There is nothing

visible in the documentation of by examining the existing methods or

properties by conventional methods.

b_0.JSON
'{"curveRings":
  [[[300020.0449000001,4999999.0010000002],
    [300020,4999999],
    [300010.0449000001,4999999.0004999992],
    [300010,4999999.0006000008],
{"c":[ [300009,5000000],
      [300009.29277758137,4999999.2929531736]]},
      [300009,5000010],
{"c":[[300010,5000011],
       [300009.29287232383,5000010.7071276763]]},
      [300020,5000011],
{"c":[[300021,5000010],
       [300020.70712767617,5000010.7071276763]]},
      [300021,5000000],
{"c":[[300020.0449000001,4999999.0010000002],
      [300020.72280301224,4999999.3089512894]]}]],
"spatialReference":{"wkid":2146,"latestWkid":2951}}'

JSON to the rescue?  Everyone loves those curly bracket things.

 

A bit more information like the curve ring thing.  But still no extra groups of points.

b_0svg = b_0.__getSVG__()

b_0svg
'<path fill-rule="evenodd"
fill="#66cc99" stroke="#555555"
stroke-width="2.0" opacity="0.6"
d="
M 300020.0449000001,4999999.001
L 300020,4999999
  L 300010.0449000001,4999999.000499999
... snip
L 300009,5000000
L 300009,5000010
L 300009.0048088113,5000010.098019366
... snip
L 300010,5000011
L 300020,5000011
L 300020.09801936575,5000010.9951911885
...snip
L 300021,5000010
L 300021,5000000
L 300020.9954546047,4999999.904777056
... snip
L 300020.0449000001,4999999.001 z" />'


# ---- now take the same corner as an array

svg_arr
print(arc)                   ID    degrees
[[20.        , 11.        ],   0
[20.09801937, 10.99519119],   1
[20.19509527, 10.98079709],   2    78.75
[20.29029264, 10.95695635],   3
[20.38269452, 10.92389861],   4    67.5
[20.47141085, 10.8819423 ],   5
[20.55558711, 10.83149154],   6    56.3
[20.63441247, 10.7730323 ],   7
[20.70712768, 10.70712768],   8    45.0
[20.7730323 , 10.63441247],   9
[20.83149154, 10.55558711],  10
[20.8819423 , 10.47141085],  11
[20.92389861, 10.38269452],  12
[20.95695635, 10.29029264],  13
[20.98079709, 10.19509527],  14    11.25
[20.99519119, 10.09801937],  15
[21.        , 10.        ]]) 16

SVG to the rescue!

 

Huge swaths of coordinates for each

corner.  

 

Miraculously, they have appeared by

some hidden magic that we are not made party to

 

So, What do the corners look like

As expected.  I created the rounded corners for the numpy-based Geo array.  In the example below, I plotted the SVG corner points and the Geo array points.  An angle of about 5 degrees (give or take) is used by both.  Smooth enough and you can account for area calculations if you know number of outward buffer corners, the total angular passage and the number n-gon shape.  Or just use arcpy's Polygon methods to get the area.

To leave you with some final thoughts. 

You can buffer 

  • with a uniform distance around the shape
  • an offset buffer, where the shape is expanded parallel to the forming lines (see above)
  • a chordal buffer, (Buffer Geometry figure)
  • buffer by increasing the area or perimeter (for polygons obviously).

 

This is an example of the last which is based on the offset buffer, but using area as the parameter.

So remember what goes into the buffer and that not all information that you would expect to see is shown in all forms of geometry representation.

Geometry ...

 

Previously...

Arcpy shapes... viewing in Spyder 

 

I really think it is a bit of overkill to create a FeatureClass just to see some geometry I have created or changed.

As of ArcGIS Pro 2.5, you can view a geometry object inside of Jupyter notebooks, Jupyter-lab or any other thing that supports SVG.  This missive shows how to view geometry from arcpy and also how to deal with geometry without having to go to arcpy to do it.

 

If you don't create or edit geometry with python, you can go now.

 

                                        

Start with geometry

 

Begin with some polygons.

In [ 1]: polys
In [or]: print(polys)
In [or]: polys.__repr__()

Out[1]:
[<Polygon object at 0x25941289d30[0x2594152f288]>,<Polygon object at 0x25941289cf8[0x259414a4bc0]>,
<Polygon object at 0x25941289c88[0x259414a4b98]>,<Polygon object at 0x25941289cc0[0x25940a241c0]>,
<Polygon object at 0x25941289c18[0x25940a240d0]>]

 

As In [ 1]: shows, you can get a 'representation' of the python geometry from within python.  It tells you that it is a polygon and then gives you the memory stuff (I presume).  The former you probably already knew and the latter you probably don't care about.  Pretty useless. So the quest continues.

 

Every object in python has a string representation, so lets try there.

 

In [ 2]: str(p_0)
In [or]: p_0.__str__()

Out[2]:
'<geoprocessing describe geometry object object at 0x000002594152F288>'

 

Heart be still! In [ 2]: is even more cryptic, but the memory location idea was spot on.

Let's see what 'dir' reveals.  I have snipped out a lot of stuff, and highlighted the more useful formats.  Make sure you explore the various geometry classes on your own.


In [3]: dir(p_0)

Out[3]:
['JSON', 'WKB', 'WKT', ... snip ...'__geo_interface__',
'__getSVG__', ... the new addition ....snip ...
'_fromGeoJson', ... snip ... '_repr_svg_', ... snip ...
...all the other properties and methods
]

 

Interesting.. maybe if a geometry contains more than one thing, slicing might real more.

 

In [4]: p_0[0]

Out[4]:
<Array
    [<Point (300010.0, 5000010.0, #, #)>, <Point (300010.0, 5000000.0, #, #)>,
    ... snip ...
    None,
          ... snip ...
    <Point (300002.0, 5000008.0, #, #)>, <Point (300001.0, 5000009.0, #, #)>,
    <Point (300001.0, 5000008.0, #, #)>, <Point (300002.0, 5000008.0, #, #)>]>

 

Perfect!  Polygons are made up of Array objects which are made up of Point objects.  Parts of arrays are separated by None, so you can have the polygons with multiple parts and/or holes in the parts.

Your education has been confirmed.

 

SVG

What is that __getSVG__ thing?

 

In [5]: p_0.__getSVG__()

Out[5]:
'<path fill-rule="evenodd"
fill="#66cc99"
stroke="#555555"
stroke-width="2.0"
opacity="0.6"
d=" M 300010,5000010 L 300010,5000000 L 300001.5,5000001.5 L 300000,5000010 L 300010,5000010
    M 300003,5000009 L 300003,5000003 L 300009,5000003 L 300009,5000009 L 300003,5000009
    M 300002,5000007 L 300001,5000007 L 300002,5000005 L 300002,5000007
    M 300002,5000008 L 300001,5000009 L 300001,5000008 L 300002,5000008
z"
/>'

 

Well they sure look like coordinates.  Time to go off and research SVG construction.

since __repr__ was revealing perhaps _repr_svg_ will be too.

 

In [6]: p_0._repr_svg_()

Out[6]:
'<svg xmlns="http://www.w3.org/2000/svg"
xmlns:xlink="http://www.w3.org/1999/xlink"
width="100.0" height="100.0"
viewBox="299999.6 4999999.6 10.800000000046566 10.800000000745058"
preserveAspectRatio="xMinYMin meet">
<g transform="matrix(1,0,0,-1,0,10000010.0)">
<path fill-rule="evenodd"
fill="#66cc99"
stroke="#555555"
stroke-width="0.21600000001490116"
opacity="0.6"
d=" M 300010,5000010 L 300010,5000000 L 300001.5,5000001.5 L 300000,5000010 L 300010,5000010
    M 300003,5000009 L 300003,5000003 L 300009,5000003 L 300009,5000009 L 300003,5000009
    M 300002,5000007 L 300001,5000007 L 300002,5000005 L 300002,5000007
    M 300002,5000008 L 300001,5000009 L 300001,5000008 L 300002,5000008
z"
/>
</g></svg>'

 

 

 

                                        

 

On to the Display


 Inside of Spyder, which uses an IPython console,

 

I used my handy function for representing a numpy-based array as geometry... which I was working with. 

 

You get a quick peek at what it looks like. 

The alternative???

  •  use my other handy functions and make a featureclass,
  • crank up Pro and
  • add it to a map

 

Saving the SVG

 

Pretty easy.  Just right-click on it and you can save the SVG as an image or a *.svg file for use in such programs as Word.

 

 

 

 

 

 

 

 

 

 

 

 

 

The Code

 

I put it in a gist at...

 

 svg display for numpy geometries

 

If you are interested in alternatives and/or supplements to the various geometry packages, I have been working on

 

 ... npGeom ... which provides the basis for my

 

... Free Tools ...collection of tools normally restricted to the Advanced or Standard license for ArcGIS Pro.

 

So if you like geometry, between arcpy, python, numpy and the various display environments, you can find a match for the job at hand.

Try your hand with a triangle and a square...

t = np.array([[ 0.00,  0.00], [ 0.50,  1.00], [ 1.00,  0.00], [ 0.00,  0.00]])
r = np.array([[ 0.00,  2.00], [ 0.00,  1.00], [ 1.00,  1.00],
              [ 1.00,  2.00], [ 0.00,  2.00]])