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

迷南。 2022-05-10 11:10 216阅读 0赞

### 美丽的滴水湖 ###

美丽的滴水湖坐落在上海的东南角，濒临东海，风景秀丽，安静舒适，是旅游、恋爱的绝佳去处。笔者有幸去过一回，对那儿的风土人情留下了深刻的印象，如果有机会，笔者还会多去几次！  
  滴水湖是个神奇的地方，神奇之处在于它的外形是一个正圆形，这源于城市规划者对临港新城的美好设想。每次路过这个美丽的湖时，笔者总会想：这个湖到底多大呢？  
  本文将会谈到如何如何得到滴水湖的水面面积。是手动测量？是地质勘测？No,No,No,我们还是借用我们熟悉的工具，那就是Python！什么，Python还能用来测量滴水湖的水面面积？你确定不是在逗我？别急，听笔者娓娓道来~

### 准备工作 ###

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

![dishuihu.png][]

![lingang\_to\_dishuihu.png][lingang_to_dishuihu.png]

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

### 滴水湖的直径 ###

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

输出的结果如下：

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

检测到的圆如下如所示：

![霍夫变换检测到的圆][1240]

事实上，这个圆基本与滴水湖外圈所在的公路的圆重合。因此，我们得到的这个圆是完全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.1\*2\*4.11,约为2686.77米，那么，滴水湖的面积为566.96万平方米。额，这样算出来的会对吗？来看看百度百科的说法吧。

![百度百科关于滴水湖的介绍][1240 1]

这简直棒呆了！  
  下一步工作，就是求解滴水湖的水面面积~

### 滴水湖的水面面积 ###

什么是滴水湖的水面面积？那就是湖的水面的面积，也就是图片中的蓝颜色部分。可是，水面的形状看上去如此不规则，该怎么计算呢？别担心，山人自有妙计。  
  第一步，去掉环形公路的外部部分，将其颜色设置为白色，避免将环形公路外的水面加进来。处理的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)

处理后的图片如下：

![去掉环形公路外面的部分][1240 2]

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

# 在灰度图片中，去掉环形公路外面的部分 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()

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

![灰度图片的颜色直方图][1240 3]

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

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

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

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

![二值化处理后的图片][1240 4]

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

第四步，计算滴水湖的水面面积。处理的代码如下（接上面代码）：

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万平方米。

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

![湖面的平均半径][1240 5]

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

![三个小岛的面积][1240 6]

因此，手动估算湖面的面积为 π ∗ 1200 ∗ 1200 / 10000 − 23.5 − 6 − 14 \\pi\*1200\*1200/10000-23.5-6-14 π∗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)

[dishuihu.png]: /images/20220510/6d590a8da42c4bd984b1c65427d0b8ed.png
[lingang_to_dishuihu.png]: /images/20220510/4ce79d0cb9404ee2a1d4c83c28a891c9.png
[1240]: /images/20220510/2a8e2d53f8664491b9e8c92d8410910d.png
[1240 1]: /images/20220510/10314cf26ffa47f5b58784dd2ea53b5c.png
[1240 2]: /images/20220510/51a1e46e3e5e4191bd2cdbf8e8997da8.png
[1240 3]: /images/20220510/f42ca0193c254eefac9b541ca038b31b.png
[1240 4]: /images/20220510/f99e89b8843c4f87b3cfaee2caf469d0.png
[1240 5]: /images/20220510/3e1c78b2b85941859d9d45ae13e994da.png
[1240 6]: /images/20220510/09e18f1839234161b24a1ab9867aa5b0.png