遥感影像中的知识点
1 安装环境
有幸遇到一个机会,接触到遥感影像,将其中遇到的所有问题记录一下。
1.1 ubuntu py3 GDAL环境
安装GDAL库
apt-get install -y python3-tifffile libgdal-dev gdal-bin python3-gdal
【注意】这样就在ubuntu下安装好了Python3所需要的GDAL库
安装python3所需库
pip3 install imagecodecs-lite tifffile shapely matplotlib
pip3 install imagecodecs tifffile shapely matplotlib
1.2 win10 py3 GDAL环境
进gisinternals官网选择指定的MSVC x64版本,点击release-1911-x64-gdal-3-0-4-mapserver-7-4-3.zip进入页面,分别下载:
release-1911-x64-gdal-3-0-4-mapserver-7-4-3.zip
GDAL-3.0.4.win-amd64-py3.7.msi
将release-1911-x64-gdal-3-0-4-mapserver-7-4-3.zip解压到C:\JAVA下,然后配置环境变量到PATH:
C:\JAVA\release-1911-x64-gdal-3-0-4-mapserver-7-4-3\bin
再双击GDAL-3.0.4.win-amd64-py3.7.msi安装对应的python3.7库中去。
再验证是否安装成功:
from osgeo import ogr, osr, gdal
1.3 win10 ArcGIS环境
经发现,ArcGIS可以直接打开tiff格式的遥感影像数据(一般一张tiff是G级别),可观察tiff全图,安装方法可百度。
1.4 将mask写入shp
这部分可参考Python-GDAL教程:面矢量数据的写入-中级
1.5 TIFF的切割
这部分可参考Python 实现遥感影像波段组合的示例代码,我这部分实现的部分代码如下:
# coding=utf-8
# __author__=whaozl
import os
import numpy
from osgeo import gdal
from tqdm import tqdm
import shutil
"""
参考 https://www.jb51.net/article/166876.htm
"""
class TIFReader:
def __init__(self, filename):
self.filename = filename
self.proj = None
self.geotrans = None
self.data = None
# 读图像文件
def read_img(self):
ds = gdal.Open(self.filename) # 打开文件
self.im_width = ds.RasterXSize # 栅格矩阵的列数
self.im_height = ds.RasterYSize # 栅格矩阵的行数
self.geotrans = ds.GetGeoTransform() # 仿射矩阵
self.proj = ds.GetProjection() # 地图投影信息
self.data = ds.ReadAsArray(0, 0, self.im_width, self.im_height) # 将数据写成数组,对应栅格矩阵
del ds
# 写文件,以写成tif为例
def write_img(self, filename, im_proj, im_geotrans, im_data):
# gdal数据类型包括
# gdal.GDT_Byte,
# gdal .GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32,
# gdal.GDT_Float32, gdal.GDT_Float64
# 判断栅格数据的数据类型
if 'int8' in im_data.dtype.name:
datatype = gdal.GDT_Byte
elif 'int16' in im_data.dtype.name:
datatype = gdal.GDT_UInt16
else:
datatype = gdal.GDT_Float32
# 判读数组维数
if len(im_data.shape) == 3:
im_bands, im_height, im_width = im_data.shape
else:
im_bands, (im_height, im_width) = 1, im_data.shape
# 创建文件
driver = gdal.GetDriverByName("GTiff") # 数据类型必须有,因为要计算需要多大内存空间
ds = driver.Create(filename, im_width, im_height, im_bands, datatype)
ds.SetGeoTransform(im_geotrans) # 写入仿射变换参数
ds.SetProjection(im_proj) # 写入投影
if im_bands == 1:
ds.GetRasterBand(1).WriteArray(im_data) # 写入数组数据
else:
for i in range(im_bands):
ds.GetRasterBand(i + 1).WriteArray(im_data[i])
del ds
def split_tif(self, seg_root:str, size:int):
"""不重采样,直接分割数据"""
if os.path.exists(seg_root):
shutil.rmtree(seg_root)
else:
os.makedirs(seg_root)
channel, width, height = self.data.shape
print("(channel, width, height)=(%s, %s, %s)" % (channel, width, height))
for i in tqdm(range(width // size)): # 切割成256*256小图
for j in range(height // size):
cur_image = self.data[:, i * size:(i + 1) * size, j * size:(j + 1) * size]
self.write_img(os.path.join(seg_root, '{}_{}.tif'.format(i, j)), self.proj, self.geotrans, cur_image)
def split_tif_resample(self, seg_root:str, size:int, step:int):
"""重采样
:param seg_root : 分割后的数据保存目录
:param size : 分割后的tiff尺寸大小 这里约定size为1024
:param step : 分割的步长,这里约定步长为512"""
if os.path.exists(seg_root):
shutil.rmtree(seg_root)
else:
os.makedirs(seg_root)
channel, width, height = self.data.shape
print("(channel, width, height)=(%s, %s, %s)" % (channel, width, height))
width_count = width // step # 38588 // 512 = 75 (75.3671875)
height_count = height // step # 68067 // 512 = 132 (132.943359375)
for i in tqdm(range(width_count)):
for j in range(height_count):
width_begin = i * step
width_end = width_begin + size
if width_end > width:
width_begin = width - size
width_end = width
height_begin = j * step
height_end = j * step + size
if height_end > height:
height_begin = height - size
height_end = height
cur_image = self.data[:, width_begin: width_end, height_begin : height_end ]
self.write_img(os.path.join(seg_root, '{}_{}_x={}_y={}.tif'.format(i, j, width_begin, height_begin)), self.proj, self.geotrans, cur_image) # x,y是左上角顶点坐标
if __name__ == '__main__':
tiff_path = "/data/L19-prj.tif"
ob = TIFReader(tiff_path)
2 提取道路
这部分思路是,先将tiff切割为很小很小的比如1024x1024的图片,然后再用深度学习模型识别。
深度学习这部分可参考基于深度学习的实现影像地图道路提取CVPR2018的挑战赛任务中有提取道路。
这部分具体实现代码可见下载链接
3 提取水面
提取水面的思路很直接,就是先使用argis将tiff直接导出,然后发现通过阈值可以直接确定水面,具体可直接参考:
基于python3+opencv3遥感影像的湖泊边界提取,当然,若提供了其他颜色通道,也可以使用水体公式NDWI。
具这部分具体实现代码可见下载链接
Acknowledge
感谢 @杰 @康 两位以前朋友的技术上的咨询和请教。方可在短短几天内完成以上的机会尝试。以及感谢GitHub上一些朋友的问题及时回复。