遥感图像定位 python python遥感库

遥感影像中的知识点

• 1 安装环境

• 1.1 ubuntu py3 GDAL环境
• 1.2 win10 py3 GDAL环境
• 1.3 win10 ArcGIS环境
• 1.4 将mask写入shp
• 1.5 TIFF的切割

• 2 提取道路
• 3 提取水面
• Acknowledge

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.zipGDAL-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上一些朋友的问题及时回复。