一、起因
只因为一句话，刺点太累了。刺点的工作除了选点之外，其余十分机械化，重复性很强。 而且人工操作效率低、精度也不一定会高。 这也为电脑辅助刺点提供了可能。 因此就想能不能做一个自动化工具来辅助刺点， 即使是半自动也可以，也会省下很多时间。 基于目前掌握的技术实现了半自动刺点辅助工具。 在下面会详细介绍。
二、整体思路
目前实现的是半自动刺点辅助工具。程序工作流程如下图所示。
首先用户在ENVI中人工选择地物点，然后获取到其对应的经纬度坐标。
然后将经纬度坐标输入程序，程序会自动从Google地图上获取以该点为中心的某一范围内的影像，
同时将一些显著的地物点提取出来。用户在提取出来的特征点中找到对应的点并点击，
便可以自动获取到Google地图上对应的真实经纬度和高程。这样整个刺点过程就完成了。
图中蓝色流程图是用户需要进行的操作，橙色部分是当用户输入经纬度后，程序处理的流程。
三、技术要点
1.由经纬度计算瓦片行列号
由于使用的是标准WMTS服务，如这里使用的是Google WMTS服务，也可以换成其它的如天地图、Bing等。 因此可以很方便地计算出包含某个经纬度的在某一缩放等级下的瓦片行列号。这里就不介绍原理了， 直接给出计算公式。
\[x = \left [ 2^{z-1} \cdot \left ( \frac{\lambda}{180}+1 \right )\right ]\] \[y = \left [ 2^{z-1}\cdot \left ( 1-\frac{ln\left [ tan\left ( \frac{\pi\cdot \varphi }{180} \right )+sec\left ( \frac{\pi\cdot\varphi}{180} \right ) \right ]}{\pi} \right ) \right ]\]其中λ为经度、φ为纬度，单位都为度。z为缩放级数，范围为0-18。
在z级数下，x、y方向都有2^z个瓦片，x、y的取值范围为0-(2^z-1)。
这样便能算出包含该点的在z缩放级数下的瓦片的行列号。
整个瓦片的划分如下图所示：
在Google WMTS服务中卫星地图/地形图返回的是jpg格式，地图返回的是png格式。
2.瓦片拼接
由于考虑到有卫星影像自带的RPC模型算出来的经纬度有误差(其实刺点就是为了验证误差到底有多大)，
所以不能只取包含该点的一个瓦片进行特征点寻找，否则很有可能找不到。
因此在本程序中采用的是以特征点所在的瓦片为中心，构造3×3的瓦片阵列进行寻找，如下图所示。
其中蓝色表示根据用户输入的λ0、φ0计算得到的瓦片(x,y)，
而实际通过用户识别发现其真实对应的地物点应该是图中红色五角星的位置(λ1,φ1)。
所以我们就能得到对于该点影像的精度，其偏离真实值的为(λ1-λ0,φ1-φ0)。
当然考虑到不同地图的偏移等问题，其实可以构造更大的范围进行寻找，如5×5。但在本程序中只是以3×3为例，
它构造好了其它尺寸的也都没有什么问题了，只是数字不同而已。
在获取到了9个瓦片以后，采用Numpy的矩阵操作进行拼接，主要使用了np.hstack()和np.vstack()两个函数，
将它们拼接成一张768×768的影像(每张影像的大小都是256×256)。
3.提取特征点
提取特征点采用的是“冗余策略”，即尽可能把范围内的特征点都提供给用户，然后让用户自己去选择最好的。 因为本程序并没有用到匹配等算法，所以只能通过人工选择来解决这个问题。 提取特征点的操作使用的是OpenCV中的Harris算子。这里就不详细介绍Harris算子的流程和原理了， 以后有时间会专门学习一下。这里只是使用OpenCV的代码调用一下即可。
4.用户交互
由于程序中采取的策略是让用户在提取了特征点的图像上点击，然后输出该特征点的信息。
因此需要能获取用户在图像中点击的坐标的功能，这基于OpenCV实现。
其次考虑到精度问题，用户可能点击的并不是特征点，而是特征点周围，或者点击的不是很准确。
因此并不能直接拿用户点击的坐标来当作特征点。
在程序中采取的策略是，计算用户点击的点与其它所有提取出来的特征点的距离，
然后得到距离最近的那个特征点，我们即认为是用户想要选择的特征点。
这样不管用户是否是在特征点周围或是边缘，都会返回那个点，如下图所示。
红色和蓝色圆点分别代表两个不同的特征点，用户需要的是红色的特征点，
绿色的鼠标不论在哪里点击，只要其到红色点的距离较近，程序都会返回红色特征点的坐标。
由于欧式距离计算相对复杂，因此在程序中采用的是曼哈顿距离(\(dist((x_{1},y_{1}),(x_{2},y_{2}))=\left | x_{1}-x_{2} \right |+\left | y_{1}-y_{2} \right |\))。
5.由瓦片坐标反算经纬度
在上一步获取到了特征点在图像上的像素坐标(x,y)，下一步就需要根据像素坐标反算经纬度。公式如下：
\[\lambda = \left [ 2^{1-z}\cdot \left ( x+\frac{m}{256} \right )-1\right ]\cdot 180\] \[\varphi = \frac{360\cdot arctan\left ( e^{\left [ 1-2^{1-z}\cdot \left ( y+\frac{n}{256} \right ) \right ]\cdot \pi} \right )}{\pi}-90\]其中z表示缩放级数，WMTS中规定取值范围为0-18，18为最大缩放级数。 x、y分别为该瓦片的行列号信息。 m、n表示该点在该瓦片中的像素坐标。m表示横向，n表示竖向。 像素坐标以左上角为原点，x轴向右，y轴向下。 计算出的经纬度单位为度。
6.获取高程
这里使用Google Maps API中的Google Maps Elevation API，更多信息可点击这里。 使用这个API和百度、高德地图等API一样是需要注册获取Key的，不过注册很简单。 通过HTTP接口访问Google Maps Elevation API，以网址字符串形式构建请求，并利用纬度/经度坐标标识位置。 服务器返回json格式的字符串，对其进行解析即可获得对应点的高程信息。如下是一个返回的json字符串示例：
{
"results" : [
{
"elevation" : 1608.637939453125,
"location" : {
"lat" : 39.73915360,
"lng" : -104.98470340
},
"resolution" : 4.771975994110107
}
],
"status" : "OK"
}
四、实现代码
解决了以上几点技术问题后，便可以编写代码实现了。代码如下：
# coding=utf-8
import urllib2 as ulb
import numpy as np
import PIL.ImageFile as ImageFile
import cv2
import math
import json
import win32clipboard as w
import win32con
# 写入剪切板
def setText(aString):
w.OpenClipboard()
w.EmptyClipboard()
w.SetClipboardData(win32con.CF_TEXT, aString)
w.CloseClipboard()
# 获取瓦片影像数据
def getTile(url):
# 用urllib2库链接网络图像
response = ulb.Request(url)
fp = ulb.urlopen(response) # 打开网络图像文件句柄
p = ImageFile.Parser() # 定义图像IO
while 1: # 开始图像读取
s = fp.read(1024)
if not s:
break
p.feed(s)
im = p.close() # 得到图像
arr = np.array(im) # 将图像转换成numpy矩阵
pic = np.zeros(arr.shape, np.uint8)
pic[:, :, 0] = arr[:, :, 2]
pic[:, :, 1] = arr[:, :, 1]
pic[:, :, 2] = arr[:, :, 0]
return pic
# 获取最近特征点
def getNearestPoint(event, x, y, flags, param):
distance = 1000
index = -1
if event == cv2.EVENT_LBUTTONDOWN:
for inx in range(harrisPoints.__len__()):
temp = abs(harrisPoints[inx][0] - x) + abs(harrisPoints[inx][1] - y)
if temp < distance:
distance = temp
index = inx
dx = harrisPoints[index][0] / 256
dy = harrisPoints[index][1] / 256
xp = (X - 1) + dx
yp = (Y - 1) + dy
m = harrisPoints[index][0] - dx * 256
n = harrisPoints[index][1] - dy * 256
latitude, longitude = calcLatLon(xp, yp, 18, m, n)
ele = getElevation(latitude, longitude)
print 'The location of nearest point:(' + latitude.__str__() + "," + longitude.__str__() + ')'
print 'The elevation of nearest point:' + ele.__str__() + " m"
result = latitude.__str__() + "\t" + longitude.__str__() + "\t" + ele.__str__()
setText(result)
print '------copy finished------'
print 'Close the window to continue or press ESC to exit.'
return i
# 计算瓦片行列号
def calcXY(lat, lon):
x = math.floor(math.pow(2, 17) * ((lon / 180.0) + 1))
tan = math.tan(lat * math.pi / 180.0)
sec = 1.0 / math.cos(lat * math.pi / 180.0)
log = math.log(tan + sec)
y = math.floor(math.pow(2, 17) * (1 - log / math.pi))
return int(x), int(y)
# 计算经纬度
def calcLatLon(x, y, z, m, n):
lon = (math.pow(2, 1 - z) * (x + m / 256.0) - 1) * 180.0
lat = (360 * math.atan(math.pow(math.e, (1 - math.pow(2, 1 - z) * (y + n / 256.0)) * math.pi))) / math.pi - 90
return lat, lon
# 获取高程
def getElevation(lat, lon):
url = 'https://maps.googleapis.com/maps/api/elevation/json?locations=' + lat.__str__() + ',' + lon.__str__() + '&key=AIzaSyA1c9KcHycinBjvnmXhEulkORsT1IEF3Dc'
response = ulb.urlopen(url)
response = response.read()
s = json.loads(response)
elevation = float(s['results'][0]['elevation'])
return elevation
# 将字符串度分秒转换为度
def cvtStr2Deg(deg, min, sec):
result = int(deg) + int(min) / 60.0 + float(sec) / 3600.0
return result
# 注意字符编码问题
# 在使用PyInstaller转换后，是不支持'°'这个符号的。因为它不支持utf-8编码，
# 只能使用ascii编码。所以一运行到这里便会报错。但是在Python环境下运行正常，
# 因为指定了utf-8编码。所以如果是在Python环境下可以用这个函数，如果不是，
# 请使用getNum2函数。并且在ENVI里设置经纬度显示方式为十进制。经纬度输入
# 格式是与ENVI保持统一的。格式如下：
# getNum 49°36'53.18"N,45°40'48.05"E
# getNum2 49.61570775N,45.66817893E
def getNum(str):
split = str.split(',')
du = split[0].split('°')[0]
fen = split[0].split('°')[1].split('\'')[0]
miao = split[0].split('°')[1].split('\'')[1].split('"')[0]
split1 = cvtStr2Deg(du, fen, miao)
du = split[1].split('°')[0]
fen = split[1].split('°')[1].split('\'')[0]
miao = split[1].split('°')[1].split('\'')[1].split('"')[0]
split2 = cvtStr2Deg(du, fen, miao)
return split1, split2
def getNum2(str):
split = str.split(',')
split1 = float(split[0].split('N')[0])
split2 = float(split[1].split('E')[0])
return split1, split2
# 获取用户在图像上点击的坐标
def getXY(event, x, y, flags, param):
if event == cv2.EVENT_LBUTTONDOWN:
print x, y
while 1:
# 输入经纬度，ENVI格式
raw = raw_input("Input lat & lon(e.g. 30.52N,114.36E):\n")
# 转换获取经纬度
lat, lon = getNum2(raw)
# 由经纬度计算得到该点所在瓦片行列号（中心点瓦片行列号）
X, Y = calcXY(lat, lon)
# 拼图大小，默认3×3
row = 5
col = 5
# 构建图像用于显示和处理
combine = np.zeros((row * 256, col * 256, 3), np.uint8)
combine2 = np.zeros((row * 256, col * 256, 3), np.uint8)
harrisPoints = []
# 计算索引范围
col_half_range = int((col - 1) / 2)
row_half_range = int((row - 1) / 2)
# 建立相关变量
base = 'https://mt2.google.cn/vt/lyrs=s&hl=zh-CN&gl=CN&x='
urls = []
tiles = []
lines = []
# 循环构造每个瓦片的url地址并放进list
for i in range(-row_half_range, row_half_range + 1, 1):
for j in range(-col_half_range, col_half_range + 1, 1):
url = base + (X + j).__str__() + '&y=' + (Y + i).__str__() + '&z=18'
urls.append(url)
# 对每个item进行获取，并将结果放在tiles中
for item in urls:
tiles.append(getTile(item))
# 先按照水平方向将各瓦片拼起来
for k in range(0, tiles.__len__(), col):
h_tuple = tuple(tiles[k:k + col])
line = np.hstack(h_tuple)
lines.append(line)
# 按照竖直方向将每行拼起来
v_tuple = tuple(lines)
combine = np.vstack(v_tuple)
# 复制新的用于显示
combine2[:, :, 0] = combine[:, :, 0]
combine2[:, :, 1] = combine[:, :, 1]
combine2[:, :, 2] = combine[:, :, 2]
gray = cv2.cvtColor(combine, cv2.COLOR_BGR2GRAY)
gray = np.float32(gray)
dst = cv2.cornerHarris(gray, 2, 9, 0.04)
# Threshold for an optimal value, it may vary depending on the image.
combine[dst > 0.01 * dst.max()] = [0, 0, 255]
# result is dilated for marking the corners, not important
dst = cv2.dilate(dst, None)
combine2[dst > 0.01 * dst.max()] = [0, 0, 255]
for i in range(combine.shape[0]):
for j in range(combine.shape[1]):
if combine[i, j, 0] == 0 and combine[i, j, 1] == 0 and combine[i, j, 2] == 255:
harrisPoints.append((i, j))
print harrisPoints.__len__().__str__() + ' harris points were found in total.'
cv2.namedWindow("pick up window")
cv2.setMouseCallback("pick up window", getNearestPoint)
cv2.imshow("pick up window", combine2)
k = cv2.waitKey(0) & 0xff
if k == 27:
cv2.destroyAllWindows()
break
else:
cv2.destroyAllWindows()
continue
五、测试效果
首先利用ENVI打开遥感影像，并选择一个特征点，如下：
然后将经纬度坐标输入程序中，并按回车，程序会自动找到其对应的区域并提取出特征点。
在控制台中会显示一共找到了多少特征点。
然后点击和我们选择对应的那个特征点，在控制台中会自动输出该特征点相关信息，并且自动将内容复制到粘贴板中。
在需要粘贴的地方粘贴即可。
关闭”pick up window”便可以继续在控制台中输入新点。
利用本程序辅助刺点，可以节省很多时间，也更容易了。
本程序省去了在Google Earth中浏览寻找目标点，再新建Marker获取其坐标、高程，再复制的步骤。
用户只需要在识别出的特征点中选择对应的那个点，程序便会自动返回坐标、高程，以及自动复制信息到粘贴板。
用户需要做的就是把它粘贴到指定地方即可。
最后，将代码打包成了exe，并在平板Win10上进行了测试，如下图。
可以正常运行。打包的exe文件在百度云，点击下载，密码：acfh。
六、问题与改进
1.特征提取
利用Harris算子有时并不能很好提出所有我们想要的特征点。有时会出现程序提取的所有待选点中， 没有我们选择的那个点的情况。这样就没法继续进行下去了。因此还需要继续完善，或改进算法或改变操作方式。
2.地图偏移
在使用过程中发现，在某些区域通过程序获得的候选范围中根本找不到对应特征点，因为偏移很大。 造成这种情况的原因有两个。一个是有可能影像本身的精度很差，和真实地物偏差很大，所以导致跑偏。 第二种情况就是Google地图本身的地图加偏造成的。这两个问题都比较“客观”，不太容易从代码中解决。 目前想到的可能解决办法就是扩大候选区域窗口。但这样可能势必会增加算法的计算量。 所以可能还需要其它更好的办法。
3.选点算法优化
在使用时也发现，有时选出的点并不是想要的点。这就与选点的算法有关。 目前的算法计算量也比较大。如果提取了上千个点，那选一次点就需要进行上千次的距离计算。 应该寻找一种更高效的算法来解决这个问题。
4.瓦片404
在使用过程中发现，获取某些区域瓦片时会报404错误。在检查了拼接的请求url没有问题之后， 初步判断是服务器上没有这片区域18级的影像。这在之前使用天地图的时候就经常遇到。 尤其是在山区等地方，天地图的影像有时只到16级，再放大就没有了。
5.图像分块
为了更快速寻找特征点，可以考虑图像分块，每块限定特征点个数。 如分4块，每块400个点。判断用户点击对应的是哪一块，直接在那一块对应的400个点中进行搜索。 这样可提高效率，解决第三点提到的问题。要不然每次都需要在1600个点中搜索，相对慢一些。 同时这也能保证特征点的分布更加均匀一些，而不是集中在图像的某个区域。
七、瓦片拼接小程序
把上述代码中获取某一范围瓦片的代码单独抽取出来，形成了这样一个小功能。
# coding=utf-8
import urllib2 as ulb
import numpy as np
import PIL.ImageFile as ImageFile
import math
import cv2
# 获取瓦片影像数据
def getTile(url):
# 用urllib2库链接网络图像
response = ulb.Request(url)
fp = ulb.urlopen(response) # 打开网络图像文件句柄
p = ImageFile.Parser() # 定义图像IO
while 1: # 开始图像读取
s = fp.read(1024)
if not s:
break
p.feed(s)
im = p.close() # 得到图像
arr = np.array(im) # 将图像转换成numpy矩阵
pic = np.zeros(arr.shape, np.uint8)
pic[:, :, 0] = arr[:, :, 2]
pic[:, :, 1] = arr[:, :, 1]
pic[:, :, 2] = arr[:, :, 0]
return pic
# 计算瓦片行列号
def calcXY(lat, lon):
x = math.floor(math.pow(2, 17) * ((lon / 180.0) + 1))
tan = math.tan(lat * math.pi / 180.0)
sec = 1.0 / math.cos(lat * math.pi / 180.0)
log = math.log(tan + sec)
y = math.floor(math.pow(2, 17) * (1 - log / math.pi))
return int(x), int(y)
# 计算经纬度
def calcLatLon(x, y, z, m, n):
lon = (math.pow(2, 1 - z) * (x + m / 256.0) - 1) * 180.0
lat = (360 * math.atan(math.pow(math.e, (1 - math.pow(2, 1 - z) * (y + n / 256.0)) * math.pi))) / math.pi - 90
return lat, lon
# 将字符串度分秒转换为度
def cvtStr2Deg(deg, min, sec):
result = int(deg) + int(min) / 60.0 + float(sec) / 3600.0
return result
# 注意字符编码问题
# 在使用PyInstaller转换后，是不支持'°'这个符号的。因为它不支持utf-8编码，
# 只能使用ascii编码。所以一运行到这里便会报错。但是在Python环境下运行正常，
# 因为指定了utf-8编码。所以如果是在Python环境下可以用这个函数，如果不是，
# 请使用getNum2函数。并且在ENVI里设置经纬度显示方式为十进制。经纬度输入
# 格式是与ENVI保持统一的。格式如下：
# getNum 49°36'53.18"N,45°40'48.05"E
# getNum2 49.61570775N,45.66817893E
def getNum(str):
split = str.split(',')
du = split[0].split('°')[0]
fen = split[0].split('°')[1].split('\'')[0]
miao = split[0].split('°')[1].split('\'')[1].split('"')[0]
split1 = cvtStr2Deg(du, fen, miao)
du = split[1].split('°')[0]
fen = split[1].split('°')[1].split('\'')[0]
miao = split[1].split('°')[1].split('\'')[1].split('"')[0]
split2 = cvtStr2Deg(du, fen, miao)
return split1, split2
def getNum2(str):
split = str.split(',')
split1 = float(split[0].split('N')[0])
split2 = float(split[1].split('E')[0])
return split1, split2
# 输入经纬度，ENVI格式
raw = raw_input("Input lat & lon(e.g. 30.52N,114.36E):\n")
# 输入需要获取的瓦片层数，默认为最高层
z = 18
z = input("Level(0-18):\n")
# 拼图大小，默认3×3
row = 3
col = 3
row = input("Window size row(odd number):\n")
col = input("Window size column(odd number):\n")
# 如果用户输入的是偶数，加1变成奇数
if row % 2 == 0:
row += 1
if col % 2 == 0:
col += 1
print "The final image is " + (row * 256).__str__() + " * " + (col * 256).__str__()
# 缩放因子
scale_factor = 1
scale_factor = input("Input the scale factor for image to show(0-1):\n")
# 保存路径
save_path = ""
save_path = raw_input("Input the output path of image(input \'no\' means don't save):\n")
# 转换获取经纬度
lat, lon = getNum2(raw)
# 由经纬度计算得到该点所在瓦片行列号（中心点瓦片行列号）
X, Y = calcXY(lat, lon)
# 计算索引范围
col_half_range = int((col - 1) / 2)
row_half_range = int((row - 1) / 2)
# 基本变量
base = 'https://mt2.google.cn/vt/lyrs=s&hl=zh-CN&gl=CN&x='
urls = []
tiles = []
lines = []
# 循环构建url
for i in range(-row_half_range, row_half_range + 1, 1):
for j in range(-col_half_range, col_half_range + 1, 1):
url = base + (X + j).__str__() + '&y=' + (Y + i).__str__() + '&z=' + z.__str__()
urls.append(url)
# 循环获取瓦片
for item in urls:
tiles.append(getTile(item))
# 构建图像用于显示
final = np.zeros((row * 256, col * 256, 3), np.uint8)
# 按水平行进行瓦片拼接
for k in range(0, tiles.__len__(), col):
h_tuple = tuple(tiles[k:k + col])
line = np.hstack(h_tuple)
lines.append(line)
v_tuple = tuple(lines)
# 按竖直方向进行行拼接
final = np.vstack(v_tuple)
# 判断是否需要保存
if save_path != "no":
cv2.imwrite(save_path, final)
# 判断是否需要缩放
if scale_factor != 1:
final = cv2.resize(final, None, fx=scale_factor, fy=scale_factor, interpolation=cv2.INTER_AREA)
# 显示图像
cv2.imshow("final", final)
cv2.waitKey(0)
打开程序，依次按照提示输入经纬度、层数、行数、列数、缩放因子、输出路径，按回车便可显示对应范围瓦片。
但是如果使用很多次之后就会发现，会遇到HTTP Error 403: Forbidden。
这其实是Python爬虫时常碰到的问题。原因是Google地图服务器拒绝了我们的请求。
因为我们在短时间内直接使用Get获取大量数据，会被服务器认为是在对它进行攻击，所以就自动把电脑IP封了。
解决这个问题有两种方法。一是将我们的请求加以包装，变成浏览器请求模式，而不再是“赤裸裸”的请求。
二是使用代理IP。具体内容之后再继续学习。
本文作者原创,未经许可不得转载,谢谢配合
