什么是密铺？密铺，什么是密铺，可以简单理解成铺瓷砖。 就是在一个平面上不留缝隙，把平面全部铺上，虽然咋一听上去不是什么大事，但是密铺是一门非常古老的研究，从古希腊就有记载（可能当时也是研究铺地板怎么好看吧），同时密铺在数学界有着严谨的分类和定义。 我们这里说的密铺包括下文的结论均指平面多边形单密铺。 平面多边形单密铺指使用一种凸多边形图形实现密铺，直接说一些基本结论吧：
所以只需要知道有多少种五边形能密铺，就宣告人类发现了所有的多边形单密铺，这个历程跨越百年；直到2015年，华盛顿大学的教授 Casey Mann 及其妻子 Jennifer McLoud，再加上学生 David Von Derau 发现第15种可以平面密铺的五边形。 尽管如此，人们还是不清楚能够密铺的五边形到底还有没有，毕竟这次发现距离上次已经过了30年，没有人能保证下一个30年不会发现第16种。 直达2017年，法国数学家 Michaël Rao 通过计算机穷举最后得出了结论并证明，那就是能够平面密铺的五边形只有这15种。 自此，人类已经完全发现了能够实现平面多边形单密铺的所有图形。 密铺与 GIS与GIS 的关系那么和 GIS 有什么关系呢？ 密铺和 GIS 的关系太多了，我们平时就有很多使用到密铺的地方，比如 创建渔网、格网 等，简单来说就是方形格网，就像下面这样，常见的还有六边形网格。 这些网格作用很大，可用于很多方面，比如栅格数据按指定形状重采样、区域划分等等。 Python 制作密铺形状ArcGIS 可以创建矩形的密铺网格，10.3之后的版本可以创建六边形密铺网格。 QGIS 3.16 版本除了可以创建矩形、六边形外，还可以创建菱形密铺。 但是呢，我还没有发现五边形密铺，所以我使用 Python 写了两种密铺五边形的创建脚本，结合 ArcPy 可以创建矢量文件，效果如下：   
如何创建这种形状的几何呢？我们以第一个形状为例： 首先分析单一图形，上面两条长边是下面短边长的两倍，最上面的角是60°，剩下的都是120°。知道了这些基础信息，就能在坐标系中准确定位了，只有在坐标系中确定了点的位置，才能继而构建线或者面几何。
坐标系使用的是笛卡尔坐标系（就是最常见的数学坐标系），向上为 Y 轴正方向，向右为 X 轴正方向。在这个坐标系下，现在可以轻松的计算出A、B、C、D、E点（示意图1）的坐标 from __future__ import absolute_importfrom __future__ import division from math import sin, radians, pow, sqrt # 短边 length = 300 # 角度 angle01 = 60 # 垂直距离 leng = length * sin(radians(angle01)) # 300*sin60° pta = (oX + length * 2, oY) ptb = (pta[0] + length / 2, oY + leng) ptc = (pta[0], oY + leng * 2) ptd = (oX + length, ptc[1]) pte = (oX + length / 2, oY + leng * 3) 其中 oX、oY 表示初始点的 x、y 坐标值； pta 表示 A点（示意图1）坐标，后面依次对应。 当然不可能每个五边形的点都算一次，找到一个有规律的整体，然后平移再平移就行了。 示意图1中的1、2、3、4、5、6块五边形组成了一个整体，我们可以轻松的获取到这六个五边形各个点的坐标，然后直接平移整体即可完美覆盖所有区域。 所以现在搞清楚平移距离和方向以及平移次数就完成了。 o 点到 M 点（示意图1绿色线），即整体向正 Y 方向平移，o 点到 N 点（示意图1绿色线），即整体向正 X 方向平移；这种平移可以使用 numpy 进行坐标加减，效率很高，所以这两个脚本的运行效率几乎可以媲美 ArcGIS 的原生工具。 好的，基本思路就是这样，下面上完整代码，感兴趣的读者可以仔细看看，有空我会更新之前的样式箱，加入这几种创建五边形密铺的小工具，到时候会发布更新文章： # -*- coding:utf-8 -*-# ------------------------------------------- # Name: tile # Author: Hygnic # Created on: 2021/7/22 11:11 # Version: # Reference: """ Description: Usage: """ # ------------------------------------------- from __future__ import absolute_import from __future__ import division import os import arcpy import numpy as np from math import sin, radians, pow, sqrt """--------------------------------------""" """--------基本方法------""" def check_extent(input): """ 获取输入图层本身的尺寸 :param input: 输入图层地址 :return: 原点，宽，高 """ lyr = arcpy.mapping.Layer(input) lyr_extent = lyr.getExtent() _origin = (lyr_extent.XMin, lyr_extent.YMin) if lyr.isRasterLayer: # arcpy.env.extent = lyr_extent # arcpy.env.mask = lyr pass return _origin, lyr_extent.width, \ lyr_extent.height, lyr_extent.spatialReference def tile_creator(array_obj, featurecalss): """ 将 array_obj 中的点去除然后生成面 :param array_obj: :param featurecalss: 矢量文件 :return: """ _rows = arcpy.da.InsertCursor(featurecalss, "SHAPE@") for _ii in array_obj: _rows.insertRow([_ii]) del _rows """--------基本方法------""" """--------------------------------------""" """--------------------------------------""" """--------基本属性------""" ws = os.path.abspath(os.getcwd()) arcpy.env.workspace = ws arcpy.env.overwriteOutput = True # 坐标原点 # lyr_o, lyr_w, lyr_h, sr=check_extent("data/grid.shp") lyr_o, lyr_w, lyr_h, sr=check_extent("tiff.tif") print "featureclass width:{}".format(lyr_w) print "featureclass height:{}".format(lyr_h) origin = lyr_o oX = origin[0] oY = origin[1] print "origin X:{}".format(oX) print "origin Y:{}".format(oY) # 五边形短边长度 length = 300 # 角度 angle01 = 60 # 垂直距离 leng = length * sin(radians(angle01)) # 300*sin60° """--------基本属性------""" """--------------------------------------""" """--------------------------------------""" """--------基本坐标------""" #一组五边形在坐标四个象限内的坐标 pta = (oX + length * 2, oY) ptb = (pta[0] + length / 2, oY + leng) ptc = (pta[0], oY + leng * 2) ptd = (oX + length, ptc[1]) pte = (oX + length / 2, oY + leng * 3) quadrant1 = [origin, pta, ptb, ptc, ptd, pte] # 二象限 pta2 = (oX - length * 2, pta[1]) ptb2 = (pta2[0] - length / 2, ptb[1]) ptc2 = (pta2[0], ptc[1]) ptd2 = (oX - length, ptd[1]) pte2 = (oX - length / 2, pte[1]) quadrant2 = [origin, pta2, ptb2, ptc2, ptd2, pte2] # 三象限 ptb3 = (ptb2[0], oY - leng) ptc3 = (ptc2[0], oY - leng * 2) ptd3 = (ptd2[0], ptc3[1]) pte3 = (pte2[0], oY - leng * 3) quadrant3 = [origin, pta2, ptb3, ptc3, ptd3, pte3] # 四象限 ptb4 = (pta[0] + length / 2, oY - leng) ptc4 = (pta[0], oY - leng * 2) ptd4 = (oX + length, ptc4[1]) pte4 = (oX + length / 2, oY - leng * 3) quadrant4 = [origin, pta, ptb4, ptc4, ptd4, pte4] pts = [origin, pta, ptb, ptc, ptd, pte, origin, pta2, ptb2, ptc2, ptd2, pte2, origin, pta2, ptb3, ptc3, ptd3, pte3, origin, pta, ptb4, ptc4, ptd4, pte4] """--------基本坐标------""" """--------------------------------------""" """--------------------------------------""" """--------构建要素------""" # 重要参数和属性 cfm = arcpy.CreateFeatureclass_management shpfile = cfm(ws, "PentagonTile", "polygon", spatial_reference=sr) # 左上方向的偏移距离 offset_x = -length * 3/2 offset_y = 5 * leng # 右下方向的偏移距离 offset_x2 = length*4.5 offset_y2 = -leng # y轴方向循环次数 loop_y = int(lyr_h/(6*leng)) array_pt = [] # 用于存放一整列的五边形 for i in xrange(int(loop_y*1.6)): # 向上偏移距离 # new_pts = [(_[0] - length * 3 / 2 * i, _[1] + 5 * leng * i) for _ in pts] new_pts = [(_[0]+offset_x*i, _[1]+offset_y*i) for _ in pts] A = new_pts[:5] B = new_pts[4:6] + [new_pts[11], new_pts[10], new_pts[0]] C = new_pts[6:11] D = new_pts[12:17] E = new_pts[16:18] + [new_pts[-1], new_pts[-2], new_pts[0]] F = new_pts[-6:-1] array_pt.append(A) array_pt.append(B) array_pt.append(C) array_pt.append(D) array_pt.append(E) array_pt.append(F) np_array = np.array(array_pt) print "Len:{}".format(len(np_array)) # 396 print "Size:{}".format(np_array.size) print "Shape:{}".format(np_array.shape) print "Ndim:{}".format(np_array.ndim) loop_x = int(lyr_w/(4*length)) for _ in xrange(int(loop_x*1.4)): tile_creator(np_array, shpfile) np_array = np_array+(offset_x2, offset_y2) """--------构建要素------""" """--------------------------------------""" 结束语
这篇文章的代码有三点价值：
|
|
|
