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
>>> Stephan Hügel<ursch...@gmail.com> 18/10/2014 3:33 a.m. >>>
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
Community@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
_______________________________________________
Community mailing list
Community@lists.gispython.org
http://lists.gispython.org/mailman/listinfo/community
