GIS 与密铺与铺瓷砖的紧密关系
前言：铺瓷砖是一门古老的技术，与我们息息相关，其中蕴含的数学道理更与 GIS 紧密结合...什么是密铺？
密铺，什么是密铺，可以简单理解成铺瓷砖。就是在一个平面上不留缝隙，把平面全部铺上，虽然咋一听上去不是什么大事，但是密铺是一门非常古老的研究，从古希腊就有记载（可能当时也是研究铺地板怎么好看吧），同时密铺在数学界有着严谨的分类和定义。我们这里说的密铺包括下文的结论均指平面多边形单密铺。平面多边形单密铺指使用一种凸多边形图形实现密铺，直接说一些基本结论吧：任何三角形都可以密铺整个平面；任何凸四边形（包括正方形，矩形）都可以密铺整个平面；正五边形不能密铺；只有正六边形能密铺平面；任何凸 n 边形（n>=7)都没有合适的单密铺方式；...所以只需要知道有多少种五边形能密铺，就宣告人类发现了所有的多边形单密铺，这个历程跨越百年；直到2015年，华盛顿大学的教授 Casey Mann 及其妻子 Jennifer McLoud，再加上学生 David Von Derau 发现第15种可以平面密铺的五边形。
发现的第15种五边形密铺尽管如此，人们还是不清楚能够密铺的五边形到底还有没有，毕竟这次发现距离上次已经过了30年，没有人能保证下一个30年不会发现第16种。直达2017年，法国数学家 Michaël Rao 通过计算机穷举最后得出了结论并证明，那就是能够平面密铺的五边形只有这15种。自此，人类已经完全发现了能够实现平面多边形单密铺的所有图形。
全部15种五边形密铺与 GIS
与GIS 的关系那么和 GIS 有什么关系呢？密铺和 GIS 的关系太多了，我们平时就有很多使用到密铺的地方，比如 创建渔网、格网 等，简单来说就是方形格网，就像下面这样，常见的还有六边形网格。这些网格作用很大，可用于很多方面，比如栅格数据按指定形状重采样、区域划分等等。
渔网（格网）
蜂窝六边形Python 制作密铺形状ArcGIS 可以创建矩形的密铺网格，10.3之后的版本可以创建六边形密铺网格。QGIS 3.16 版本除了可以创建矩形、六边形外，还可以创建菱形密铺。但是呢，我还没有发现五边形密铺，所以我使用 Python 写了两种密铺五边形的创建脚本，结合 ArcPy 可以创建矢量文件，效果如下：
第一种自定义的五边形
栅格数据分区统计效果图
第二种自定义的五边形
栅格数据分区统计效果图如何创建这种形状的几何呢？我们以第一个形状为例：首先分析单一图形，上面两条长边是下面短边长的两倍，最上面的角是60°，剩下的都是120°。知道了这些基础信息，就能在坐标系中准确定位了，只有在坐标系中确定了点的位置，才能继而构建线或者面几何。
角&边长
示意图1坐标系使用的是笛卡尔坐标系（就是最常见的数学坐标系），向上为 Y 轴正方向，向右为 X 轴正方向。在这个坐标系下，现在可以轻松的计算出A、B、C、D、E点（示意图1）的坐标from __future__ import absolute_importfrom __future__ import divisionfrom 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_importfrom __future__ import divisionimport osimport arcpyimport numpy as npfrom 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 = lyrpassreturn _origin, lyr_extent.width, \lyr_extent.height, lyr_extent.spatialReferencedef 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 = wsarcpy.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_ooX = 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_managementshpfile = cfm(ws, "PentagonTile", "polygon", spatial_reference=sr)# 左上方向的偏移距离offset_x = -length * 3/2offset_y = 5 * leng# 右下方向的偏移距离offset_x2 = length*4.5offset_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)) # 396print "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)"""--------构建要素------""""""--------------------------------------"""结束语
这篇文章的代码有三点价值：使用 ArcPy 创建几何；使用 numpy 做简单的矩阵运算；实现了五边形的密铺，提供了更多选择，不仅仅局限于矩形等。荟GIS精粹，关注公众号：GIS荟欢迎交流，更多文章请使用搜索