需求：
使用Hilbert 曲线对遥感影像瓦片数据进行编码，获取某个区域的编码值即可
Hilbert 曲线编码方式
思路
大致可以对四个方向的数据进行归类
- 左下
- 左上
- 右上
- 右下
这个也对应着编码的顺序
思考在不同Hilbert深度（阶）情况下的四个区域的取值范围
-
左下
[ 0 , 4 n − 1 ) [0,4^{n-1}) [0,4n−1)
-
左上
[ 4 n − 1 , 2 × 4 n − 1 ) [4^{n-1}, \ 2\ \times \ 4^{n-1}) [4n−1, 2 × 4n−1)
-
右上
[ 2 × 4 n − 1 , 3 × 4 n − 1 ) [2\ \times\ 4^{n-1},3 \times4^{n-1}) [2 × 4n−1,3×4n−1)
-
右下
[ 3 × 4 n − 1 , 4 n ) [3 \ \times \ 4^{n-1}, 4^n) [3 × 4n−1,4n)
因此可以使用递归的方式来进行生成数据
主要问题
左上角和右上角不需要做旋转操作
// 左上角，待判定点不需要做任何调整，直接递归新子区域即可
case (true, false) =>
recursiveHilbertEncode(point, (Point(lat_half, lon_half), Point(bounds._1.lat,bounds._2.lon)), maxLevel - 1 ) + Math.pow(4, maxLevel-2).toInt
// 右上角
case (false, false) =>
recursiveHilbertEncode(point, (Point(lat_half, lon_half), bounds._2), maxLevel - 1) + 2 * Math.pow(4, maxLevel - 2).toInt
左下角和左上角需要考虑旋转问题
可以使用举例子的方式，来找到旋转后的关系
假设空间区域的左下角坐标为（10，20），右上角的坐标为（20，30），中心坐标为(15, 25)
同样以二阶图为例
这里点位置变更主要针对的是待判定的位置属于哪个区域！
Hilbert 编码区域左下区域(左侧为原始点位置，右侧为旋转后的点位置)
（18，22.5）=> （12.5，28）
（17.5，22） => (12，27.5)
Hilbert 编码区域右下区域
（12，27.5）=>（12.5, 28）
（12.5, 22） => （18, 27.5）
左下旋转后坐标生成方式
原始点位置-中心点位置 = （a,b）
结果 = 中心点位置 + （b,a)
例：（18，22.5） - （15，25） = （3，-2.5）
结果 = (15,25 ) + (-2.5, 3) = (12.5, 28)
右下旋转后坐标生成方式
原始点位置-中心点位置 = （a,b）
结果 = 中心点位置 -（b,a)
例： （12， 27.5） - （15，25） = （-3，2.5）
结果 = （15，25） -（2.5，-3） = （12.5，28）
思路捋清之后，可以考虑代码
// 左下角，
case (true, true) =>
recursiveHilbertEncode(Point(lat_half + lon_reverse,lon_half + lat_reverse), (bounds._1, Point(lat_half, lon_half)),maxLevel - 1)
// 右下角
case (false, true) =>
recursiveHilbertEncode(Point(lat_half - lon_reverse,lon_half - lat_reverse), ((Point(lat_half,bounds._1.lon),Point(bounds._2.lat,lon_half))), maxLevel - 1) + 3 * Math.pow(4, maxLevel-2).toInt
需要对待判定点重新生成，但区域的经纬坐标可以不变，我们不关系点坐标怎么变换，我们只是需要最终的编码结果，保证递归正常运行即可
Hilbert 源码生成
input：
- 正方形边界两个对角顶点的坐标
- 最大Hilbert 层级（阶、深度）
- 待判定的区域两点坐标（或者是一个点）
output：
该区域（区域也是计算区域的中心点位置）、点所在的Hilbert 编码值
递归函数
recursiveHilbertEncode
def recursiveHilbertEncode(point: Point, bounds:(Point,Point), maxLevel:Int): Int = {
if (maxLevel == 1)
return 0
// 边界中心点
val lat_half = (bounds._1.lat + bounds._2.lat) / 2 //x
val lon_half = (bounds._1.lon + bounds._2.lon) / 2 //y
// 反转差值
val lat_reverse = point.lat - lat_half
val lon_reverse = point.lon - lon_half
(point.lat < lat_half, point.lon < lon_half) match {
// 右上角 , p1,p2 不需要动，需要动的是边界和
case (false, false) =>
recursiveHilbertEncode(point, (Point(lat_half, lon_half), bounds._2), maxLevel - 1) + 2 * Math.pow(4, maxLevel - 2).toInt
// 右下角 向左旋转
case (false, true) =>
recursiveHilbertEncode(Point(lat_half - lon_reverse,lon_half - lat_reverse), ((Point(lat_half,bounds._1.lon),Point(bounds._2.lat,lon_half))), maxLevel - 1) + 3 * Math.pow(4, maxLevel-2).toInt
//左上
case (true, false) =>
recursiveHilbertEncode(point, (Point(lat_half, lon_half), Point(bounds._1.lat,bounds._2.lon)), maxLevel - 1 ) + Math.pow(4, maxLevel-2).toInt
// 左下角 向右旋转
case (true, true) =>
recursiveHilbertEncode(Point(lat_half + lon_reverse,lon_half + lat_reverse), (bounds._1, Point(lat_half, lon_half)),maxLevel - 1)
}
}
完整代码
测试代码没整理，有两个main函数，一个是零散的编码区域，一个是从0-15区域的编码结果
import HilbertCurve.getTileCode
object HilbertCurve {
case class Point(lat: Double, lon: Double)
// 递归查看Hilbert编码
/**
* 瓦片中心点坐标
* 新bounds边界（其实也是用来定义边界的中心点的）
* maxLevel
* @return
*/
def recursiveHilbertEncode(point: Point, bounds:(Point,Point), maxLevel:Int): Int = {
if (maxLevel == 1)
return 0
// 边界中心点
val lat_half = (bounds._1.lat + bounds._2.lat) / 2 //x
val lon_half = (bounds._1.lon + bounds._2.lon) / 2 //y
//
// println("边界中心点 :(" + lat_half + ", " + lon_half + ")")
//
// println("瓦片中心点：(" + point.lat + ", " + point.lon + ")")
// point 为瓦片的中心点
// 反转差值
val lat_reverse = point.lat - lat_half
val lon_reverse = point.lon - lon_half
(point.lat < lat_half, point.lon < lon_half) match {
// 右上角 , p1,p2 不需要动，需要动的是边界和
case (false, false) =>
recursiveHilbertEncode(point, (Point(lat_half, lon_half), bounds._2), maxLevel - 1) + 2 * Math.pow(4, maxLevel - 2).toInt
//TODO: 问题
// 右下角 向左旋转
// 此处需要也对point 点进行旋转吗？，好像不需要，这里是沿着中轴线旋转180度
case (false, true) =>
recursiveHilbertEncode(Point(lat_half - lon_reverse,lon_half - lat_reverse), ((Point(lat_half,bounds._1.lon),Point(bounds._2.lat,lon_half))), maxLevel - 1) + 3 * Math.pow(4, maxLevel-2).toInt
//左上
case (true, false) =>
recursiveHilbertEncode(point, (Point(lat_half, lon_half), Point(bounds._1.lat,bounds._2.lon)), maxLevel - 1 ) + Math.pow(4, maxLevel-2).toInt
//TODO： 左下角也有问题
// 左下角 向右旋转
case (true, true) =>
recursiveHilbertEncode(Point(lat_half + lon_reverse,lon_half + lat_reverse), (bounds._1, Point(lat_half, lon_half)),maxLevel - 1)
}
}
// 第一步函数：获取瓦片编码
def getTileCode(p1: Point, p2: Point, bounds: (Point, Point), maxLevel: Int): Int = {
// TODO：此处也默认 bounds 第一个点位于左下角
// TODO: 此处默认p1 位于p2的左下角
// 定义中心点 x，y值
val x_half = (p1.lat + p2.lat)/2
val y_half = (p1.lon + p2.lon)/2
// 边界的中心点
val encoder = recursiveHilbertEncode(Point(x_half,y_half), bounds, maxLevel)
encoder
}
/**
* 测试左下角的16个区域
* @param args
*/
def main2(args: Array[String]): Unit = {
val bounds = (Point(10.0, 10.0), Point(20.0, 20.0))
val maxLevel4 = 4
// TODO： 左下角
// 15 Error
// 最左下角的四个，但是此时如果是旋转过，再去查看右上角和左下角，就会有旋转问题
//
val p1 = Point(10,10)
val p2 = Point(11.25,11.25)
val tileCode0 = getTileCode(p1,p2,bounds,maxLevel4)
println(s"TIle Code0 is given region :$tileCode0")
val p3 = Point(10,11.25)
val p4 = Point(11.25,12.5)
val tileCode1 = getTileCode(p3,p4,bounds,maxLevel4)
println(s"TIle Code1 is given region :$tileCode1")
val p5 = Point(11.25,11.25)
val p6 = Point(12.5,12.5)
val tileCode2 = getTileCode(p5,p6,bounds,maxLevel4)
println(s"TIle Code2 is given region :$tileCode2")
val p7 = Point(11.25,10)
val p8 = Point(12.5,11.25)
val tileCode3 = getTileCode(p7, p8, bounds, maxLevel4)
println(s"TIle Code3 is given region :$tileCode3")
println("==================")
val a = Point(12.5,10)
val b = Point(13.75,11.25)
val tileCode4 = getTileCode(a, b, bounds, maxLevel4)
println(s"TIle Code4 is given region :$tileCode4")
val t1 = Point(13.75,10)
val t2 = Point(15.0, 11.25)
val tileCode5 = getTileCode(t1,t2,bounds,maxLevel4)
println(s"Tile Code5 is given region:$tileCode5")
val c = Point(13.75,11.25)
val d = Point(15,12.5)
val tileCode6 = getTileCode(c,d,bounds,maxLevel4)
println(s"TIle Code6 is given region :$tileCode6")
val b3 = Point(12.5,11.25)
val b4 = Point(13.75,12.5)
val tileCode7 = getTileCode(b3,b4,bounds,maxLevel4)
println(s"TIle Code7 is given region :$tileCode7")
println("=================")
val e = Point(12.5 ,12.5)
val f = Point(13.75, 13.75)
val tileCode8 = getTileCode(e,f,bounds,maxLevel4)
println(s"TIle Code8 is given region :$tileCode8")
val g = Point(13.75, 12.5)
val h = Point(15,13.75)
val tileCode9 = getTileCode(g, h,bounds, maxLevel4)
println(s"TIle Code9 is given region :$tileCode9")
val c3 = Point(12.5,13.75)
val c4 = Point(13.75,15)
val tileCode11 = getTileCode(c3,c4, bounds, maxLevel4)
println(s"TIle Code11 is given region :$tileCode11")
// 左下角的左上角的右上角有问题
println("==========")
val c1 = Point(11.25,13.75)
val c2 = Point(12.5,15)
val tileCode12 = getTileCode(c1,c2,bounds,maxLevel4)
println(s"TIle Code12 is given region :$tileCode12")
val d1 = Point(11.25,12.5)
val d2 = Point(12.5,13.75)
val tileCode13 = getTileCode(d1,d2,bounds,maxLevel4)
println(s"TIle Code13 is given region :$tileCode13")
val d3 = Point(10,12.5)
val d4 = Point(11.25,13.75)
val tileCode14 = getTileCode(d3,d4,bounds,maxLevel4)
println(s"TIle Code14 is given region :$tileCode14")
val d5 = Point(10,13.75)
val d6 = Point(11.25,15)
val tileCode15 = getTileCode(d5,d6,bounds,maxLevel4)
println(s"TIle Code15 is given region :$tileCode15")
// TODO； 右下角TODO Error
// 58
// val p74 = Point(15,10)
// val p84 = Point(16.25,11.25)
// val tileCode44 = getTileCode(p74, p84, bounds, maxLevel4)
// println(s"Tile Code for given region: $tileCode44")
}
def main(args: Array[String]): Unit = {
val bounds = (Point(10.0, 10.0), Point(20.0, 20.0))
// 右上角 42
val p14 = Point(18.75, 18.75)
val p24 = Point(20, 20)
val maxLevel4 = 4
val tileCode14 = getTileCode(p14, p24, bounds, maxLevel4)
println(s"Tile Code for given region: $tileCode14")
// 左上角 //21
val p34 = Point(10, 18.75)
val p44 = Point(11.25, 20)
val tileCode24 = getTileCode(p34, p44, bounds, maxLevel4)
println(s"Tile Code for given region: $tileCode24")
// 5
val p54 = Point(13.75, 10)
val p64 = Point(15, 11.25)
val tileCode34 = getTileCode(p54, p64, bounds, maxLevel4)
println(s"Tile Code for given region: $tileCode34")
// TODO 15 Error
val p545 = Point(10, 13.75)
val p645 = Point(11.25, 15)
val tileCode344 = getTileCode(p545, p645, bounds, maxLevel4)
println(s"Tile Code for given region: $tileCode344")
// 右下角
// TODO Error
// 58
val p74 = Point(15,10)
val p84 = Point(16.25,11.25)
val tileCode44 = getTileCode(p74, p84, bounds, maxLevel4)
println(s"Tile Code for given region: $tileCode44")
// 34
val p1 = Point(16.25, 16.25)
val p2 = Point(17.5, 17.5)
val maxLevel = 3
val tileCode = getTileCode(p1, p2, bounds, maxLevel4)
println(s"Tile Code for given region: $tileCode")
// 10
val p3 = Point(13.75, 13.75)
val p4 = Point(15.0, 15.0)
val tileCode2 = getTileCode(p3, p4, bounds, maxLevel4)
println(s"Tile Code for given region: $tileCode2")
// 8
val p5 = Point(12.5, 12.5)
val p6 = Point(13.75, 13.75)
val tileCode3 = getTileCode(p5, p6, bounds, maxLevel4)
println(s"Tile Code for given region: $tileCode3")
// 32
val p7 = Point(15.0, 15.0)
val p8 = Point(16.25, 16.25)
val tileCode4 = getTileCode(p7, p8, bounds, maxLevel4)
println(s"Tile Code for given region: $tileCode4")
}
}
