[Community] Hexagon grid using shapely

[Stephan Hügel]

Apologies for the general question, but I?m wondering if anyone?s got a recipe to generate a grid of hexagon polygons (or any tessellating shape!) for given coordinate bounds, using Shapely.
Essentially I want to emulate the MMQGIS ?Create Grid Layer? function?

Pointers and suggestions gladly accepted.

[tõnis kärdi]

Hi. Not sure if this will do exactly what You need, but maybe try creating
a grid of points and then buffer each point using a different resolution (
http://toblerity.org/shapely/manual.html#object.buffer) than the default 16.

Post by Stephan Hügel

from shapely.geometry import Point
node = lambda x,y: (x, y)
coords = [[node(x, y) for x in range(10)] for y in range(10)]
hexgrid = [Point(*c).buffer(0.1, resolution=2) for c in coords]
import json
json.dumps(a[0].__geo_interface__)

'{"type": "Polygon", "coordinates": [[[0.1, 0.0], [0.07071067811865482,
-0.07071067811865471], [1.6155445744325868e-16, -0.1],
[-0.0707106781186546, -0.07071067811865492], [-0.1,
-3.2310891488651737e-16], [-0.07071067811865504, 0.07071067811865446],
[-4.624589118372729e-16, 0.1], [0.07071067811865436, 0.07071067811865515],
[0.1, 0.0]]]}'

You'd have to work out the proper range for grid definition and figure out
the right buffer distance aswell.

Hope this helps.

All the best,
T?nis

Post by Stephan Hügel
Apologies for the general question, but I?m wondering if anyone?s got a
recipe to generate a grid of hexagon polygons (or any tessellating shape!)
for given coordinate bounds, using Shapely.
Essentially I want to emulate the MMQGIS ?Create Grid Layer? function?
Pointers and suggestions gladly accepted.
--
s

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.gispython.org/pipermail/community/attachments/20141017/6119a9a4/attachment.htm>

[tõnis kärdi]

.. and for whatever reason I was thinking that a hexagon has 8 sides. So
just discard my previous message, hexa can't be achieved with a buffer as
easily. :S

sorry for the mixup.

t;

Post by tõnis kärdi
Hi. Not sure if this will do exactly what You need, but maybe try creating
a grid of points and then buffer each point using a different resolution (
http://toblerity.org/shapely/manual.html#object.buffer) than the default 16.

Post by Stephan Hügel

from shapely.geometry import Point
node = lambda x,y: (x, y)
coords = [[node(x, y) for x in range(10)] for y in range(10)]
hexgrid = [Point(*c).buffer(0.1, resolution=2) for c in coords]
import json
json.dumps(a[0].__geo_interface__)

'{"type": "Polygon", "coordinates": [[[0.1, 0.0], [0.07071067811865482,
-0.07071067811865471], [1.6155445744325868e-16, -0.1],
[-0.0707106781186546, -0.07071067811865492], [-0.1,
-3.2310891488651737e-16], [-0.07071067811865504, 0.07071067811865446],
[-4.624589118372729e-16, 0.1], [0.07071067811865436, 0.07071067811865515],
[0.1, 0.0]]]}'
You'd have to work out the proper range for grid definition and figure out
the right buffer distance aswell.
Hope this helps.
All the best,
T?nis

Post by Stephan Hügel
Apologies for the general question, but I?m wondering if anyone?s got a
recipe to generate a grid of hexagon polygons (or any tessellating shape!)
for given coordinate bounds, using Shapely.
Essentially I want to emulate the MMQGIS ?Create Grid Layer? function?
Pointers and suggestions gladly accepted.
--
s

[Robert Sanson]

Here is some code I wrote a while ago to create a Shape file of
hexagons using Python + OGR:

import math
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
from pyproj import Proj

#input lower left and upper right bounds and side length
sl = 5000.0 #side length
origx = 1772193.50482601
origy = 5495586.50355601
endx = 1892367.49517399
endy = 5617030.49644399

#output format
outpf = "esri"
#outpf = "wkt"

#calculate coordinates of the hexagon points
p = sl * 0.5 #sin(30)
b = sl * math.cos(math.radians(30))
w = b * 2
h = 2 * sl

#offsets for moving along and up rows
xoffset = b
yoffset = 3 * p

#Proj.4 coordinate system
tmpr = Proj('+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000
+y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs')

#Create the output file object
fp = open("hex5km.txt", "w")
if outpf == "esri":
fp.write("Polygon\n")

#set up OGR stuff
driver = ogr.GetDriverByName('ESRI Shapefile')
newShp = driver.CreateDataSource('hex5km.shp')
layer = newShp.CreateLayer('hex1',geom_type=ogr.wkbPolygon)
field = ogr.FieldDefn('hex_id',ogr.OFTInteger)
field.SetWidth(9)
layer.CreateField(field)
out_ref = osr.SpatialReference()
out_ref.ImportFromEPSG(2193)

startx = origx
starty = origy
row = 1
counter = 0

while starty < endy:
if row % 2 == 0:
startx = origx + xoffset
else:
startx = origx
while startx < endx:
p1x = str(startx)
p1y = str(starty + p)
p2x = str(startx)
p2y = str(starty + (3 * p))
p3x = str(startx + b)
p3y = str(starty + h)
p4x = str(startx + w)
p4y = str(starty + (3 * p))
p5x = str(startx + w)
p5y = str(starty + p)
p6x = str(startx + b)
p6y = str(starty)
wkt = "POLYGON(("+p1x+" "+p1y+","+p2x+" "+p2y+","+p3x+"
"+p3y+","+p4x+" "+p4y+","+p5x+" "+p5y+","+p6x+" "+p6y+","+p1x+"
"+p1y+"))"
feature = ogr.Feature(layer.GetLayerDefn())
feature.SetField(0,counter)
geom = ogr.CreateGeometryFromWkt(wkt)
feature.SetGeometryDirectly(geom)
geom.AssignSpatialReference(out_ref)
layer.CreateFeature(feature)
if outpf == "wkt":
fp.write("POLYGON(("+p1x+" "+p1y+","+p2x+" "+p2y+","+p3x+"
"+p3y+","+p4x+" "+p4y+","+p5x+" "+p5y+","+p6x+" "+p6y+","+p1x+"
"+p1y+"))\n")
else:
fp.write(str(counter)+" 0\n")
fp.write("0 "+p1x+" "+p1y+"\n")
fp.write("1 "+p2x+" "+p2y+"\n")
fp.write("2 "+p3x+" "+p3y+"\n")
fp.write("3 "+p4x+" "+p4y+"\n")
fp.write("4 "+p5x+" "+p5y+"\n")
fp.write("5 "+p6x+" "+p6y+"\n")
fp.write("6 "+p1x+" "+p1y+"\n")
counter = counter + 1
startx = startx + w
starty = starty + yoffset
row = row + 1
if outpf == "esri":
fp.write("END\n")
print str(counter) + " Polygons created"
fp.close()
newShp.Destroy()

Regards,

Robert
Apologies for the general question, but I*m wondering if anyone*s
got a recipe to generate a grid of hexagon polygons (or any tessellating
shape!) for given coordinate bounds, using Shapely.
Essentially I want to emulate the MMQGIS *Create Grid Layer*
function*

Pointers and suggestions gladly accepted.
--
s

_______________________________________________
Community mailing list
***@lists.gispython.org
http://lists.gispython.org/mailman/listinfo/community


This email and any attachments are confidential and intended solely for the addressee(s). If you are not the intended recipient, please notify us immediately and then delete this email from your system.

This message has been scanned for Malware and Viruses by Websense Hosted Security.
www.websense.com

[Stephan Hügel]

Robert,
That got me essentially all the way there. Herewith a gist: https://gist.github.com/urschrei/17cf0be92ca90a244a91

This simply calculates a grid of hexagon coordinates given LL and UR x, y values and side length. The returned list can easily generate Shapely polygons, from which the sky is the limit.

Here’s what the result looks like for London, using a 1000m side length:
Loading Image...
--
s

Post by Robert Sanson
Here is some code I wrote a while ago to create a Shape file of
import math
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
from pyproj import Proj
#input lower left and upper right bounds and side length
sl = 5000.0 #side length
origx = 1772193.50482601
origy = 5495586.50355601
endx = 1892367.49517399
endy = 5617030.49644399
#output format
outpf = "esri"
#outpf = "wkt"
#calculate coordinates of the hexagon points
p = sl * 0.5 #sin(30)
b = sl * math.cos(math.radians(30))
w = b * 2
h = 2 * sl
#offsets for moving along and up rows
xoffset = b
yoffset = 3 * p
#Proj.4 coordinate system
tmpr = Proj('+proj=tmerc +lat_0=0 +lon_0=173 +k=0.9996 +x_0=1600000
+y_0=10000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs')
#Create the output file object
fp = open("hex5km.txt", "w")
fp.write("Polygon\n")
#set up OGR stuff
driver = ogr.GetDriverByName('ESRI Shapefile')
newShp = driver.CreateDataSource('hex5km.shp')
layer = newShp.CreateLayer('hex1',geom_type=ogr.wkbPolygon)
field = ogr.FieldDefn('hex_id',ogr.OFTInteger)
field.SetWidth(9)
layer.CreateField(field)
out_ref = osr.SpatialReference()
out_ref.ImportFromEPSG(2193)
startx = origx
starty = origy
row = 1
counter = 0
startx = origx + xoffset
startx = origx
p1x = str(startx)
p1y = str(starty + p)
p2x = str(startx)
p2y = str(starty + (3 * p))
p3x = str(startx + b)
p3y = str(starty + h)
p4x = str(startx + w)
p4y = str(starty + (3 * p))
p5x = str(startx + w)
p5y = str(starty + p)
p6x = str(startx + b)
p6y = str(starty)
wkt = "POLYGON(("+p1x+" "+p1y+","+p2x+" "+p2y+","+p3x+"
"+p3y+","+p4x+" "+p4y+","+p5x+" "+p5y+","+p6x+" "+p6y+","+p1x+"
"+p1y+"))"
feature = ogr.Feature(layer.GetLayerDefn())
feature.SetField(0,counter)
geom = ogr.CreateGeometryFromWkt(wkt)
feature.SetGeometryDirectly(geom)
geom.AssignSpatialReference(out_ref)
layer.CreateFeature(feature)
fp.write("POLYGON(("+p1x+" "+p1y+","+p2x+" "+p2y+","+p3x+"
"+p3y+","+p4x+" "+p4y+","+p5x+" "+p5y+","+p6x+" "+p6y+","+p1x+"
"+p1y+"))\n")
fp.write(str(counter)+" 0\n")
fp.write("0 "+p1x+" "+p1y+"\n")
fp.write("1 "+p2x+" "+p2y+"\n")
fp.write("2 "+p3x+" "+p3y+"\n")
fp.write("3 "+p4x+" "+p4y+"\n")
fp.write("4 "+p5x+" "+p5y+"\n")
fp.write("5 "+p6x+" "+p6y+"\n")
fp.write("6 "+p1x+" "+p1y+"\n")
counter = counter + 1
startx = startx + w
starty = starty + yoffset
row = row + 1
fp.write("END\n")
print str(counter) + " Polygons created"
fp.close()
newShp.Destroy()
Regards,
Robert
Apologies for the general question, but I*m wondering if anyone*s
got a recipe to generate a grid of hexagon polygons (or any tessellating
shape!) for given coordinate bounds, using Shapely.
Essentially I want to emulate the MMQGIS *Create Grid Layer*
function*
Pointers and suggestions gladly accepted.
--
s
_______________________________________________
Community mailing list
http://lists.gispython.org/mailman/listinfo/community
This email and any attachments are confidential and intended solely for the addressee(s). If you are not the intended recipient, please notify us immediately and then delete this email from your system.
This message has been scanned for Malware and Viruses by Websense Hosted Security.
www.websense.com
_______________________________________________
Community mailing list
http://lists.gispython.org/mailman/listinfo/community