美丽的滴水湖

  美丽的滴水湖坐落在上海的东南角,濒临东海,风景秀丽,安静舒适,是旅游、恋爱的绝佳去处。笔者有幸去过一回,对那儿的风土人情留下了深刻的印象,如果有机会,笔者还会多去几次!

  滴水湖是个神奇的地方,神奇之处在于它的外形是一个正圆形,这源于城市规划者对临港新城的美好设想。每次路过这个美丽的湖时,笔者总会想:这个湖到底多大呢?

  本文将会谈到如何如何得到滴水湖的水面面积。是手动测量?是地质勘测?No,No,No,我们还是借用我们熟悉的工具,那就是Python!什么,Python还能用来测量滴水湖的水面面积?你确定不是在逗我?别急,听笔者娓娓道来~

准备工作

  我们的准备工作很简单,打开百度地图,不需要考虑地图的比例尺,然后截两张图,一张截图(dishuihu.png)里面含有滴水湖的全部,一张截图(lingang_to_dishuihu.png)的一个角是临港大道地铁站,另一个角是滴水湖地铁站。如下图:

利用Python测量滴水湖的水面面积-LMLPHP

利用Python测量滴水湖的水面面积-LMLPHP

  dishuihu.png的大小为720*694, lingang_to_dishuihu.png的大小为444*451. 我们将要用到的工具是Python的图像处理工具模块cv2,它是OpenCV的Python接口。

滴水湖的直径

  首先我们需要将滴水湖的外形,即滴水湖所在的圆,检测出来,我们用的方法是霍夫变换(Hough Transform)圆检测。处理的Python代码如下:

import cv2

imagepath = 'F://dishuihu.png'
image = cv2.imread(imagepath, 3)
# 把图片转换为灰度模式
gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY) # 霍夫变换圆检测
circles= cv2.HoughCircles(gray, cv2.HOUGH_GRADIENT, 1, 100, param1=150, param2=50, minRadius=5, maxRadius=500)
# 输出霍夫变换圆检测的返回值
print('霍夫变换圆检测的信息:\n', circles)
# 输出检测到圆的个数
print('圆的个数为%d个.'%len(circles[0])) #根据检测到圆的信息,画出每一个圆
for circle in circles[0]:
#坐标行列
x=int(circle[0])
y=int(circle[1])
#半径
r=int(circle[2])
# 圆的基本信息
print('圆的中心为(%s, %s), 半径为%s.' % (x, y, r))
#在原图用指定颜色标记出圆的位置
image = cv2.circle(image, (x,y), r, (0,0,255), 1) #显示新图像
cv2.imshow('res', image)
cv2.waitKey(0)

输出的结果如下:

霍夫变换圆检测的信心:
[[[ 351.5 341.5 327.1000061]]]
圆的个数为1个.
圆的中心为(351, 341), 半径为327.

检测到的圆如下如所示:

利用Python测量滴水湖的水面面积-LMLPHP

事实上,这个圆基本与滴水湖外圈所在的公路的圆重合。因此,我们得到的这个圆是完全OK的。

  接下来,我们要测量这个圆的直径。

  那么怎么测呢?这时候,就要派上大将lingang_to_dishuihu.png了。这张图片能够帮助我们得到地图的比例尺!But how? 其实呢,很简单。临港大道地铁站到滴水湖地铁站是一条直线,在百度地图上,利用测距工具(在地图右上方的工具箱能找到)得到这两个地铁站的距离是2.6公里,当然,这只是近似距离,但够用了。有了这个距离,我们就能得到图片中的一个像素代表的实际距离,利用以下Python代码(接上面代码)得到:

from math import sqrt
imagepath = 'F://lingang_to_dishuihu.png'
image2 = cv2.imread(imagepath, 3)
# 把图片转换为灰度模式
gray2 = cv2.cvtColor(image2, cv2.COLOR_BGR2GRAY) # 图片的宽度与高度
w,h = gray2.shape
# 图片的对角线长度
diag = sqrt(w**2+h**2)
map_scale = 2600/diag
print('图片的一个像素单位,代表实际距离:%.2f米'%map_scale)

输出结果如下:

图片的一个像素单位,代表实际距离:4.11米

  在刚才的圆检测中,我们得到滴水湖所在的圆的半径为327.1个像素单位,因此,滴水湖的直径为327.124.11,约为2686.77米,那么,滴水湖的面积为566.96万平方米。额,这样算出来的会对吗?来看看百度百科的说法吧。

利用Python测量滴水湖的水面面积-LMLPHP

  这简直棒呆了!

  下一步工作,就是求解滴水湖的水面面积~

滴水湖的水面面积

  什么是滴水湖的水面面积?那就是湖的水面的面积,也就是图片中的蓝颜色部分。可是,水面的形状看上去如此不规则,该怎么计算呢?别担心,山人自有妙计。

  第一步,去掉环形公路的外部部分,将其颜色设置为白色,避免将环形公路外的水面加进来。处理的Python代码(接上面的代码)如下:

height, width = gray.shape
for i in range(width):
for j in range(height):
if sqrt((i-x)**2+(j-y)**2) >= r:
image[j,i] = [255, 255, 255] #显示新图像
cv2.imshow('res', image)
cv2.waitKey(0)

处理后的图片如下:

利用Python测量滴水湖的水面面积-LMLPHP

  第二步,绘制刚才图片的灰色直方图,为了求出图片的颜色阈值,以方便后续对图片做二值化处理。处理的代码(接上面代码)如下:

# 在灰度图片中,去掉环形公路外面的部分
height, width = gray.shape
for i in range(width):
for j in range(height):
if sqrt((i-x)**2+(j-y)**2) >= r:
gray[j,i] = 255 # 绘制灰度图片的颜色直方图
from matplotlib import pyplot as plt plt.hist(gray.ravel(), 256, [0,256])
plt.show()

得到颜色直方图分布如下:

由此可以分析得到,湖水的颜色应该在214附近。

  第三步,先对刚才的灰度图片进行中值滤波,然后对其进行二值化(黑白)处理,设定阈值为220,处理的代码(接上面代码)如下:

# 中值滤波
blur = cv2.medianBlur(gray, 17) # 二值化
ret,thresh = cv2.threshold(blur, 220, 255, cv2.THRESH_BINARY) cv2.imshow('res', thresh)
cv2.waitKey(0)

得到的二值化后的图片如下:

利用Python测量滴水湖的水面面积-LMLPHP

可以发现,这张图片中的黑色部分与滴水湖的水面基本上是重合的。因此,我们可以用图片的黑色部分来代替滴水湖的水面,为此,我们只需要计算黑色像素的个数,这样,我们就能求出滴水湖的水面面积了喂!

  第四步,计算滴水湖的水面面积。处理的代码如下(接上面代码):

import numpy as np

black_pixel_number = np.sum(thresh==0)
water_area = black_pixel_number*map_scale**2/10000
print('滴水湖的水面面积为%.2f万平方米。'%water_area)

输出的结果如下:

滴水湖的水面面积为429.24万平方米。

额,这样计算出来的水面面积会靠谱吗?下面,我们手动地来估算一下滴水湖的水面面积。

利用Python测量滴水湖的水面面积-LMLPHP

我们假设湖面的平均半径为1.2公里,再减去湖中三个小岛的面积,就是水面面积了。

利用Python测量滴水湖的水面面积-LMLPHP

因此,手动估算湖面的面积为\(\pi*1200*1200/10000-23.5-6-14\),约为408.89万平方米,当然,我们还漏算了正上方的那个小湖。

  总的来说,我们计算出来的滴水湖的水面面积,在允许的误差范围内,还是靠谱的。

问题来源

  笔者在乘车去阳澄湖的时候,查了一下阳澄湖的百度百科介绍,讲到了阳澄湖的蓄水量。之后又查了湖泊的蓄水量的计算方法,大概是利用微积分的思想求解的。当然,在地图上,是无法求出湖泊蓄水量的,但是,估算湖泊的水面面积应该还是可以的。因此,笔者找了之前去过的滴水湖作为实验样本,发现还是能估算出来的。

总结

  以上内容纯属自娱自乐,不能作为严谨地测量湖泊水面面积的方法。如有任何问题,欢迎大家交流~

  最后附上本文图片处理过程中完整的Python代码,供大家参考。

附录

  完整的Python代码如下:

import cv2

imagepath = 'F://dishuihu.png'
image = cv2.imread(imagepath, 3)
# 把图片转换为灰度模式
gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY) # 霍夫变换圆检测
circles= cv2.HoughCircles(gray, cv2.HOUGH_GRADIENT, 1, 100, param1=150, param2=50, minRadius=5, maxRadius=500)
# 输出霍夫变换圆检测的返回值
print('霍夫变换圆检测的信息:\n', circles)
# 输出检测到圆的个数
print('圆的个数为%d个.'%len(circles[0])) #根据检测到圆的信息,画出每一个圆
for circle in circles[0]:
#坐标行列
x=int(circle[0])
y=int(circle[1])
#半径
r=int(circle[2])
# 圆的基本信息
print('圆的中心为(%s, %s), 半径为%s.' % (x, y, r))
#在原图用指定颜色标记出圆的位置
image = cv2.circle(image, (x,y), r, (0,0,255), 1) #显示新图像
# cv2.imshow('res', image)
# cv2.waitKey(0) from math import sqrt
imagepath = 'F://lingang_to_dishuihu.png'
image2 = cv2.imread(imagepath, 3)
# 把图片转换为灰度模式
gray2 = cv2.cvtColor(image2, cv2.COLOR_BGR2GRAY) # 图片的宽度与高度
w,h = gray2.shape
# 图片的对角线长度
diag = sqrt(w**2+h**2)
map_scale = 2600/diag
print('图片的一个像素单位,代表实际距离:%.2f米'%map_scale) radius = map_scale*r
print('滴水湖的直径为:%.2f米。'%(2*radius)) import math area = math.pi*radius**2/10000
print('滴水湖的面积为:%.2f万平方米。'%area) height, width = gray.shape
for i in range(width):
for j in range(height):
if sqrt((i-x)**2+(j-y)**2) >= r:
image[j,i] = [255, 255, 255] #显示新图像
#cv2.imshow('res', image)
#cv2.waitKey(0) # 在灰度图片中,去掉环形公路外面的部分
height, width = gray.shape
for i in range(width):
for j in range(height):
if sqrt((i-x)**2+(j-y)**2) >= r:
gray[j,i] = 255 #cv2.imshow('res', gray)
#cv2.waitKey(0) # 绘制灰度图片的颜色直方图
#from matplotlib import pyplot as plt #plt.hist(gray.ravel(), 256, [0,256])
#plt.show() # 中值滤波
blur = cv2.medianBlur(gray, 17) # 二值化
ret,thresh = cv2.threshold(blur, 220, 255, cv2.THRESH_BINARY) #cv2.imshow('res', thresh)
#cv2.waitKey(0) import numpy as np black_pixel_number = np.sum(thresh==0)
water_area = black_pixel_number*map_scale**2/10000
print('滴水湖的水面面积为%.2f万平方米。'%water_area) # 手动估算湖面面积 approximately_water_area = (math.pi*1200**2/10000)-23.5-6-14
print('手动估算的滴水湖的水面面积为%.2f万平方米。'%approximately_water_area)
05-02 08:47