第五次作业20210516
写在前面
- 请打开简书链接审阅此份作业！第五次作业-简书 强烈建议用简书查看，对文字的渲染和排版比邮件强很多很多
- 本文是用markdown写的，发送邮件时使用了markdown here插件作为html格式发送，文中所有图像均使用外链嵌入。如果出现排版严重错误，看不到图片等情况，请立刻用邮件联系我。
- 文中截图为了保证完整性，尺寸较大，而且是在21：9的屏幕中截取，请左右拖动查看。
- 邮件正文已经包括了作业的全部内容。
第一题
1.改造淹没分析代码（B13346_08_03-flood），对srtm.asc数据进行淹没分析计算，要求：
（1）修改代码中的淹没种子点位置坐标参数，从像素坐标改为地理坐标进行计算（即加入通过地理坐标计算得到像素坐标的过程）；
（2）适当修改淹没分析参数，包括种子点位置和淹没高程，不同组合实验5次，并截图提交；
（3）淹没高程参数从种子点参数读出，可根据淹没种子点位置的高程增加一定的高度值求出。
源程序简单分析
"""
Crawls a terrain raster from a starting
point and "floods" everything at the same
or lower elevation by producing a mask
image of 1, 0 values.
"""
# http://git.io/v3fSg
import numpy as np
from linecache import getline
def floodFill(c, r, mask):
"""
Crawls a mask array containing
only 1 and 0 values from the
starting point (c=column,
r=row - a.k.a. x, y) and returns
an array with all 1 values
connected to the starting cell.
This algorithm performs a 4-way
check non-recursively.
"""
# cells already filled
filled = set()
# cells to fill
fill = set()
fill.add((c, r))
width = mask.shape[1]-1
height = mask.shape[0]-1
# Our output inundation array
flood = np.zeros_like(mask, dtype=np.int8)
# Loop through and modify the cells which
# need to be checked.
while fill:
# Grab a cell
x, y = fill.pop()
if y == height or x == width or x < 0 or y < 0:
# Don't fill
continue
if mask[y][x] == 1:
# Do fill
flood[y][x] = 1
filled.add((x, y))
# Check neighbors for 1 values
west = (x-1, y)
east = (x+1, y)
north = (x, y-1)
south = (x, y+1)
if west not in filled:
fill.add(west)
if east not in filled:
fill.add(east)
if north not in filled:
fill.add(north)
if south not in filled:
fill.add(south)
return flood
source = "terrain.asc"
target = "flood.asc"
print("Opening image...")
img = np.loadtxt(source, skiprows=6)
print("Image opened")
a = np.where(img < 70, 1, 0)
print("Image masked")
# Parse the headr using a loop and
# the built-in linecache module
hdr = [getline(source, i) for i in range(1, 7)]
values = [float(h.split(" ")[-1].strip()) for h in hdr]
cols, rows, lx, ly, cell, nd = values
xres = cell
yres = cell * -1
# Starting point for the
# flood inundation
sx = 2582
sy = 2057
print("Beginning flood fill")
fld = floodFill(sx, sy, a)
print("Finished Flood fill")
header = ""
for i in range(6):
header += hdr[i]
print("Saving grid")
# Open the output file, add the hdr, save the array
with open(target, "wb") as f:
f.write(bytes(header, 'UTF-8'))
np.savetxt(f, fld, fmt="%1i")
print("Done!")
这个程序非常清晰，简单，从上往下看主要分几步：
- 建立函数模型：定义函数floodFill进行具体的淹没计算
- 建立遮罩，设定洪水淹没的最高高度
- 定义淹没起始点（sx，sy），注意，这里使用的是图像坐标
- 进行具体计算
- 导出淹没后的模型
重新构建程序
在理解了源程序的结构之后，分析一下我们作业的要求，简单概括来说：
- 使用经纬度作为淹没起始点
- 多次修改淹没参数（起始点，遮罩高程）并得出多个结果
- 从淹没起始点的高程计算淹没高程参数
感觉每次分析前再去看淹没起始点高程很麻烦，我就直接把功能3作为默认功能封装在函数里，这样就比较方便。
程序功能介绍与代码分析
整个程序功能分为三个函数实现
import部分
import numpy as np
from linecache import getline
import gdal as gd
以下部分np就是numpy，gd就是gdal
1.地理坐标转换图像坐标函数geo2image（完成要求1）
def geo2image(filepath, x, y):
"""
根据GDAL的参数模型将给定的投影或地理坐标转为影像图上坐标（行列号）
:param filepath: 地理文件
:param x: 投影或地理坐标x，经度
:param y: 投影或地理坐标y，纬度
:return: 影坐标或地理坐标(x, y)对应的影像图上行列号(row, col)
"""
dataset = gd.Open(filepath)
trans = dataset.GetGeoTransform()
a = np.array([[trans[1], trans[2]], [trans[4], trans[5]]])
b = np.array([x - trans[0], y - trans[3]])
# 使用numpy的linalg.solve进行二元一次方程的求解
resultX, resultY = np.linalg.solve(a, b)
return int(resultX + 0.5), int(resultY + 0.5)
逻辑还是很清晰，用gdal的功能和简单的投影计算进行一个坐标的转换
我这里最后加了一句四舍五入的逻辑，因为毕竟像素点是整数，函数返回的也应是整数结果
2.淹没分析函数floodFill
def floodFill(c, r, mask):
"""
从起始点(c=column,r=row - a.k.a. x, y)
开始遍历只包含0和1的遮罩数组
然后返回和起始单元格相连并且值为1的数组
该算法会进行非递归的、4个方向上的检查
:param c: column
:param r: row
:param mask: 遮罩
:return: 淹没后图像
"""
# 已填充的单元格
filled = set()
# 待填充的单元格
fill = set()
fill.add((c, r))
width = mask.shape[1] - 1
height = mask.shape[0] - 1
# 输出的淹没数组
flood = np.zeros_like(mask, dtype=np.int8)
# 遍历和修改需要检查的单元格。当填充时获取一个单元格
while fill:
# 获取单元格
x, y = fill.pop()
if y == height or x == width or x < 0 or y < 0:
# 不填充
continue
if mask[y][x] == 1:
# 填充
flood[y][x] = 1
filled.add((x, y))
# 检查相邻单元格
west = (x - 1, y)
east = (x + 1, y)
north = (x, y - 1)
south = (x, y + 1)
if west not in filled:
fill.add(west)
if east not in filled:
fill.add(east)
if north not in filled:
fill.add(north)
if south not in filled:
fill.add(south)
return flood
这个也是非常的清晰，进行了一个非递归的遍历
这里我的实现和给出的示例代码是一样的，坐标系转换的逻辑用下一个函数实现
3.保存淹没后的asc文件到硬盘的saveFloodAsc（完成要求3）
def saveFloodAsc(source, target, lon, lat, above):
"""
保存淹没分析的图片，直接输入比淹没的位置高多少
:param source: 淹没底图，格式asc
:param target: 淹没成果图，格式asc
:param lon: 淹没起始点经度
:param lat: 淹没起始点纬度
:param above: 淹没的水最高比起始点高多少
:return: void
"""
print("打开图片中，可能耗时较长")
img = np.loadtxt(source, skiprows=6)
print("图片已载入")
# 测算淹没起始点图像坐标
sx, sy = geo2image(source, lon, lat)
height = img[sx, sy]
print("淹没起始点坐标：经度{lonF} 纬度{latF} 对应图像坐标（{sxF}，{syF}）起始点高度{h}".format(lonF=lon, latF=lat, sxF=sx, syF=sy, h=height))
# 遮罩高程，可以近似理解为，设定洪水的最高淹没高度
# 这里直接搞一个，最高比起始点高多少的算法，非常简单
a = np.where(img < height + above, 1, 0)
print("遮罩已生成")
# 使用linecache解析头部信息
hdr = [getline(source, i) for i in range(1, 7)]
values = [float(h.split(" ")[-1].strip()) for h in hdr]
print("开始淹没分析")
fld = floodFill(sx, sy, a)
print("淹没分析完成")
header = ""
for i in range(6):
header += hdr[i]
print("保存图像中")
# 打开输出文件并保存
with open(target, "wb") as f:
f.write(bytes(header, 'UTF-8'))
np.savetxt(f, fld, fmt="%1i")
print("全部完成")
程序中间加了个左边转换的逻辑，就是把起算点的坐标转成图像坐标，然后读出高程，再用高程+above作为遮罩的阈值，最后保存图像，非常简单直观
五个示例结果
先用gdal查看图片，选五个幸运点作为淹没的起始点
图 1
- 71.564563,36.776074
图 2 - 71.81531,36.35192
图 3 - 70.4150774,35.8929393
图 4 - 72.267072,39.188082
图 5 - 74.17957,37.22421
图 7
直接进行一个运算，使用代码
source = "srtm.asc"
saveFloodAsc(source, "flood1.asc", 71.564563, 36.776074, 50)
saveFloodAsc(source, "flood2.asc", 71.81531, 36.35192, 40)
saveFloodAsc(source, "flood3.asc", 70.4150774, 35.8929393, 30)
saveFloodAsc(source, "flood4.asc", 72.267072, 39.188082, 20)
saveFloodAsc(source, "flood5.asc", 74.17957, 37.22421, 10)
其实算的挺久的，图片有点大，一个几百M，要导入挺久，写到这想到可以先把图像写到内存里，这样应该快一点
但是也不会快太多，因为最慢的还是淹没分析那段的遍历，真的很慢
这个地方如果用MATLAB搞可以很方便的GPU并行计算，MATLAB只需要一个函数就能调用GPU资源，真的方便，python或者c++用cuda搞并行就有点麻烦。很想用MATLAB处理这个题。
结果截图
关于这个渲染方式，我设置了一下，就是把淹没区域搞成透明，黑色区域是没淹到的
-
第一个点，淹没点上的50M区域.其实感觉淹的地方不多，应该是这个点本身比较低
图 1 -
第二个点，淹40m，这直接都淹没了，全是透明区域
图 2 -
这就没淹多少，一大半地方都没被淹，起始点太低了吧，感觉这个dem图就是东高西低
图 3 -
这也没淹多少，淹没高度低了
图 5 -
这个淹了好多啊，虽然只设置了10m，还是起点高
图 6
附录：第一题的完整代码
"""
Crawls a terrain raster from a starting
point and "floods" everything at the same
or lower elevation by producing a mask
image of 1, 0 values.
"""
# http://git.io/v3fSg
import numpy as np
from linecache import getline
import gdal as gd
def geo2image(filepath, x, y):
"""
根据GDAL的参数模型将给定的投影或地理坐标转为影像图上坐标（行列号）
:param filepath: 地理文件
:param x: 投影或地理坐标x，经度
:param y: 投影或地理坐标y，纬度
:return: 影坐标或地理坐标(x, y)对应的影像图上行列号(row, col)
"""
dataset = gd.Open(filepath)
trans = dataset.GetGeoTransform()
a = np.array([[trans[1], trans[2]], [trans[4], trans[5]]])
b = np.array([x - trans[0], y - trans[3]])
# 使用numpy的linalg.solve进行二元一次方程的求解
resultX, resultY = np.linalg.solve(a, b)
return int(resultX + 0.5), int(resultY + 0.5)
def floodFill(c, r, mask):
"""
从起始点(c=column,r=row - a.k.a. x, y)
开始遍历只包含0和1的遮罩数组
然后返回和起始单元格相连并且值为1的数组
该算法会进行非递归的、4个方向上的检查
:param c: column
:param r: row
:param mask: 遮罩
:return: 淹没后图像
"""
# 已填充的单元格
filled = set()
# 待填充的单元格
fill = set()
fill.add((c, r))
width = mask.shape[1] - 1
height = mask.shape[0] - 1
# 输出的淹没数组
flood = np.zeros_like(mask, dtype=np.int8)
# 遍历和修改需要检查的单元格。当填充时获取一个单元格
while fill:
# 获取单元格
x, y = fill.pop()
if y == height or x == width or x < 0 or y < 0:
# 不填充
continue
if mask[y][x] == 1:
# 填充
flood[y][x] = 1
filled.add((x, y))
# 检查相邻单元格
west = (x - 1, y)
east = (x + 1, y)
north = (x, y - 1)
south = (x, y + 1)
if west not in filled:
fill.add(west)
if east not in filled:
fill.add(east)
if north not in filled:
fill.add(north)
if south not in filled:
fill.add(south)
return flood
def saveFloodAsc(source, target, lon, lat, above):
"""
保存淹没分析的图片，直接输入比淹没的位置高多少
:param source: 淹没底图，格式asc
:param target: 淹没成果图，格式asc
:param lon: 淹没起始点经度
:param lat: 淹没起始点纬度
:param above: 淹没的水最高比起始点高多少
:return: void
"""
print("打开图片中，可能耗时较长")
img = np.loadtxt(source, skiprows=6)
print("图片已载入")
# 测算淹没起始点图像坐标
sx, sy = geo2image(source, lon, lat)
height = img[sx, sy]
print("淹没起始点坐标：经度{lonF} 纬度{latF} 对应图像坐标（{sxF}，{syF}）起始点高度{h}".format(lonF=lon, latF=lat, sxF=sx, syF=sy, h=height))
# 遮罩高程，可以近似理解为，设定洪水的最高淹没高度
# 这里直接搞一个，最高比起始点高多少的算法，非常简单
a = np.where(img < height + above, 1, 0)
print("遮罩已生成")
# 使用linecache解析头部信息
hdr = [getline(source, i) for i in range(1, 7)]
values = [float(h.split(" ")[-1].strip()) for h in hdr]
print("开始淹没分析")
fld = floodFill(sx, sy, a)
print("淹没分析完成")
header = ""
for i in range(6):
header += hdr[i]
print("保存图像中")
# 打开输出文件并保存
with open(target, "wb") as f:
f.write(bytes(header, 'UTF-8'))
np.savetxt(f, fld, fmt="%1i")
print("全部完成")
# 以下为主程序
source = "srtm.asc"
saveFloodAsc(source, "flood1.asc", 71.564563, 36.776074, 50)
saveFloodAsc(source, "flood2.asc", 71.81531, 36.35192, 40)
saveFloodAsc(source, "flood3.asc", 70.4150774, 35.8929393, 30)
saveFloodAsc(source, "flood4.asc", 72.267072, 39.188082, 20)
saveFloodAsc(source, "flood5.asc", 74.17957, 37.22421, 10)
第二题
改造遥感数据分类计算程序（B13346_06_05-classify），对GF1.jpg进行分类计算，分别按5，10，15个类别进行分类计算，并适当修改颜色列表参数（可随机生成），将结果截图表达。
源程序分段分析
这个源程序第一句话就无法运行......应该是python版本问题，新版的gdal_array不在gdal里面，反而在osgeo包里，改成
from osgeo import gdal_array
才能运行。最可气的是...中文版教材里译者已经做出了订正，官方的GitHub都没有改...实在是...
这个程序核心是这个二层循环
for i in range(len(classes)):
mask = gdal_array.numpy.logical_and(start <= srcArr, srcArr <= classes[i])
for j in range(len(lut[i])):
rgb[j] = gdal_array.numpy.choose(mask, (rgb[j], lut[i][j]))
start = classes[i] + 1
就是根据lut数组内定的类别对图像分类
对源程序进行修改
实现要求：修改颜色列表
题目要求修改颜色参数列表，先实现一个随机生成num种分类的lut数组函数吧
def lut_random_maker(num):
"""
制作随机种类数的lut
:param num: 种类数
:return: lut数组
"""
lut = []
i = 0
while i < num + 1:
# 这里加一，是因为颜色查找表的记录数必须为len(classes)+1
i = i + 1
lut.append(np.random.randint(0, 255, 3))
return lut
这个颜色其实就是一个1*3数组，就直接用numpy这个功能生成，非常简单。这里注意循环的时候要写num+1，实际生成的lut数组应比种类数多1
实现要求：分5 10 15类各做一遍
因为这个代码重复三遍很麻烦，我就把程序写成了一个函数
def img_classify(src, num, lut):
"""
根据参数分类图像
:param src: 原图像
:param num: 分类的类别个数
:param lut: RGB元组
:return: void
"""
因为很长，这里不把全部贴出来，函数头就是这样。程序运行这样写
my_src = "GF1.jpg"
img_classify(my_src, 5, lut_random_maker(5))
img_classify(my_src, 10, lut_random_maker(10))
img_classify(my_src, 15, lut_random_maker(15))
就能完成功能、当然，后面我会贴出本题使用的完整程序
输出图像
图 1
效果感觉是正常的，升随机颜色边数挺大
附录：第二题的完整代码
"""遥感影像分类"""
# https://github.com/GeospatialPython/Learn/raw/master/thermal.zip
from osgeo import gdal_array
import numpy as np
def lut_random_maker(num):
"""
制作随机种类数的lut
:param num: 种类数
:return: lut数组
"""
lut = []
i = 0
while i < num + 1:
# 这里加一，是因为颜色查找表的记录数必须为len(classes)+1
i = i + 1
lut.append(np.random.randint(0, 255, 3))
return lut
def img_classify(src, num, lut):
"""
根据参数分类图像
:param src: 原图像
:param num: 分类的类别个数
:param lut: RGB元组
:return: void
"""
# 使用gdal库加载图片到numpy库
srcArr = gdal_array.LoadFile(src)
# 根据类别数目将直方图分割成20个颜色区间
classes = gdal_array.numpy.histogram(srcArr, bins=num)[1]
# 分类初始值
start = 1
# 创建一个RGB颜色的JPEG输出图片
rgb = gdal_array.numpy.zeros((3, srcArr.shape[0],
srcArr.shape[1],), gdal_array.numpy.float32)
# 处理所有类并声明颜色
for i in range(len(classes)):
mask = gdal_array.numpy.logical_and(start <= srcArr, srcArr <= classes[i])
for j in range(len(lut[i])):
rgb[j] = gdal_array.numpy.choose(mask, (rgb[j], lut[i][j]))
start = classes[i] + 1
# 保存图片
output = gdal_array.SaveArray(rgb.astype(gdal_array.numpy.uint8), str(num) + "分类后的" + src, format="JPEG")
output = None
my_src = "GF1.jpg"
img_classify(my_src, 5, lut_random_maker(5))
img_classify(my_src, 10, lut_random_maker(10))
img_classify(my_src, 15, lut_random_maker(15))
